NOVEL EFFICIENT INTEGRAL-BASED TECHNIQUES FOR CHARACTERIZATION OF PLANAR MICROWAVE STRUCTURES by Kazem Sabetfakhri A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy (Electrical Engineering) in The University of Michigan 1995 Doctoral Committee: Professor Linda P.B. Katehi, Chairperson Assistant Professor Kamal Sarabandi Professor Thomas B.A. Senior Professor Herbert G. Winful Associate Professor Andrew E. Yagle

K. Sabetfakhri 1995 All Rights Reserved

To my beloved wife Delband. To my sweet little Niki. And to my ever-inspiring parents.

ACKNOWLED GEMENTS The completion of this dissertation and my entire doctoral study at the University of Michigan would not have been possible without the support of many people. I would like to take this opportunity to express my sincere gratitude to all those who helped me in one way or another during these prolific years. First and foremost, I would like to thank my advisor Dr. Linda Katehi, whose insight, expertise, encouragement and understanding have been reflected throughout this work. I am also grateful for her continuous financial support and for providing me with the opportunity to attend numerous conferences to present my research work. I also thank the rest of my doctoral committee for their time and consideration. Special thanks go to Dr. Kamal Sarabandi, from whom I have constantly benefited through a great deal of advice and fruitful discussions. Many friends and colleagues at the Radiation Laboratory, either in the past or recently, have contributed to this work through endless technical or non-technical discussions. Special mention is deserved by Messrs. Jong-Gwan Yook and Thomas M. Weller, and Drs. Emilie van Deventer and Nihad I. Dib. In summary, all Rad Lab friends, too many to mention by name, created a very enjoyable environment for work and study. The major part of my study and research at the University of Michigan was sponsored by the Army Research Office under contract DAAL03-91-G-0116. The iii

continuous enthusiasm of this office in my research work is cordially appreciated. For the last year of my study, I received a Rackham Predoctoral Fellowship from the University of Michigan. This honor is also sincerely acknowledged. Last but certainly not least, I wish to thank my wife Delband, who has been a constant source of hope and inspiration during all these years. She proved the best companion one may ever long for throughout the long journey I started years ago. I also express my deepest gratitude to my parents Golnaz and Arjang Sabetfakhri, who sowed in me from the early years the seeds of aspiration to ascend to the highest levels in every aspect of life, especially in education and knowledge. Finally, I wish to thank my younger siblings Tannaz and Ali, who have always believed in me throughout my entire education. iv

TABLE OF CONTENTS DEDICATION................................ ii ACKNOWLEDGEMENTS................. iii LIST OF FIGURES............................... vii LIST OF TABLES................................ xi LIST OF APPENDICES............................ xii CHAPTER I. INTRODUCTION........................... 1 1.1 Motivation and Background................... 1 1.2 Integral Formulations and The Method of Moments...... 5 1.3 Objectives and Overview.................... 8 II. AN INTEGRAL TRANSFORM TECHNIQUE FOR PLANAR DIELECTRIC STRUCTURES............... 11 2.1 Introduction........................... 11 2.2 A Primer on Higher-Order Boundary Conditions....... 13 2.3 General Methodology...................... 18 2.4 The Modified Green's Functions................ 26 2.4.1 Dielectric Slab Waveguide............. 30 2.4.2 Dielectric Strip Waveguide.............. 31 2.4.3 Leaky Dielectric Waveguides............. 33 2.4.4 Coupled Dielectric Waveguides............ 36 2.5 Numerical Solution Using the Method of Moments...... 38 2.6 Numerical Results and Discussion............... 43 2.7 Concluding Remarks....................... 55 III. MULTIRESOLUTION EXPANSIONS FOR MOMENT METHOD SOLUTION OF ELECTROMAGNETIC PROBLEMS.... 57 v

3.1 Introduction.......................... 57 3.2 One-Dimensional Multiresolution Analysis...........60 3.3 The Fast Wavelet Algorithm................. 66 3.4 Two-Dimensional Multiresolution Analysis........... 68 3.5 Construction of Moment Method Expansions......... 71 3.6 Wavelets as Electromagnetic Radiators............. 78 IV. SPECTRAL-DOMAIN WAVELET ANALYSIS OF INTEGRATED PLANAR WAVEGUIDES............... 82 4.1 Introduction.......................... 82 4.2 Formulation of Integral Equations............... 84 4.3 Multiresolution Expansion of Transformed Currents.....91 4.4 The Moment Method Solution and Numerical Considerations 93 4.5 Numerical Results........................ 97 4.6 Concluding Remarks....................... 110 V. SPACE-DOMAIN WAVELET ANALYSIS OF THREE-DIMENSIONA PLANAR MICROWAVE STRUCTURES............ 112 5.1 Introduction.......................... 112 5.2 Preliminary Numerical Considerations............. 116 5.3 A Wavelet-Based Study of Dielectric Resonators....... 122 5.3.1 Integral Formulation................. 122 5.3.2 Numerical Procedure and Results.......... 126 5.4 A Wavelet-Based Study of Microstrip Patch Antennas.... 133 5.4.1 Integral Formulation................. 137 5.4.2 The Microstrip Patch as a Resonating Structure..142 5.4.3 The Microstrip Patch as a Radiating Element.. 145 5.5 Concluding Remarks....................... 153 VI. CONCLUSION............................. 156 6.1 Summary of Achievements................. 156 6.2 Recommendations for Future Work...............158 APPENDICES..................... 160 BIBLIOGRAPHY.............................. 184 vi

LIST OF FIGURES Figure 1.1 A typical integrated planar microwave structure............ 2 2.1 Geometry of interface between two different dielectric media.... 14 2.2 Geometry of a planar dielectric waveguide embedded in a general layered background structure...................... 19 2.3 Geometry of a dielectric strip waveguide................ 32 2.4 Contour deformation in the complex plane for leaky dielectric waveguides................................... 35 2.5 Geometry of an array of coplanar dielectric waveguides........ 37 2.6 A typical component of the transform-domain unknown J.(() and its inverse Fourier transform........................ 40 2.7 Orthogonal Hermite-Gauss functions.................. 42 2.8 Normalized propagation constants of dominant EI and E1f modes of a rectangular dielectric slab as a function of normalized frequency B..................................... 46 2.9 Convergence of the normalized propagation constants of dominant modes of a rectangular dielectric slab with respect to the number of basis functions.............................. 47 2.10 Normalized propagation constants of higher-order modes of a rectangular dielectric slab as a function of normalized frequency B... 48 2.11 Normalized propagation constant of the dominant modes of a highcontrast dielectric slab as a function of normalized frequency B... 49 vii

2.12 Normalized propagation constants of dominant symmetric and antisymmetric modes of a coupled dielectric slab waveguide as a function proximity ratio s/a........................... 50 2.13 Normalized propagation constants of dominant Efl and E1 modes of a non-leaky dielectric strip waveguide as a function of normalized strip thickness.............................. 52 2.14 Normalized propagation constants of the dominant Efl and Ez modes of a leaky dielectric strip waveguide as a function of strip thickness [in mm] at f = 30 GHz......................... 53 2.15 Attenuation constant of the leaky ET1 mode of the dielectric strip waveguide shown in Fig. 2.14 as a function of strip thickness [in mm] at f = 30 GHz....................... 54 3.1 Scaling function and mother wavelet of the Haar multiresolution analysis.................................... 61 3.2 The cubic spline Battle-Lemarie scaling function.......... 63 3.3 The cubic spline Battle-Lemarie mother wavelet........... 63 3.4 Fourier transform of the cubic spline Battle-Lemarie scaling function. 65 3.5 Fourier transform of the cubic spline Battle-Lemarie mother wavelet. 65 3.6 The wavelet localization lattice.................... 66 3.7 Signal flow graph for the 1-D FWA................. 68 3.8 Signal flow graph for the 2-D FWA................ 71 3.9 Typical locations of 2-D multiresoluiton basis functions...... 77 4.1 Geometry of an integrated planar waveguide structure...... 86 4.2 Normalized propagation constants of dominant Ef1 and El1 modes of a non-leaky single-strip dielectric waveguide with a single-layer substrate as a function of normalized strip thickness....... 99 4.3 The structure of moment matrix of Fig. 4.2 after applying a threshold of 1%.................................. 101 viii

4.4 Normalized propagation constant of the dominant mode of a singlestrip dielectric waveguide with a double-layer substrate as a function of frequency................................ 103 4.5 Normalized propagation constant of the dominant mode of a doublestrip dielectric waveguide with a single-layer substrate as a function of frequency............................... 104 4.6 Normalized propagation constant of a microstrip line printed over a double-layer substrate as a function of frequency............ 106 4.7 The structure of moment matrix of Fig. 4.6 after applying a threshold of 1%................................... 107 4.8 Normalized propagation constant of a hybrid microstrip-dielectric waveguide with a double-layer substrate as a function of frequency.. 108 4.9 The structure of moment matrix of Fig. 4.8 after applying a threshold of 1%................................... 109 5.1 Geometry of a rectangular dielectric resonator............ 123 5.2 Variation of E-field magnitude at the center of the resonator with the frequency of excitation field.................... 134 5.3 The structure of moment matrix of a dielectric resonator after applying a threshold of 1%...................... 135 5.4 Variation of resonant frequencies of dominant modes of a dielectric resonator as a function of the aspect ratio h/a............. 136 5.5 Geometry of a rectangular microstrip patch printed over a grounded substrate................................. 138 5.6 Variation of resonant frequency of a microstrip patch antenna as a function of patch width......................... 144 5.7 Patch current distribution along the direction of current....... 146 5.8 Patch current distribution normal to the direction of current..... 146 5.9 The structure of moment matrix of a microstrip patch after performing a threshold of 1%........................... 147 ix

5.10 Geometry of a coaxial-probe-fed microstrip patch antenna..... 149 5.11 Variation of patch input impedance as a function of frequency. Increment: 5MHz (increasing clockwise)................. 152 5.12 Variation of patch input impedance as a function of frequency. Increment: 10MHz (increasing clockwise)................. 154 A.1 Geometry of a double-layer grounded substrate............ 168 B.1 Branch cut for the logarithmic function in the extraction of surface wave poles................................. 172 x

LIST OF TABLES Table 2.1 A comparison between the integral transform technique and other methods in view of numerical efficiency................ 44 4.1 Effect of thresholding on the moment matrix............. 100 5.1 Resonant frequencies of different modes of a rectangular dielectric resonator................................. 132 C.1 Characteristic coefficients of the Battle-Lemarie multiresolution analysis.................................... 177 xi

LIST OF APPENDICES Appendix A. DYADIC GREEN'S FUNCTIONS FOR PLANAR LAYERED GEOM ETRIES................................ 161 A.1 The Unbounded Space...................... 163 A.2 Unbounded Half-space...................... 164 A.3 Single-layer Conductor-backed Substrate............ 165 A.4 Double-layer Conductor-backed Substrate........... 167 B. METHOD OF EXTRACTION OF SINGULARITIES......... 170 B.1 Extraction of Branch Point Singularities............ 171 B.2 Extraction of Surface Wave Poles................ 171 C. THE BATTLE-LEMARIE MULTIRESOLUTION ANALYSIS.... 174 D. METHOD OF STATIONARY PHASE................. 178 E. THE BI-CONJUGATE GRADIENT METHOD............ 181 xii

CHAPTER I INTRODUCTION 1.1 Motivation and Background The rapid growth of microwave and millimeter-wave integrated circuit technologies in recent years has initiated a wide range of new planar passive devices with high levels of complexity. These structures are composed of various planar metallic and dielectric components, which are fabricated collectively on the same substrate and are subsequently integrated into a larger circuit containing active devices. Fig. 1.1 shows a typical integrated planar geometry. In microwave and lower millimeterwave applications, planar circuits are mostly designed in the form of an assembly of microstrip-type, coplanar-waveguide-type or other planar metallized lines and segments, which are printed on a dielectric substrate through an appropriate fabrication process such as lithography. As the operating frequency increases into the higher millimeter-wave and submillimeter-wave regions, dielectric components are needed to substitute their metallic counterparts to avoid substantial conductor losses which would otherwise deteriorate the device performance in these frequency regions. Planar dielectric structures are also of central importance in the field of integrated optics. The electromagnetic characteristics of integrated planar structures often originate from complex field distributions, which do not render themselves to the simple ap1

2::.::::/::..,. -......:::. / ~~~~~~~~';":::0 ad /::::::~:: ~,~.] /~~~~~~~~~~ s _M F. E | _ i _ Or I:I..'.' / ~~~~~~~~~~~~~~~i: r::if_ l l -iii-ro ii:il:;iiiii /.<...;................:. —. -_............................. Figure LI:: A typical integrated planar microwave structure.

3 proximate models available for conventional microwave devices. The study of complex dielectric structures requires a rigorous full-wave analysis with a high degree of accuracy and versatility. Due to the lack of exact analytical solutions, this type of problems demand a robust numerical treatment, which does not suffer from unrealistic approximations. As the complexity of the geometry under study slightly increases, the approximate methods usually are not even able to provide a crude estimate of the solution. On the other hand, due to the complex nature of the underlying physical phenomena, the numerical modeling of integrated structures can easily turn into a time-consuming and expensive computational task with increasing the size and complexity of the problem. Thus, the availability of rigorous numerical techniques which rely upon modest computing resources without sacrificing the accuracy is a key factor in the efficient design of planar devices. During the past two decades, various approximate and potentially exact methods for the study of microwave circuit problems have been proposed which can adequately treat a sizable variety of two-dimensional geometries with a moderate degree of complexity [1-10]. Reference [11] presents a comprehensive survey of the available numerical techniques for microwave and millimeter-wave structures. Of this diverse repertory, however, the majority are limited in scope and applicability to 2-D or extremely simple 3-D problems, and very few have the capability of handling large-scale problems [12, 13, 14]. In particular, the current trend toward CAD-oriented modeling calls for flexible methods of a general nature, which can easily adapt to a wide variety of geometries without major modifications. The general-purpose full-wave techniques are based on the discretization of the geometry under study into an appropriate mesh, whereupon a rigorous numerical scheme for the solution of the related boundary-value problem is implemented. Upon discretization, the original formulation is normally re

4 duced to the solution of a linear system of algebraic equations. The accuracy of these techniques largely depends on the fineness of the discretization mesh. Hence. in large and complex problems, especially those involving three-dimensional geometries, to achieve a high degree of accuracy would require a very fine meshing, which in turn introduces a very large number of unknowns. The numerical processing of such big systems can easily exceed the capability of simple workstations and demands largescale computing resources with huge memory space and fast processors. It might be possible to manage a powerful computer to solve a big and time-consuming problem for validation or visualization purposes, but this approach certainly does not lend itself to the microwave design procedures, which usually involve iterative computations. Another major difficulty arising from large systems is the condition number of underlying matrices, which can severely deteriorate with the increasing size of matrices and may cause serious numerical problems. The computational difficulties due to the largeness of systems become even more pronounced in the case of open structures, which are of primary interest in the present work. Unlike shielded structures where the numerical modeling is confined to the interior of the enclosing box, for open geometries, the entire unbounded space must be modeled in an accurate and realistic manner. To this end, techniques like finite element and finite difference methods discretize a certain portion of the space surrounding the device under study, and then resort to some type of mesh termination by incorporating an appropriate absorbing boundary condition or placing an artificial absorber around the geometry [15, 16, 17]. In addition to the resulting major increase in the mesh size and hence the intensity of the computational work, the proper implementation of the mesh termination is another critical issue, which can dramatically complicate the numerical modeling of planar microwave devices.

5 By contrast, another class of full-wave techniques, which are based on the integral formulation of the related boundary value problem, are naturally appropriate for the treatment of open structures. In such techniques, only finite regions of the geometry where there is a substantial wave activity, and not the entire space, are discretized. In the following, we reflect upon this class of techniques in more detail. 1.2 Integral Formulations and The Method of Moments The integral equation technique is one of the oldest rigorous approaches for the numerical solution of electromagnetic problems [18]. In this technique, an integral equation is derived for an unknown field or current, which is usually sustained in a finite region of the geometry under study. The fields are postulated to be radiated by actual or equivalent sources. The effect of the entire background geometry including the unbounded sections are incorporated into the Green's function of the problem. This function is indeed the solution of the associated boundary value problem due to an elementary point source. Therefore, for the numerical solution of the integral equation, only the finite regions containing the radiating sources need be discretized. This amounts to a significant reduction in the number of unknowns and hence a better numerical efficiency. An attractive feature of integral formulations is that they provide some physical insight into the problem under investigation. This is particularly useful in the study of some special phenomena that are characteristic of open structures such as space and surface wave leakage effects. Pure numerical techniques like those mentioned earlier are naturally deprived of this privilege. However, a major limitation of the integral equation technique is the availability of analytical closed forms or efficient computation schemes for the relevant Green's functions. Yet,

6 this technique provides a potentially exact full-wave analysis of problems with welldocumented Green's functions. Integral-based techniques have been extensively used in the numerical modeling of planar passive microwave structures [5,19-23]. In spite of all the advantages listed above, the application of integral-based techniques has been limited mostly to small-scale or medium-scale electromagnetic problems. For large and complex structures, numerically intensive differential-based techniques such as finite element and finite difference methods have traditionally been preferred. This scenario stems for the most part from the difficulties associated with the numerical implementation of integral formulations for large-scale problems. The numerical solution of the integral equations encountered in electromagnetic problems is usually achieved through the use of method of moments (MoM) [24]. In this method, the unknown field or current is expanded in a set of discrete basis functions. The resulting discretized equations are then tested using a proper scheme such as Galerkin's technique to obtain a system of linear algebraic equations involving the unknown field or current amplitude coefficients. The numerical solution of this linear system leads to determination of the field profile in the entire structure. The characteristic parameters of the device under study are then computed from the field or current distributions. For small- to medium-scale problems with reasonable matrix sizes, this method offers a straightforward and efficient numerical solution. The difficulty arises when the complexity of the geometry and subsequently the number of unknowns increase resulting in very large matrices. The bottleneck is that all conventional types of basis functions traditionally used in the method of moments generate full moment matrices. This is due to the fact that these basis functions interact with each other very strongly and are basically good radiators. Moreover, even though the elements of moment matrices drop in magnitude as one goes away from the diagonal,

the bad condition numbers do not allow for any sort of thresholding. In large-scale problems such as complex 3-D geometries, the storage and manipulation of large, densely populated matrices turns into a formidable or even impossible computational task, which easily rules out the practicality of MoM-based approaches in competition with the sparse differential-based techniques. In the beginning of the present decade, it was discovered that the wavelet expansion of certain types of integral operators generates highly sparse linear systems [25]. First introduced in the closing years of the last decade, orthonormal wavelet expansions soon began to find a variety of interesting applications in such diverse fields as functional analysis, numerical analysis and signal processing [26, 27, 28, 30, 29, 31]. During the past two years, the concepts of multiresolution analysis theory, which is the mathematical foundation of wavelet expansions, were successfully applied to the numerical solution of a number of electromagnetic problems. It was reported by several researchers including the present author that the use of wavelet expansions in conjunction with the method of moments leads to highly sparse moment matrices [32, 33, 22, 34]. These exciting findings opened new horizons for the application of integral-based techniques to large-scale electromagnetic problems. It was established in effect that the traditional drawback of the method of moments in view of fullness of moment matrices is not an inherent feature of this technique, but it is a direct consequence of the type of basis functions used for the discretization of the underlying integral operators.

8 1.3 Objectives and Overview From the discussion of the preceding section, it can be concluded that thes t p applicability of the integral equation technique for the study of complex planar structures is contingent upon developing schemes or tools for effective reduction of the size of moment matrices. This goal may be achieved in two different ma tways: (1) by formulating novel integral equations which naturally involve smaller numbers of unknowns for a given required accuracy of the solution, or (2) by developing new moment method schemes which do not suffer from the computational deficiencies mentioned earlier. This work explores both of these two possibilities. The present research effort has produced two major outcomes: 1. A new efficient integral formulation, named the integral transform technique, has been developed which by using higher-order boundary conditions effectively reduces the dimensionality of the relevant boundary-value problem. In this technique, by introducing an appropriate integral transform, equivalent oneor two-dimensional transform-domain integral equations are derived for twoor three-dimensional electromagnetic problems, respectively. This reduction in dimensionality results in a significant reduction in the number of unknowns involved in the numerical solution of the modified integral equations. 2. A novel moment method formulation has been developed which exploits multiresolution expansions constructed based on the newly developed theory of orthonormal wavelets. It has been demonstrated that this approach leads to the generation of highly sparse moment matrices for two- and three-dimensional electromagnetic problems. Although the reported numerical results are mainly concerned with planar microwave structures, the basic methodology is presented

9 in a very general context, which can easily be applied to other types of electromagnetic problems or even other computational disciplines. The present work is one of the pioneering research efforts in the application of wavelet theory to computational electromagnetics. In particular, the results concerning three-dimensional structures are being reported for the first time to the best knowledge of the author. The present dissertation has been organized in six chapters and five appendices. In chapter 2, the basic methodology of the integral transform technique for planar dielectric structures is presented. Application of this technique to several two-dimensional dielectric waveguides is illustrated in detail and its numerical implementation through the method of moments is discussed. This chapter is concluded by emphasizing some practical limitations of the integral transform technique, which in reality led us to the second stage of this research effort concerning wavelet-based formulations. Chapter 3 is devoted to a review of multiresolution analysis theory, where the mathematical foundation of wavelet expansions is established. Some of the basic and useful properties of this type of expansions that will be used extensively in later formulations are explored in this chapter. The last two sections of chapter 3 discuss special considerations for the application of multiresolution expansions to electromagnetic problems. Chapter 4 presents a spectral-domain wavelet-based MoM formulation for integrated planar waveguides. Extensive numerical results are presented for a variety of practical waveguide configurations, and comparison is made with other available techniques for the validation of the results. In chapter 5, we develop a space-domain wavelet-based formulation for three-dimensional microwave structures. This formulation takes advantage of a very efficient numerical tool called the fast wavelet algorithm. Due to the large sizes of the problems under study, it is necessary to resort to special sparse

10 matrix techniques for the treatment of the resulting linear systems. It is shown that the use of the bi-conjugate gradient method for this purpose amounts to a dramatic improvement in the computational efficiency. A 3-D dielectric structure consisting of a rectangular dielectric resonator embedded in the free space and a 3-D metallic structure consisting of a rectangular patch antenna printed over an infinite grounded substrate are studied in detail. Finally, chapter 6 presents a brief summary and conclusion of the entire work, and gives recommendations for future work. Five appendices describe some of the mathematical procedures used in various chapters of this dissertation.

CHAPTER II AN INTEGRAL TRANSFORM TECHNIQUE FOR PLANAR DIELECTRIC STRUCTURES 2.1 Introduction Millimeter- and submillimeter-wave integrated circuits are usually made up of several metallic and dielectric planar sections embedded in a layered substrate structure. Depending on the complexity of the geometry, the full-wave analysis of these structures may become a difficult computational task. It should be noted that, unlike metallic strips which are modeled using planar conduction currents, dielectric sections are represented by equivalent volume polarization currents, whose dimensionality is higher than planar currents. In other words, the volume currents require a discretization along the normal to the substrate layers. In the past, some attempts have been made to derive surface current formulations for dielectric structures. Some of these works are based on equivalence theorems, and make use of equivalent electric and magnetic surface currents defined over an appropriate closed surface [35]. In a recent work, higher-order boundary conditions were employed to formulate an equivalent planar integral equation for shielded dielectric waveguides [36, 37]. In this formulation, the volume polarization current is integrated along the normal direction and the resulting average, which now has 11

12 the dimensionality of a planar current, is placed on an appropriate plane inside the region occupied by the dielectric material. This method, which is approximate by nature, provides good results for shielded step-index dielectric waveguides exciting LSE modes, but its application to open structures is limited to very thin dielectric layers. Within this limit, the first few low-order boundary conditions suffice and the inclusion of higher-order ones not only does not improve the results, but introduces spurious non-physical solutions [38]. In this chapter, a two-dimensional integral transform technique is developed for the analysis of open, or closed, planar dielectric waveguides of any size and medium parameters with a piecewise homogeneous layered substrate geometry. This method allows for an arbitrary refractive index profile within the guiding region. Thus, multilayered dielectric strips or embedded channels can easily be accommodated by regarding a piecewise-constant index profile. Although open geometries are of primary interest in this work, shielded structures can also be treated by using the relevant Green's functions. The integral transform technique provides a potentially exact formulation of the boundary value problem. First a transform domain is defined by introducing an appropriate integral transform in the direction of the substrate layers. The integral transform is chosen in accordance with the spectral representation of the Green's function of the substrate structure. Then, using the concept of the ceequivalent polarization current as a source in the guiding region, new transform-domain vector unknowns are introduced as integral transforms of this volume current with respect to kernels which are specified by the substrate Green's function. The fields in the source region are expanded along the normal to the substrate layers in a Taylor series at the planar boundaries, where the natural and higher-order boundary conditions are applied.

13 This leads to the formulation of second-kind Fredholm integral equations of reduced dimensionality for the transform-domain unknowns. Thus, two- or three-dimensional boundary-value problems are reduced to one- or two-dimensional transform-domain integral equations, respectively. For the numerical solution of the reduced integral equations, the method of moments is employed, whereupon the transform-domain unknowns are expanded in a complete set of orthonormal entire-domain basis functions. It is shown that the use of orthogonal Hermite-Gauss functions as the expansion basis proves to be numerically very efficient, and retaining only a few expansion terms leads to very satisfactory results. In the following, first the general methodology of the two-dimensional integral transform technique is presented in a form that can easily be extended to threedimensional geometries. Then, the implementation of the method of moments and the choice of basis functions are contemplated. For the validation of the technique, various planar dielectric waveguides are treated, and the numerical results are compared with other methods. 2.2 A Primer on Higher-Order Boundary Conditions In the methodology which we will develop in this chapter, in order to reduce the dimensionality of integral equations, the electric fields inside dielectric regions are expanded in Taylor series. Integration of such power series can be carried out analytically. Taylor expansions, however, involve field derivatives at the boundaries between dielectric regions and their background media. These derivatives must be related to the fields (and their derivatives) in the exterior regions, where the Green's functions are evaluated. To this end, higher-order boundary conditions are enforced

14 Z Figure 2.1: Geometry of interface between two different dielectric media. in conjunction with the Taylor expansions of the fields. Before proceeding to the formulation of reduced integral equations, in this section we present a brief review of these boundary conditions, which relate the electric field and its derivatives across the interface between two dielectric media. Suppose that the plane z = 0 is the interface between dielectric region 1 occupying the half-plane z < 0 with relative permittivity cr1 and dielectric region 2 occupying the half-plane z > 0 with relative permittivity Cr2, as shown in Fig. 2.1. The natural (zeroth-order) boundary conditions at z = 0 are stated as follows: 1 0 0 E(2)(x,y,O +)= 0 0 E()(, y, O-), C\ r2 (2.1) where E(i)(r), i = 1,2, are the electric fields in regions 1 and 2, respectively, and r is the position vector. Note that equation (2.1) holds even if the two dielectric media are inhomogeneous, and in this case, one should use the values of xr and Er2 at z = 0. This equation can be differentiated n times with respect to variables x and y to relate the nth-order field derivatives across the interface along these directions.

15 From Maxwell's equations, one can easily derive the first-order boundary conditions along the normal to the interface in the following form: Az O -i 9 TZ Oz a 0 0 -Z. (2.2) where S, = 1 - L is the relative index contrast. To derive higher-order boundary conditions, use is made of the wave equation in the two dielectric media, which is written in the following form: 2+ 2 2 +,lk (2) = 0, z < (a2 d22 + 2 + k2 E r 2) 0, z > 0. (2.3) +X2 8y2 0Z2 To relate the nth-order z-derivatives of the fields in the two regions, equations (2.3) are differentiated (n -2) times with respect to z and then subtracted from each other. Introducing the following notations: Sl = z ((E2) -E')), =x,= S: = a(E(2)-E ))X), Sn Zn (DZ Dz ) R - E v =x, v y, and the differential operator: 2 2 2 (2 R = + E2+k~, (2.5) n azn d 2

16 the following recurrence relation is obtained: S =- DS2-(E2-Er)k R_2, X n > 2, (2.6) with SO and S' being known from equations (2.1)-(2.2) (see [37]). In definitions (2.4), D(i)(r), i = 1,2, denote the electric displacement vectors in the two regions, respectively. The solution of equation (2.6) through a recursive procedure leads to the nth-order boundary condition at z = 0, which can be expressed in the following dyadic form: an 2 E(2)(x, y, 0+) Ln. E()(x, y, 0-), (2.7) where Ln is a dyadic differential operator. If we assume that the index profiles of the two dielectric media are constant or very slowly-varying functions in the immediate vicinity of the interface, then Ln is given by Tn 0 -OnQn x Ln= O Tn -OnQny l 0 0 Tn+OnQn a (2.8) where the differential operators Tn and Qn are given by: an [n/2] n-21 Znn aZn-21 = - ko2(Er2 - crl) E (-D)-' - — =1 Qn = (-Z)[n/2], (2.9) with On = 0 for even n, n, = 1 for odd n, and [n/2] denoting the largest integer less than or equal to n/2. For index profiles with rapid variations near the interface, Ln will of course contain derivatives of Er:.

17 For fields of exponential form, the dyadic operator of (2.8) reduces to a very simple form. To this end, assume that the electric fields in the two regions have the following spectral plane wave forms in the vicinity of the interface: E()(r) = E-i)e -jki. i = 1,2, (2.10) where k2 = k' + k2 + k, = eik2. Then, we have D = k22, the differential operator Tn becomes k2 [n/2] Tn = (-jkzl)n [1 k2 (Or2 -r1) k21, L. ZI, 1=0 and (2.8) reduces to the following form: (-jkl)n"-2[/2] 0 jkxsOn L= (-k )[n/2] 0 (-jkzl)n-2[n/2 jko,9a' 0 0 - (-jkzl,)n-2[n/2] _ jkz, nS (2.11) Separating the even and odd orders, and comparing with equations (2.1)-(2.2), it is found that L2m = (-kz2)mL 0 L2m+i = (-kz2)m L1. (2.12) In other words, the even- and odd-order boundary conditions are proportional to the zeroth- and first-order conditions, respectively.

18 2.3 General Methodology In this section, the two-dimensional integral transform technique is developed for a general class of planar dielectric waveguides. Fig. 2.2 shows the geometry of a y-directed dielectric waveguide which consists of a finite planar dielectric region of relative permittivity Crg(r), embedded in an infinite, piecewise homogeneous, layered, background structure. The relative permittivity of the substrate layer immediately surrounding the guiding region is denoted by eb. In the geometry of Fig. 2.2, y = 0 is the waveguide transverse plane, and the domains of the guiding region along the xand z-axes in this plane are denoted by D: and PD,, respectively. To include losses, the permittivities are assumed to be complex in general. The background structure can be any layered substrate configuration with or without conducting ground planes, or it can be an enclosed partially-filled waveguide geometry. We begin with the definition of a generalized integral transform. Given a function f (x), its integral transform with respect to the kernel g({x), abbreviated g-transform, is defined as f() - T [f(x) ] - j f (x)g(x)dx, (2.13) with the inverse transform f(x) = T'1 [f(s)] = K, J/ f(~)g*(x)d, (2.14) where Kg is a constant depending on the kernel g, Coo and C~ are appropriate contours, g* denotes the complex conjugate of g, and the integrals are assumed to be convergent [39]. The appropriate g-transform is chosen in accordance with the geometry of the problem and the Green's function of the background structure. For example, two- or three-dimensional planar geometries with piecewise rectangular cross

19 Erb - y x Figure 2.2: Geometry of a planar dielectric waveguide embedded in a general layered background structure. background structure.

20 sections prompt one- or two-dimensional Fourier transforms, respectively, while the suitable transform for a three-dimensional planar structure with circular cylindrical symmetry is a generalized Hankel transform. The designated g-transform must have a convolution property in the form of: f(x')f2(x- x')dx' I& c fl( )f2(I), (2.15) coo with a Parseval relation in the following manner: f (x) f2(x)dx = Kp, f(() f()d, (2.16) 00oQ. 00 where Kc and Kp are appropriate constants. Note that depending on the kernel g, the integrals of (2.16) may require a proper weight function, which is suppressed here without loss of generality. Adopting a harmonic time variation of the form exp(jwt), the electric and magnetic fields are expressed as E(r) = E(x,z)e-jP, H(r) = H(x,z)e-jy, (2.17) where a propagating mode of propagation constant /3 traveling along the positive yaxis has been assumed. Then, Maxwell's curl equations can be written in the following form: V x E(r) = -jkoZoH(r), V x H(r) = Jp(r) + jebkOYOE(r), (2.18) where V is the gradient operator, ko and Zo = 1/Yo are the free-space propagation constant and intrinsic impedance, respectively, and Jp(r) is an equivalent volume

21 polarization current defined as jkoYo [cg(r)- Cb]E(r), rE V, Jp(r) - (2.19) 0, elsewhere, with V being the volume of the dielectric guiding region. From equations (2.18), one can derive the wave equation for the electric field in the following way: V x V x E(r) - erbk2 E(r) = -jkoZo Jp(r) (2.20) along with the boundary conditions of the background structure. Using the dyadic form of Green's identities, it can be shown that the solution of the differential equation (2.20) in terms of the source function Jp(r) is given by E(r)= -jkoZOJJ Ge(r I r'). Jp (r') dv + {[~Vx,(,)] [An'x G]-[ h X E(r')]. [' x Ge]})ds', (2.21) where Ge(r I r') denotes the electric-field dyadic Green's function of the background -= T structure, G, is its dyadic transpose, V' is the gradient operator with respect to the i primed coordinates, and So is the boundary surface with in' being its outward unit normal [40]. In the case of an open geometry, So is chosen to be a circular cylinder of infinite radius surrounding the waveguide, whereupon the radiation condition forces the surface integral of (2.21) to become zero. For a shielded dielectric waveguide, So is chosen to coincide withthe enclosing metallic structure, where a Dirichlet condition again makes this integral vanish.

22 Since the waveguide geometry is of infinite extent along the y-axis, equation (2.21) can be integrated analytically with respect to y. Then, taking a g-transform of this equation with respect to the spatial variable x, i.e. in the direction of the substrate layers, and in view of property (2.15) of the g-transform, the following transformdomain equation is obtained: E(, D, z) = -jkoZo JD Ge(, Z, z'). Jp),, z)dz', (2.22) z where Ge( 3, z \ z') denotes the g-transform of the two-dimensional dyadic Green's function, with ( being the transform variable. The dyadic Green's function of background structures with piecewise homogeneous layered substrate geometries can be found by introducing Sommerfeld potentials or other choices of vector potentials [11]. A general plane-wave spectral eigenfunction expansion is presented in [41] for the Green's function of a planar multilayered medium with an arbitrary number of layers. From an inspection of these Green's functions, it is seen that its g-transform in the p-th substrate layer outside the source region can be expressed in the following general form: MP =p Ge(~,,z I Z')- 4 (, (, (,pz)hj(bz'), (2.23) j=1 where (p and (b are the complex eigenvalues along the z-axis in the p-th layer and the layer enclosing the guiding region, respectively, and satisfy the consistency conditions in these layers: C2 =_ 2 + P2 _ er, q = p, b. (2.24) Equation (2.23) represents a decomposition of the transform-domain dyadic Green's function into spectral terms which are separable with respect to the observation and source coordinates z and z' normal to the substrate layers. Note that within the

23 enclosing layer itself, this equation is valid only outside the guiding region, i.e., for z < TDz. The number of spectral terms, Mp, may vary from layer to layer, but in the most general case, Mp equals 2, complying with the primary and secondary source terms in the Sommerfeld representation. 4j and hj are exponential or hyperbolic functions of z and z', respectively, with real or complex arguments. The next section provides explicit examples of such a spectral decomposition. For the sake of simplicity, we drop the scripts p and b hereinafter. At this point, we define new transform-domain vector unknowns J1j in the following manner: Jej(,P)= - Jp(f3,z,)hi((4z)dz, j= 1,...,M, (2.25) which are indeed two-dimensional integral transforms of the volume polarization current Jp,(r), including a g-transform in the x direction and another generalized integral transform in the y direction with respect to the kernels hj((z). To ensure the convergence of the integrals in equation (2.25), special care must be taken in defining the functions hj of the spectral decomposition (2.23). Outside the guiding region, equation (2.22) now becomes M E(, 3, z) = -jkoZo E j( / ),,z) Joj(z, /), z i D,. (2.26) j=1 The fields in the guiding region and its enclosing or neighboring layers are related through the boundary conditions at the planar interface between these regions. In the guiding region, we expand the transformed electric field E(Q, 3, z) in a Taylor series with respect to z about its value at its boundary. To ensure the convergence of the Taylor expansion, the guiding region is partitioned in the z direction into subregions TzDVz2,*.*. separated by planar boundaries at z = zo,z1,..., within

24 which the electric field is continuous with respect to z and has the following expansion: 00 (-Z)~ On_ E(,,z) n! nE(, Zk)) z E zk+l, k =0, 1,... n=o (2.27) In the case of multilayered strips or channels, more than one Taylor expansion may be needed. Here, the refractive index profile of the guiding region is assumed to be a continuous function of z over the entire D, and the electric field is expanded in a single Taylor series at z = zo = 0. Note that even in the absence of any field discontinuity inside the guiding region, the convergence of the Taylor expansion must be verified. For all of the structures studied in the present work, this single expansion is valid for slab thicknesses of up to several free-space wavelengths, i.e. much larger than the practical values. Even for thicker dielectric regions, the convergence of Taylor expansions can be ensured by partitioning the region into several subregions. The higher-order derivatives appearing in (2.27) can now be related to the transformed electric field and its z derivatives just outside the guiding region in the enclosing layer via higher-order boundary conditions, which can be written in the following dyadic form: On region is then expressed in terms of the transform-domain unknownsh in the following form: oo M Zn E(,, z) =-jkoZo E E n. Ln j(,/3,0-). j((,), zD,. n=O j=l (2.29)

25 The volume polarization current and the electric field inside the guiding region are related in the space domain via definition (2.19). To relate the g-transforms of these quantities, use is made of the properties of the g-transform as outlined in the beginning of this section. Since J(x, y,z) = 0 for x ~ D: and z' DZ, one can write: Jp(, /, z) = jkoYo I S(7,', z)E(', /, z)d', (2.30) where the kernel S(5,', z) is given by S(,, z) = Kp D [r9g(x,z) - erb ]g(x)g*('x)dx. (2.31) Finally in view of (2.25) and (2.29)-(2.31), the following system of coupled integral equations are obtained: M Ji((=E Sj K( I ).J83(()d z=,..., M, j=l o0 (2.32) where the 3-dependence has been dropped for convenience, and Gij(' | /') are modified transform-domain dyadic Green's functions given by 00 Gij( I') = ko E RZ(, 5, ) Ln. ) j(E',, 0), (2.33) n=O with 7', (',/;3) being analytical integrals of known functions given by Z,', /3) = (I S', z )hi((z). dz. (2.34) Equation (2.32) represents a system of one-dimensional second-kind Fredholm integral equations for the transform-domain unknowns Jij(g). This system provides an exact one-dimensional formulation of the related two-dimensional boundary-value problem in the transform domain. Using a proper numerical technique such as the method of moments, one can reduce this system to a matrix eigenvalue problem,

26 from which the propagation constant of of the waveguide is found. The field distribution inside and outside the guiding region can then be computed by taking inverse g-transforms of (2.29) and (2.26), respectively. As mentioned earlier, in the most general case, we have a system of two coupled one-dimensional integral equations. The numerical treatment of such a system is by far simpler than the two-dimensional integral equations of conventional space-domain formulations. 2.4 The Modified Green's Functions It is important to note that due to the inclusion of all Taylor expansion terms, the modified Green's functions of (2.33) contain infinite summations. If analytical expressions can be found for the integrals 7Z, then the infinite summations can be rearranged and converted into known series expansions, and subsequently, closedform expressions for the modified Green's functions can be derived. In many two and three-dimensional practical waveguide structures, the guiding region is made up of one or more rectangular dielectric regions with a uniform index profile all over the cross section, i.e. erg(r) = Erg. For such geometries, the appropriate integral transform in the x direction is naturally a Fourier transform with g(6x) = exp(jfx). An inspection of equations (2.31) and (2.34) reveals that the infinite summations encountered in the modified Green's functions depend only on the guiding region and not on the geometry of the substrate structure. This implies that many planar dielectric structures share the same infinite summations, whose conversion into closed forms can be a tedious analytical task. Therefore, this section is devoted to the derivation of closed-form expressions for the modified dyadic Green's functions of a generic rectangular uniform guiding region.

27 For an open layered background geometry, the spectral decomposition (2.23) of the dyadic Green's function can be expressed in terms of the following constituent functions: ~ (3,, Cz) =- (,I)e (8z, (2.35) and h~((z)= e~bZ', (2.36) where (b and C, are the z-directed eigenvalues in the background layer enclosing the guiding region and a neighboring substrate layer, respectively, and satisfy equation (2.24). The spectral dyadic functions a (J) depend on the substrate configuration. The dyadic differential operator L, in this case reduces to the following form: Ln(=) - (n/2] Ln-2tn/2](), n > 2, (2.37) with Lo and L1 depending on the substrate geometry, and (g being the z-directed eigenvalue in the guiding region. Assuming a rectangular guiding region of dimensions a x b with a constant permittivity rg, the integrals 17 of (2.34) can then be written as the product of two separate functions: sin(e - J1)a/2 S({, (') = (rg - rb) sin(- a/2 (2.38) rb (5 -')' and (b ~Z it R =(nC) _ e(Z dz (h ( r l-e~bbE (Go l} (2.39) I /=0 L I

28 Now, the modified Green's functions can be written the following form: G,(s I = S(i,0') [ ( A( ) + 441(, ]) A1(('), (2.40) where 00 kI4(4, ) = Z,12m +k (), k 0,1, (2.41) m=O and Ak(') = ko2 Lk('. ('), k = 0,1, (2.42) with i,j = 1,2, and the values 1 and 2 of the indices i and j corresponding to the + and - signs, respectively. In view of equations (2.39) and (2.41), it can be shown that the functions Ok involve the following infinite sums: M-1 \2m S1 lim E M-oo0 m=0 =o \1 2 [cosh b) sinh(2b) 3 =t~ lim E1 2mE (C 2 (1 b)' M-T, oo M12 2 C b=o 2M" [2 [cosh(2b) + sinh(C2b) - eClb () ] (2.43)

29 Here, a limiting process has been invoked to ensure the convergence of the infinite sums for any possible range of the parameters. Setting 1I = ~(b, (2 = a, and using the above expansions with M - oo, the following expressions for Q/ are obtained: =() = _-Fe { - (cosh Cb T: sinh ) eCbb S6 S^ (b> / 1 J 6,) _ 1-icosh'b b Sinh b e ~bb} (2.44) Equation (2.40) along with (2.42) and (2.44) give explicit closed-form expressions for the modified transform-domain Green's functions of an open planar dielectric waveguide with a rectangular homogeneous guiding region and an arbitrary layered substrate geometry. Note that the complexity of the substrate structure manifests itself through the spectral functions I (i), appearing only in Ak, and the functions Y/ containing the infinite summations are independent of the substrate geometry. In the remainder of this section, we examine the dyadic functions Ak for a number of interesting practical waveguides. However, before proceeding to this task, first we study two special cases. By letting a -+ oo, the dielectric waveguide of Fig. 2.2 reduces to a one-dimensional layered geometry with infinite extent along the x and y directions. In this case, equation (2.38) approaches the following limit: lim S((, ) = (rg - Erb) 6( -'), (2.45) and the integral equations reduce to decoupled transcendental algebraic equations, which give the exact eigenvalues of the TE and TM modes of the one-dimensional infinite dielectric layer.

30 The next special case is when the guiding region is embedded in an infinite or semi-infinite cover medium. The latter case also includes a waveguide resting on the interface between a layered substrate structure and an infinite half-space. In such geometries, the equivalent polarization currents (sources) lie in an unbounded medium. Then, the Fourier transform of the dyadic Green's function of the background geometry, inside the substrate layer neighboring the interface, consists of a single spectral term with the constituent function a ((, f,,z) and h- (iz') defined by (2.35)-(2.36), where Cc and C, are the z-directed eigenvalues in the cover medium and the neighboring substrate layer, respectively. Therefore, M = 1, and the system of integral equations (2.32) reduces to a single integral equation given by J()= f G ( I'). J8(e')d', (2.46) J-00 where the transform-domain unknown J8,() is defined in the following manner: J(5)= i Jp(,, z) e-Cczdz. (2.47). —00 Thus, the h-transform in the z direction for this geometry is a bilateral Laplace transform. All of the geometries considered in the remainder of this section fall into this special class. 2.4.1 Dielectric Slab Waveguide The first example considered is a rectangular homogeneous dielectric slab of relative permittivity,,rg which is immersed in an unbounded cover medium of relative permittivity rc,. From an inspection of the free-space Green's function given in appendix A, it is seen that the dyadic function 4 for z < 0 is given by __= (- eCcz ~FP~cf = e ~J'__\~ r~f CD - (0 AQ\

31 where the spectral function D (F,/3, rc) is defined by equation (A.11). In this case. the modified Green's function is given by equation (2.40) with i = j = 2. The functions 0o and $1 are given by (2.44) with the lower signs, and the dyadic functions Ao and Al reduce to the following forms: ko2 _ ~12,/ ko -2/rc -^3/Crc -JfC/Erc jfC/c -J/gC'/.cg (6'2 + 2)/rg k2 _ f12/ er -glP Erg -jfg2/(ergc') - Al'/ -'//~rg -j _ I22/( 2 Cb) ) = | -'3/Erg ko - /r2/ -jg /(Cr') -J/('C/Cr -J/C/Erc ('12 + /2)/C) (2.49) 2.4.2 Dielectric Strip Waveguide The next example considered in this section is a dielectric strip waveguide consisting of a rectangular uniform slab of dimensions a x b and relative permittivity erg, resting over a single-layer conductor-backed substrate of thickness d and relative permittivity cr,. The geometry of the waveguide is shown in Fig. 2.3, where the ground plane is located at z = -d. In view of the Fourier transform of the Green's function of this background geometry, as given in appendix A, the dyadic function b for z < 0 can be written in the following form: (~,~,,) 1 = (,zd 4(5/, ('z) = D (, /, Ers,). (6 3) eCd(z+d) + 2 D (, /, cs).~ (, /) e( "+d, (2.50)

32 z Erc - a "' - I a x d rrs Figure 2.3: Geometry of a dielectric strip waveguide. where the spectral functions D (, 3, 6rs) are defined by equation (A.11), and the dyadic functions g(1,/ ) are given by 0 0 DTE sinh Cad ^ W )- 0 I =9 (6)0 =DTE sinh Csd _-J(crs -erc) -J(_rs - rc) Ers DTEDTM cosh Cd DTEDTM cosh (Cd DTM cosh(sd (2.51) with DTE and DTM defined by (A.18). In this case, the dyadic functions Ao and Al reduce to the following form: 1 =+ +. + C2d Ak( )= ko L.D. g e + 2 ko Lk.-. e, k = 0,, (2.52) where L c = k 0, 1 are the zeroth- and first-order dyadic boundary conditions given by 1 0 0 Lo (') O 1 o 0 0 1 0 9rg

33:(' 0 j'(1 - -) ~L1(~^) = O o jo ((l- ) (2.53) As expected, the expressions for Ao and A1 are quite complicated, and this complexity continues to grow as the number of substrate layers are increased. 2.4.3 Leaky Dielectric Waveguides The dielectric strip waveguide of Fig. 2.3 and other similar dielectric structures with an infinite layered substrate geometry are capable of supporting leaky modes. These modes deliver energy away from the axis of the waveguide. As a result, even in the case of completely lossless dielectric materials and perfectly conducting grounds, the propagation constant of the mode becomes a complex number, implying wave attenuation in the direction of propagation [42, 43]. It has been shown that if the propagation constant of the waveguide structure is less than that of a substrate mode, then energy leaks from the guiding region into the infinite substrate. Since background structures such as the single-layer conductor-backed substrate of Fig. 2.3 always support at least one surface wave mode with a zero cut-off, there is a good chance of encountering a leaky mode in this type of structures. In this section, we limit our attention to the dielectric strip waveguide of Fig. 2.3 and only consider leakage to the surface wave modes of the background geometry. The mathematical emergence of surface leaky modes can be attributed to the

34 behavior of the surface wave poles of the Green's functions of the structure. These poles are the roots of the following characteristic equations: DTM(,P) = 0, DTE(~,) = 0, (2.54) where DTM and DTE are defined in (A.18). These poles are also identical to the TM and TE eigenvalues of an infinite conductor-backed dielectric layer (see appendix A). This infinite geometry has a dominant mode of TMo with a zero cut-off. Therefore, the non-leaky modes of the dielectric strip waveguide of Fig. 2.3 are characterized by the following condition: / > /TMo> (2.55) where /3MO is the propagation constant of the substrate dominant mode [44]. In this case, it can be shown that the surface wave poles all lie on the imaginary axis, and the real-axis integration of the Sommerfeld integrals can be carried out without any difficulty. However, when 3 < O3Mo, the surface-wave poles begin to migrate toward the real axis, as shown in Fig. 2.4. The number of real-axis poles is determined by the number of the substrate surface wave modes whose propagation constants are larger than 3. To account for the wave attenuation along the waveguide axis, a small imaginary part is added to the waveguide propagation constant (/ - P - ja). This causes an expulsion of the leaky poles off the real axis and into the complex plane. The direction of expulsion is chosen so as to satisfy the radiation condition. Now, the integration contour must be properly deformed to enclose the leaky surface wave poles. Note that all those surface wave poles in the complex plane, for which /3s < l7e(/3), are ignored in the course of contour deformation. Fig. 2.4 shows a deformed integration contour for a typical leaky strip waveguide.

35 Im(e) branch point / non-leaky pole leaky pole - ~[' 0O Re(4) Figure 2.4: Contour deformation in the complex plane for leaky dielectric waveguides.

36 From equations (2.51) and (2.52), it is seen that the TM and TE surface wave poles of the substrate geometry appear in the spectral dyadic functions Ao and Ai. The real-axis integration of these functions is carried out by using the technique of extraction of singularities, which is discussed in appendix B. 2.4.4 Coupled Dielectric Waveguides The general methodology developed in section 2.3 is based on an arbitrary refractive index profile within the guiding region of the planar structure. By considering a piecewise-constant index profile along the x direction, we can easily extend the results of the preceding section to coupled coplanar dielectric waveguides. Here again we assume a uniform index profile along the z-axis; in other words, crg(r) = crg(x). In this case, S((, J', z) of (2.31) reduces to a convolutional kernel S(< - c'), independent of z, with S(<) given by the following Fourier transform: 1 S(s) = 2F { [ E(X) - erc ] P (X), (2.56) where p, (x) is the characteristic function of domain DZ, defined to have unit magnitude over this domain and vanish outside it. Now consider an array of N = 2N' coplanar, rectangular, homogeneous, dielectric slabs or strips of dimensions a x b, as shown in Fig. 2.5, which are uniformly spaced along the x-axis and separated from each other by a distance s. The guiding region of this structure can be regarded as a single rectangular region of dimensions [(N- l)s + a] x b embedded in the background geometry and with the following piecewise-constant index profile:,rg, Ix~ (2k-1)s/2 1<a/2, k- 1,...N', erg(X) = (2.57) erc, elsewhere.

37 z Erc a a T 2 Gometryonary xi Figure 2.5: Geometry of an array of coplanar dielectric waveguides. In view of (2.56)-(2.57) and the trigonometric identity: cos(2k -1) = 2sin 2N'58 k= ss ss s the function S,(S, (') associated with the waveguide array can be expressed as the following product: Sa(, = Se(,') Fa(,'6), (2.59) where Se(~, ~') is the S-function of a single element given by (2.38), and Fa(6, 6') is an array factor given by sin(6 -')Ns/2 Fo(,')- sin( -'s/2' (2.60) Thus, the two-dimensional problem of an array of dielectric waveguides reduces again to the one-dimensional integral equation (2.46), with exactly the same modified Green's function except for the S-function now given by (2.59). A special case of the array discussed above is a coupled dielectric waveguide with N = 2. This structure, which is of central importance in the design of submillimeter

38 wave and optical directional couplers, supports symmetric and antisymmetric propagating modes. In this case, the array factor (2.60) reduces to Fa(, ) - = 2cos(~ -')s/2. (2.61) 2.5 Numerical Solution Using the Method of Moments The system of the reduced one-dimensional integral equations derived in section 2.3 is now solved numerically using the method of moments to determine the transform-domain unknowns Jji, and subsequently the field distribution in the structure. It was seen in section 2.4 that for many practical structures, the system of integral equations reduces to a single equation given by (2.46). To simplify the notation, we will limit our attention to this very case in the remainder of this chapter. To implement the method of moments, the transform-domain unknown is expanded in an appropriate complete orthonormal basis {fn(()}n=1,2,... in the following manner: n=1 Then, by applying the Galerkin technique for testing, the following system of linear matrix equations are obtained: N Kmn.an = 0, m=1,...,N, (2.63) n=l where Kmn = 1 fm()G I |)fn(l)dd - mnI, (2.64) -00 -00 with Smn being the Kronecker delta. The propagation constant P of the waveguide is found by solving the following characteristic equation: det ([K(/3)]) =0. (2.65)

39 The proper choice of the expansion basis in a moment method solution depends mainly on the boundary conditions of the problem. In two-dimensional space-domain formulations, either the polarization current or the electric field inside the fguiding region is discretized by use of a suitable subsectional basis, which explicitly delineates the boundaries of the guiding region. In the integral transform technique, the geometry of the guiding region including its domains D3, and DTz and the relevant boundary conditions are directly incorporated into the modified transform-domain Green's functions. Therefore, for the expansion of the transform-domain unknowns, one requires a complete orthonormal basis which provides a good approximation of J.(0) in the (-domain. To examine the functional behavior of the transform-domain unknown, definition (2.47) is rewritten in the following manner: JV(0) = j L Jp(2,z)eAx/-; 2' dxdz, (2.66) 0-oo -oo where a2 = 32 - crc k2. It is not difficult to show that the inverse Fourier transform of (2.66) is given by J(W) = ()] a/2 K1 (c[ (x -')2 + z'2]1/2 ) = aZ'/2 J(x',z')dx'dz', [(, - 5')2 + z'21/2 (2.67) where K (x) is the first-order modified Bessel function of the second kind. Fig. 2.6 shows typical normalized components of J1(.) and its inverse Fourier transform Js(x) for a single slab and two coupled slabs of rectangular cross sections, assuming a uniform polarization current of unit magnitude and a prescribed value of f. It is seen that J1 extends over the entire real line but decays rapidly in both <- and x-domains.

40 Geometry J^(() Js(x) a single rectangular slab a -a/2 a/2 x two coupled rectangular slabs b e,.g \ s/2 s/2 X Figure 2.6: A typical component of the transform-domain unknown J.(<) and its inverse Fourier transform

41 In the latter domain, J.(x) indeed extends well beyond the boundaries of the guiding region. The rapid decay in this domain can be attributed to the convolution with the modified Bessel function. Thus, for the expansion of J.(0) in the (-domain, we need a complete orthonormal set of entire-domain basis functions, which span the space of square-integrable functions L2(R). In the meanwhile, the right candidate must be well localized in both transform and space domains with sufficient decay properties. Another consideration in the choice of the basis functions comes from the fact that the integral transform technique makes use of higher-order boundary conditions where the fields undergo abrupt discontinuities. In order to avoid non-physical solutions due to singular field derivatives, an entire-domain basis with sufficient smoothness is required that can replace the abrupt discontinuities by smooth transitions with any degree of accuracy. Considering all these requirements, the set of orthogonal Hermite-Gauss functions has been chosen as an expansion basis for this problem. These functions are defined in the following manner: 1n(X) = Mn e-2/2 HH(x) n = 0,1,..., (2.68) where H(x), n = 0, 1,... are the Hermite orthonormal polynomials [45], and Mn is a normalization factor. These functions are the exact eigensolutions of the wave equation in a medium with an untruncated parabolic refractive index profile [46]. The first four Hermite-Gauss functions are depicted in Fig. 2.7. They exhibit an oscillatory behavior in the vicinity of the origin and rapidly decay outside this region. Moreover, 7ln(-x) = (-1)n71n(z). The recurrence relation for these orthonormal functions is given by n.. = -x H, -- 1 H n-H 1, n > 1, (2.69) n+1 +1 -

42 Ho(x)'7H(x) 7H2(X) H33(x) 0\ ~7 Figure 2.7: Orthogonal Hermite-Gauss functions. with H-_(x) = 0, and Ho(x) = 7r-~/4. This particular normalization has been adopted to avoid the computer overflow error for large n [47]. An interesting property of Hermite-Gauss functions is that their Fourier transforms are also Hermite-Gaussian: Ann(x) < V/ jn,(~). (2.70) In order to obtain a reasonable approximation of J8(0) given a fixed number of the basis functions, one should scale the Hermite-Gauss functions appropriately. Introducing a scaling parameter do, expansion (2.62) can be'written in the following form: N J = E_ an VoWn(0do). (2.71) n=O The Galerkin technique guarantees the convergence of the solution for a sufficiently large number of expansion terms. Yet, to increase the numerical efficiency, the parameter do can be properly optimized [48]. The intuitive choice of do = a/2 usually

43 leads to fast convergence. Note that since in most practical waveguide configurations, the guiding region(s) is(are) symmetric with respect to x = 0, these structures support propagating modes whose field profiles are even or odd functions of x. Hence, expansion (2.71) can be confined to the even- or odd-order Hermite-Gauss functions for the respective components of J,. 2.6 Numerical Results and Discussion As for the validation of the integral transform technique, this section presents some numerical results for the dielectric waveguides treated in section 2.4. The simplicity of these structures allows for a fair comparison between this technique and other available numerical methods. In order to investigate the range of validity of the present formulation, rectangular dielectric waveguides of different sizes, aspect ratios and material parameters have been examined. The propagation constant of the waveguide is found from a numerical solution of the characteristic equation (2.65). The doublyinfinite integrals of (2.64) are computed using either a Hermite-Gauss quadrature or a high-order Legendre-Gauss quadrature after reduction to semi-infinite integrals based on symmetry considerations [49]. Note that in the latter quadrature, the integrals must be truncated appropriately. Due to the decay properties of Hermite-Gauss functions, these integrals converge very fast. The numerical experiments indicate a very rapid convergence in the solution with only a few basis functions. This convergence efficiency may be attributed to the superb ability of the Hermite-Gauss functions to closely approximate the transformdomain unknowns. In most cases, 2-5 even- or odd-order basis functions for each component of J., requiring a total of 6-15 unknowns, provide very accurate results

44 Method Number of unknowns Type of basis or discretizatiol Marcatili [1] 3 closed-form approximation Integral Transform [this work] 6-12 entire-domain basis Domain Integral Equation [5] 45-90 2-D pulse basis Finite Difference [4] 120-210 (for a quadrant) rectangular mesh Finite Element [7] 234 (for a quadrant) triangular elements Table 2.1: A comparison between the integral transform technique and other methods in view of numerical efficiency within less than 1% error. Table 2.1 makes a comparison of different techniques in view of computational efficiency. This table shows the typical numbers of unknowns required by these techniques for the solution of the eigenvalue problem associated with a rectangular dielectric slab. It is seen that the methods based on the discretization of the cross-section of the waveguide involve huge numbers of unknowns. Thus, by reducing the dimensionality of the boundary-value problem in the transform domain and avoiding the discretization of the geometry, the integral transform technique achieves a high degree of numerical efficiency. Note that in addition to the reduction in the number of unknowns, the rapid convergence of the involving integrals also contribute to a significant reduction in the computation time. Fig. 2.8 shows the variation of the normalized propagation constants of the dominant x- and z-polarized modes versus a normalized frequency defined as B = (2b/Xo)(erg - ~rc)1/2 for a rect

45 angular dielectric slab with erg = 2.25, er = 1, and R = a/b = 2. These results are in good agreement with those of [5] based on a two-dimensional domain integral equation. This figure also shows the results of Marcatili's approximate method [1], which is seen to fail as the guide dimensions begin to shrink. Fig. 2.9 shows the convergence of the solutions versus the number of basis functions for B = 1. In Fig. 2.10, the normalized propagation constants of some of the higher-order modes are shown for the same dielectric slab as in Fig. 2.8. In the submillimeter-wave integrated circuit technology, higher dielectric constants, especially of semiconductor materials, are of more practical interest. Fig. 2.11 shows the normalized propagation constant of the dominant x-polarized mode for a typical high-contrast semiconductor slab with Erg = 13.1, c 1, and R = a/b =3. The results are compared to those of Marcatili's approximation [1] and the frequencydomain finite difference method [4], and excellent agreement is observed. The next example is a coupled dielectric slab waveguide as shown in Fig. 2.5 with N = 2. This structure supports symmetric and antisymmetric modes, whose propagation constants are denoted by Pf and Pa, respectively. Fig. 2.12 shows the variation of the normalized propagation constants of the dominant symmetric and antisymmetric x-polarized modes versus the proximity ratio s/a, for two identical dielectric slabs with Erg = 2.25, erc = 1, R = a/b = 2, and B = 1. This figure also compares our results with those of Marcatili's approximation [1]. It is seen that as the spacing between the two slabs increases, the symmetric and antisymmetric modes become degenerate and approach the dominant mode of a single isolated slab. To study the performance of the integral transform technique for substrate structures, first a non-leaky dielectric strip waveguide of the type shown in Fig. 2.3 is considered with Erg = 3.8, Ers = 1.5, Erc = 1, a/b = 2.25, and d/b = 0.5. Fig. 2.13

46 1.50........ 1.40 - aa.. Ex 1.30 1.20 This work 1.10 //.. Marcatili ASi i So Domain IE 1.00 0.00 0.50 1.00 1.50 2.00 B Figure 2.8: Normalized propagation constants of dominant E1i and EJi modes of a rectangular dielectric slab as a function of normalized frequency B.

47 1.50 1.40 >1 1.30 Ex _11 1.20 Ey11 1.10 1.00,, 0 2 4 6 8 10 N Figure 2.9: Convergence of the normalized propagation constants of dominant modes of a rectangular dielectric slab with respect to the number of basis functions.

48 1.50.......... 1.40 - a E E2 1.30 *... A... /X../....E. 1.20 - E12 1.1o 0 E:2 1.00 0.00 0.50 1.00 1.50 2.00 B Figure 2.10: Normalized propagation constants of higher-order modes of a rectangular dielectric slab as a function of normalized frequency B.

49 4.00 3.50 a 3.00 2.50 2.00 /:r/ ~ —-- This method 1.50 -. Marcatili f E Finite difference 1.00 0.00 0.50 1.00 1.50 2.00 B Figure 2.11: Normalized propagation constant of the dominant modes of a highcontrast dielectric slab as a function of normalized frequency B.

50 1.50.. I 1.40 fo symmetric 1.30 -..,..........."" ---— " 5 1.20'" v Cl,,' antisymmetric 1.10 - Marcatili A This work 1.00 1.00 1.20 1.40 1.60 1.80 2.00 s/a Figure 2.12: Normalized propagation constants of dominant symmetric and antisymmetric modes of a coupled dielectric slab waveguide as a function proximity ratio s/a.

51 shows the variation of the normalized propagation constants of the dominant x- and z-polarized modes versus normalized strip thickness kob for this planar waveguide. The results are compared to those of a frequency-domain finite difference method [8], and excellent agreement is observed. It is seen that the Ez1 mode with a vertical polarization has a zero cut-off, while the dispersion curve of the EfZ mode with a horizontal polarization starts from a cut-off and after a cross-over region takes dominance over the other mode. The last example of this section is concerned with the study of leakage effects. Consider the geometry of the dielectric strip waveguide of Fig. 2.3 with Erg = 2.62,,rs = 2.55, cr, = 1, a = 10mm, d = 3.8mm, and a varying strip thickness b. The first two dominant modes of this structure are Ejl and Ef1. When b = 0, the geometry reduces to the infinite substrate layer, and the two above-mentioned modes degenerate to the substrate TE1 and TMo modes, respectively. The El1 mode of the strip waveguide, which is TM-like, is always non-leaky, and its propagation constant is larger than that of the dominant substrate mode, i.e., the TMo mode, for any strip thickness. However, the propagation constant of the E1l mode is smaller than that of the dominant substrate mode, and this mode always leaks to the substrate TMo mode for any strip thickness. As a result, the propagation constant of the Ef1 mode is always complex with an imaginary part that represents wave attenuation along the waveguide axis. Figs. 2.14 and 2.15 show the real and imaginary parts of the normalized propagation constants of the dominant x- and z-polarized modes versus strip thickness b [mm] at a frequency of f = 30 GHz (A0 = 10 mm). The results show a good agreement with those based on a domain integral equation technique [20].

52 2.00.... a 1.80 d9 Ex11 1.60 1.40 1.20 o —- Finite Difference 0 This work 1.00.00.50 1.00 1.50 2.00 2.50 3.00 k0b Figure 2.13: Normalized propagation constants of dominant ET and Ez1 modes of a non-leaky dielectric strip waveguide as a function of normalized strip thickness.

53 1.60...... i|... |.... i... 1.... 1.60 EZ11 1.50 TMo(substrate) 1.40 1.30 TE, (substrate) - EFIE o This work 1.20.00.10.20.30.40.50.60 b [mm] Figure 2.14: Normalized propagation constants of the dominant El and E1l modes of a leaky dielectric strip waveguide as a function of strip thickness [in mm] at f = 30 GHz.

54.000E+00... -1.00 -2.00 - -3.00 - -4.00 - EFIE EFIE / Em This work -5.00..00.10.20.30.40.50.60 b [mm] Figure 2.15: Attenuation constant of the leaky E-T mode of the dielectric strip waveguide shown in Fig. 2.14 as a function of strip thickness [in mm] at f = 30 GHz.

55 2.7 Concluding Remarks In this chapter, the integral transform technique was presented as an accurate, efficient, full-wave method for the study of two-dimensional planar dielectric structures. This technique is applicable to a wide variety of layered waveguide geometries which are of practical importance in submillimeter-wave and integrated optics applications. The essence of this novel approach lies in the reduction of the dimensionality of the boundary-value problem by invoking appropriate integral transforms. It can easily be extended to three-dimensional structures with corresponding two-dimensional reduced transform-domain formulations. The numerical experiments clearly reveal the efficiency of this technique from a computational point of view. It might appear that this high degree of numerical efficiency is acquired at the expense of heavy analytical preprocessing involving higher-order boundary conditions, infinite summations, and tedious derivations of closed-form expressions for the modified Green's functions. However, note that the dyadic Green's functios of most practical layered geometries have similar spectral decompositions, generally made up of real or complex exponential functions. Therefore, the higher-order boundary conditions reduce more or less to the same simple forms discussed before, and except for the cases of anomalous dielectric inhomogeneities, the infinite summations may vary very slightly from one geometry to another. Finally, it should be noted that the present technique, with all its merits, is limited to geometries with planar boundaries. The applicability of this technique strongly depends on the derivation of closed-form expressions for the modified Green's functions. This task in not generally feasible for arbitrary geometries or index profiles. In the case of curved boundaries or more complex and irregular shapes, one still has

56 to resort to the methods based on a fine disretization of the geometry. It is therefore concluded that the integral transform technique is a highly geometry-dependent technique which offers a very high degree of computational efficiency within its limits of applicability.

CHAPTER III MULTIRESOLUTION EXPANSIONS FOR MOMENT METHOD SOLUTION OF ELECTROMAGNETIC PROBLEMS 3.1 Introduction The Method of Moments (MoM) has long been used as one of the primary tools for the rigorous study of electromagnetic problems [24]. This method in conjunction with the integral equation technique provides the most accurate and in many cases an exact formulation of the related the boundary value problem. In spite of all its advantages over the other numerically intensive techniques, the major challenge in the method of moments is the choice of appropriate basis functions that would provide the required accuracy and computational efficiency. An impeding characteristic of this approach is that the resulting matrices are full and have very poor condition numbers. This challenge becomes even bigger as the geometry under study becomes non-canonical and complex, involving large numbers of unknowns. In the previous chapter, a novel integral formulation was developed which leads to the derivation of transform-domain integral equations of reduced dimensionality for the relevant boundary value problem. This approach results in a reduction of the number of unknowns involved in the problem. The use of the method of mo57

58 ments for the numerical solution of the reduced integral equations does not pose any particular problem due to the relatively small sizes of the matrices. However, the limitation of the applicability of this technique in view of geometry dependence was also emphasized. The full-wave analysis of complex electromagnetic structures requires a fine discretization of the geometry under study, whereupon a rigorous numerical scheme is implemented to solve Maxwell's equations. As the complexity of the structure slightly increases, the numerical treatment of the problem can easily turn into a huge, expensive and time-consuming computational task. The differential-based techniques like finite element and finite difference methods are capable of handling large-scale problems thanks to their inherent sparsity properties, even though at the expense of extremely intensive computational procedures. The integral-based techniques, however, are deprived of this sparsity privilege. All the conventional basis functions normally used in the method of moments generate full moment matrices. The computational problems associated with the storage and manipulation of large, densely populated matrices easily rules out the practicality of the integral equation technique in competition with the other numerically intensive methods mentioned earlier. This bottleneck has traditionally limited the application of the method of moments to small-scale or medium-scale electromagnetic problems. In the final years of last decade, the newly developed theory of wavelet transforms introduced an entirely new class of expansions, called discrete wavelet expansions, with some extraordinary properties which were long sought after in various fields of computational sciences [50, 51]. A very salient feature of these new expansions was that the entire set of basis functions is generated by the dilation and shifting of a single function, called the mother wavelet. More prominent was the nice localization of the new expansion functions in both time and frequency domains. The potential

59 application of wavelet theory in the numerical solution of integral equations was soon exploited and led to the finding that the wavelet expansion of certain types of integral operators generates highly sparse linear systems [25]. The past two years have witnessed the successful application of the wavelet theory to various electromagnetic problems [32, 33, 23, 34]. It has been demonstrated that the use of wavelet expansions in conjunction with the method of moments leads to highly sparse moment matrices. This property is based on the fact that the resulting matrices enjoy very good condition numbers so that performing a thresholding procedure by discarding a large number of the matrix entries which fall below a certain threshold level does not degrade the accuracy of results. In essence, expansion in wavelets and the subsequent thresholding process amounts to a regularization of the system, by which ill-conditioned components are discarded at the price of loosing some details. Such a sparsity makes it possible to employ a wide range of numerical methods especially developed for the fast and efficient storage and inversion of this type of matrices. Orthonormal wavelets have some other interesting properties which can also be properly exploited to speed up the moment method computations significantly. One such feature is the Fast Wavelet Algorithm (FWA) which has been used extensively in signal processing applications [27]. In brief, this property enables one to compute the wavelet expansion coefficients at a certain resolution from those of the next higher resolution. Thus, it is possible to construct a recursive algorithm to compute all the moment matrix elements in consecutive resolution levels from those of the highest resolution. This chapter starts with a short basic account of the theory of one-dimensional Multiresolution Analysis (MRA), where the construction of a one-dimensional mul

60 tiresolution expansion for the approximation of an arbitrary function is contemplated. Then, the fast wavelet algorithm is described in its one-dimensional version. This is followed by the construction of a two-dimensional multiresolution analysis and the associated FWA algorithm. Finally, the adoption of multiresolution expansions for the moment method analysis of electromagnetic problems and the implementation of the fast wavelet algorithm for this purpose are discussed. 3.2 One-Dimensional Multiresolution Analysis A multiresolution analysis (MRA) of L2(R) is a sequence of nested approximation subspaces of L2(R), denoted by {Vm}meZ, such that all subspaces are scaled versions of a central space Vo obtained through dilation by integral powers of 2 [50]. Each subspace Vm,, which is said to have a resolution of 2-m, possesses details twice as fine as those of its predecessor, Vm-1, on this approximation ladder. In mathematical terms, the following properties must hold: *V C Vm+l, VmE Z, * UmEz Vm is dense in L2(R), * nEZ Vm = 0, * f(x) E Vm, f(2x) E Vm+l, Vm E Z, * f(x) E Vm X f(x- 2-mn) E Vm, m,n E Z. where Z is the set of integers. It has been shown that there exists a characteristic function +(x) E L2(R), called the scaling function, such that bm,n(z) = 2m/2q$(2mx - n), n E Z is an orthonormal

(il1 ~(x)i v(x) 1 _ 1 ____Ws~~~____ __ 1/2 0 1 x 0 lx Figure 3.1: Scaling function and mother wavelet of the Haar multiresolution analysis. basis of Vm. The functions qm,n are in fact properly dilated and shifted versions of the scaling function. We denote the orthogonal complement of Vm in Vm+i by Wm, and write Vm @ Wm = Vm+. (3.1) It can be shown that there exists another characteristic function 4(x) e L2(R), called the mother wavelet, such that om,n(x) = 2m/20 (2mx - n), n e Z is an orthonormal basis of the complement space Wm. From the properties listed above, it follows that ( Wm = L2(R), (3.2) mEZ implying that the set of functions {fm,n}m,lEZ is a complete orthonormal basis of L2(R). The simplest multiresolution analysis is generated by Haar functions. Fig. 3.1 shows the scaling function and mother wavelet of the Haar multiresolution analysis. Given an arbitrary square-integrable function f(x) C L2(R), one can express its approximation at a resolution of 2-m by defining a projection operator Pm(f) onto the subspace Vm in the following form: Pm(f) = E < f, m,n > )m,n((), (3.3) nEZ

62 where < f,g > denotes the inner product of the two functions f(x.r an-d g(.r). Fronm equation (3.1), one can then write: Pm+l(f) = Pm(f) + Z < f, bm,n > bm,n(X). (3.4) nEZ It is clear that at each resolution level, the approximation of the given function in the form of (3.3) results in losing some fine details, which of course can be restored in the next higher-level approximation subspace. As evidenced by equation (3.4), these lost details may be thought of as resting in the complement space l7m. In view of the above properties, it can be shown that having a crude approximation of the function f(x) at a coarse resolution level, say mi, one can improve it to a finer resolution level, say m2, by exploiting the wavelets of the intermediate levels in the following manner: m2 -1 Pm2(f)= P1 (f) + E ]E < f, %m,n > Vbm,n(zX) m2 > ml. (3.5) m=m1 nEZ Note that for the expansion functions involved in equation (3.5), the following orthogonality relation holds: < >ml,l,?Pm2,n >=0, im2 > mi, 1,n E Z. (3.6) Various multiresolution analyses with different mathematical properties have been introduced in recent years [50]. In this work, we adopt the classical Battle-Lemarie MRA, which is constructed based on polynomial spline functions. Note the Haar multiresolution analysis is the first member of this family of MRAs. Figs. 3.2 and 3.3 show the plots of the cubic spline Battle-Lemarie scaling function and mother wavelet, respectively. The construction of these functions is discussed in detail in appendix C. The Battle-Lemarie MRA has an exponential decay, which makes it very efficient in view of numerical convergence issues. For practical purposes, these functions are assumed to have an effective support (width) over [-6,6] and [-5.5,6.5], respectively

6:3 2.0 0................................... 1.50 4(X) 1.00.50.00 -.50 -1.00.... -6.00-5.00-4.00-3.00-2.00-1.00.00 1.00 2.00 3.00 4.00 5.00 6.00 x Figure 3.2: The cubic spline Battle-Lemarie scaling function. 2.00................. 1.50 v (X) 1.00.50 -1.00 ---—. -6.00-5.00-4.00-3.00-2.00-1.00.00 1.00 2.00 3.00 4.00 5.00 6.00 x Figure 3.3: The cubic spline Battle-Lemarie mother wavelet.

61 [32]. Although this kind of localization cannot match the actual coi~lipactiness of soIic other types of wavelets such as Daubechies wavelets [26], it is quite sufficient for our computational needs. In exchange, the Battle-Lemarie functions enjoy a perfect symmetry that compactly supported wavelets are inherently deprived of. Such a symmetry is a big plus in the numerical evaluation of moment integrals. In particular, when dealing with geometries which possess some sort of symmetry, using BattleLemarie functions amounts to significant computational savings. By combining the symmetry of basis functions with the inherent symmetry of the Green's functions, computation of many moment integrals which are either equal or differ only in a sign can be spared. Yet, the biggest advantage of Battle-Lemarie functions is the availability of simple closed-form expressions for their Fourier transforms, even if they are truncated. This point will become more clear later in chapter 5. As mentioned earlier, multiresolution functions are highly localized in both time and frequency domains, or more appropriately in the context of the present work, in both space and spectral domains. Figs. 3.4 and 3.5 show the Fourier transforms of the cubic spline Battle-Lemarie scaling function and mother wavelet, respectively. Other types of MRAs have a similar Fourier-domain behavior. As seen from these figures, the mother wavelet is a bandpass function with no dc component, while the scaling function is a lowpass function with its spectrum centered around the origin. In fact, the definition of a multiresolution analysis requires that roo / (x)dx = 1, -00 r00 (x)dx = O. (3.7) -00 The localization of wavelets in both domains can be visualized through the lattice shown in Fig. 3.6. In this figure, the horizontal and vertical axes represent the space

6.5 1.50... l, 1.00.50.00.. -15.00 -10.00 -5.00.00 5.00 10.00 15.00 Figure 3.4: Fourier transform of the cubic spline Battle-Lemarie scaling function. 1.50 1.00.50.00 -15.00 -10.00 -5.00.00 5.00 10.00 15.00 Figure 3.5: Fourier transform of the cubic spline Battle-Lemarie mother wavelet.

66: 2-m 22 m=l ~:e ~ ~ ~ ~ e —e ~i m=0 * * * * _ * * * 0: CO m=O m=-1 * * - _ _ _-I ~' x -2 -1 0- 2 1 1 2.* *.. * *- --- * * * -2wV region of interest Figure 3.6: The wavelet localization lattice. and spectral domains, respectively. Each point represents the center of a wavelet along the (horizontal) x-axis and the center of its spectrum along the (vertical) <axis. Horizontal levels correspond to different resolutions levels (m), and the points on each horizontal level represent different value of the shift index n. As the resolution increases, so do the concentration of shifted wavelets inside a given interval. The center of their spectrum also shifts away from the origin. 3.3 The Fast Wavelet Algorithm In the beginning of this chapter it was mentioned that the fast wavelet algorithm (FWA) is an effective tool for speeding up the computation of the moment integrals. This algorithm is founded upon the following two-scale property of the multiresolution

(i7 analysis: [50] (x) = v/2 hno (2x-n ), nEZ +(x) = X/2 En g (2x- n), (3.8) nEZ where {hn}nEZ is a zero-centered symmetric discrete sequence with fast decay, and n = (-l)n-1hl_,. (3.9) The values of the first few positive hn coefficients for a cubic spline Battle-Lemarie multiresolution analysis are given in appendix A. In view of equations (3.8), it is not difficult to derive the two following relations: < f, m,n > = hk-2n < f,/m+l,k >, kEZ < f, m,n > = g,-2n < f, >m+l,k > (3.10) kEZ where * denotes the complex conjugate. To cast into a signal flow diagram, we introduce the following expansion coefficients: roo cm,n = < f,/m,n > = J f(zx)>m,n(x)dx, J -00 dm,n = < fm,n > = f(x)bmn,n(x)dx. ~~-~ ~(3oo11) Then, equations (3.10) take the following simple forms: Cm,n = h _k-2n Cm+l,k, keZ dm,n = 5 -2n Cm+l,k, (3.12) keZ which is the mathematical statement of the fast wavelet algorithm. Fig. 3.7 shows the signal flow diagram for the fast wavelet algorithm, where the filters H and G

6(S Cm+1 H cm H m-1 H m2 G G G dm m-1 m-2 Figure 3.7: Signal flow graph for the 1-D FWA. involve a discrete convolution plus a decimation by 2. Note that according to equations (3.12), the scaling and wavelet coefficients at each resolution can be computed by a discrete convolution from the knowledge of only scaling coefficients at the next higher resolution. When several resolution levels are required, it is therefore sufficient to compute the scaling coefficients only at the resolution level immediately following the highest resolution involved using the first equation of (3.11). All the expansion coefficients at the lower resolution levels can be evaluated recursively using equations (3.12). 3.4 Two-Dimensional Multiresolution Analysis The concept of a one-dimensional multiresolution analysis can readily be extended to the space of two-dimensional square integrable functions. By analogy, a twodimensional multiresolution analysis is a ladder of successive nested approximation subspaces Vm of L2(R2). Each subspace Vm at a resolution of 2-m is constructed from the Cartesian product of two one-dimensional MRAs at the same resolution, i.e., Vm = Vm 0 Vm [50]. An orthonormal basis of Vm would be comprised of the

i9 functions mn,nt,ny(X, y) = Cm,n(Z)qPm.7.(y) 7nlr., E Z, where the one-dimensional functions in x and y are properly dilated and shifted versions of the scaling functions O(x) or 4(y), respectively, as described in section 3.2. The orthogonal complement of Vm in Vm+1 is denoted by Wm, and it can easily be shown that Wm = (Wm Vm) (Vm 0 W/) (Wmi 0 4l), (3.13) where Wm is the one-dimensional complement space of Vm. One can deduce from equation (3.13) that the orthonormal basis of Wm consists of three sets of wavelets: m,ny,,( ) = Imn,()m)n^ nny), (( y) = <m,n (z)m,n()y and In(, (y) = 0m,nx(X)0'm,ny(y), with +'(x) and'(y) being the mother wavelet. These wavelets are called vertical, horizontal and diagonal, respectively, in that they sample variations of the expanded function primarily along the corresponding directions. It should be mentioned here that this is not the only way of constructing a twodimensional multiresolution analysis. Another option, which is often called the standard construction, involves Cartesian products of only the complement spaces Wm but at different resolutions. The resulting 2-D wavelets include all the possible products bmx,nx,(X)'m.,ny(Y). It can be shown that the set {,nmx,,nx,my,ny}mi,ni,my,nyEz is a complete orthonormal set of L2(R2). Also, due to the two-scale property (3.8), it is possible to show that the two standard and non-standard constructions are related to each other. However, for reasons which will be stated in the next section, we adopt the first (non-standard) construction scheme. The approximation of a square-integrable function f(x,y) at a resolution of 2-m can be expressed in terms of a projection operator onto the subspace Vm defined in the following way: Pm(f) = E E < fm,n,,ny,n~ > (i)m,nx,ny(X,y). (3.14) nxEZ nyEZ

7( Then, the improved approximation at the next filler resolution. 2-"'-I involves the wavelets of the complement subspace Wm and is givenl by Pm+(f) = Pm(f) + E E < f, lm,n,ny > qm,n,?n,,(xy). i=v,h,d nxEZ nyEZ (3.15) The improvement of the approximation can be continued to any arbitrary resolution by including the intermediate 2-D wavelets in the following way: m2 -1 Pm (f) Pm, (f)+ E E E E < fm,nny > m,n,n (X y) m=mi i=v,h,d nxEZ nyEZ m2 > mi. (3.16) Equation (3.16) is the two-dimensional counterpart of equation (3.5). The two-dimensional fast wavelet algorithm is a direct extension of the onedimensional version. In the 2-D case we have four types of expansion coefficients, which are defined in the following form: roo r00oo Cmnxny = < fm,nx,ny > = J J f(z, y)4m,nx,ny(X, y)dxdy, -00 -00 00 00 Dm,nny = < fm,n,,ny > = /00 / f(Z Y)fm,ny,ny(X, y)dxdy, i =v, h,d, J-OO J -OO (3.17) where the superscripts v, h, and d stand for vertical, horizontal and diagonal, respectively. Then, the mathematical statement of the 2-D FWA is given by Cm,ni=,ny E 1 hk -2nhky-2nyhC*l~kT ky: kxEZ kyEZ m,nx,ny = 9kx-2nxky-2ny Cm+l,k,,ky, kEZ kyEZ Dmn,,ny = 5 C hkl-2nx9ky-2ny Cm+l,kx,ky? kxEZ kyEZ

71 Cm+1 H HH Cmi HH Cm- HH2 GH GH GH DVm DVm DVm-2 HG HG HG Dhm Dhm1 Dhm2 GG GG GG rm Dm-1 m-2 Ddm Ddmi Ddm2 Figure 3.8: Signal flow graph for the 2-D FWA. Dm,n,ny E g kx-2nxgky-2ny Cm+l,kx,ky kxEZ kyEZ (3.18) Fig. 3.8 shows the signal flow diagram for the two-dimensional fast wavelet algorithm. Here again, when different resolution levels are required, it suffices to compute only the 2-D scaling coefficients at the resolution level immediately following the highest resolution involved using the first equation of (3.17). All the expansion coefficients at the lower resolution levels can be evaluated recursively using equations (3.18). 3.5 Construction of Moment Method Expansions From the previous sections, we conclude that the multiresolution analysis theory presents a systematic approach to the approximation of square-integrable functions. This approximation is in the form of expansion in a set of orthonormal basis functions.

72 Note that the resulting expansion basis is a doubly-indexed set whliclh includes dilated and shifted wavelets and scaling functions. In the method of moments, the unknown fields or currents are expanded in anl appropriate set of linearly independent basis functions. By replacing the fields or currents with their discrete expansions, the integral equations are thus discretized. Then, using a suitable testing procedure such as the Galerkin technique, in which the test functions are chosen to be identical to the expansion functions, leads to a system of linear matrix equations. By solving this linear system, the unknown expansion coefficients and subsequently the field distribution in the structure are determined. The accuracy and efficiency of the method of moments largely depends on the proper choice of the basis functions. These functions must be able to provide a fair approximation of the unknown quantities and also satisfy the boundary conditions of the problem. In this section, we describe how to construct one- or two-dimensional multiresolution expansions for the moment method treatment of electromagnetic fields or currents. The completeness of the orthonormal set {1Im,n}m,nEZ intuitively suggests the use of a pure wavelet expansion in the method of moments. Such an expansion is theoretically capable of approximating a given square-integrable function to any arbitrary degree of accuracy by choosing very large values of m and n. However, as pointed out earlier, wavelets are bandpass functions with no low-frequency spectral content. By contrast, electromagnetic fields and currents often have considerable lowfrequency spectra, which correspond to the uniform or slowly varying distributions. Note that the term "spectrum" here is equivalent to the spatial frequency content. To restore the lowpass component of the fields or currents, we make use of equation (3.5) to construct the moment method expansion. Thus, the proper expansion has

the following form: f(X) = E Cmo,, m k() + E E dn, (3.19) kEZ do m=mo nEZ do where mo is an initial resolution level, and do is the characteristic length of the problem. This expansion contains both lowpass scaling functions and bandpass wavelets, and is therefore capable of restoring the complete spectral content of the fields or currents. Equation (3.19) may also be interpreted as follows: First the given unknown quantity is approximated with an initial resolution of 2-mo using the scaling functions at this resolution. This initial approximation, which contains all the details with resolutions equal to or greater than 2-mo, is then improved to any arbitrary degree of accuracy by bringing in the wavelets at the initial and subsequent higher resolutions. The choice of the parameter do depends on the problem. For instance, in a space-domain formulation, do is taken to be an effective wavelength. It is very important to note that many electromagnetic entities like currents are often confined to finite regions and vanish outside those regions. By contrast, regular multiresolution expansions, and in particular, the Battle-Lemarie functions, which are used in the present work, are defined over the entire real line. In order to have a good approximation of the fields or currents in the vicinity of the boundaries and discontinuities, it is inevitable to place some basis functions outside these boundaries. Then, to satisfy the boundary conditions of the problem, the expansion must be explicitly enforced to vanish in the exterior regions. This amounts to an additional set of equations which are solved simultaneously with the original integral equations. After the reduction of the integral equations into a linear system by applying the Galerkin technique, the total number of equations will of course be greater than the number of unknowns. To have a consistent system, it is therefore necessary to devise

74 a slightly modified version of the Galerkin technique, in which a smaller llIlnumber ol equations exactly equal to the number of unknowns are generated. This goal can l)e achieved by partitioning the set of basis functions into smaller subsets and then use them appropriately to test different sets of equations. The details of this technique will be illustrated in the following chapters. Note that the need for placement of some multiresolution basis functions outside the domain of the problem is not due to the infinite support of Battle-Lemarie functions, which are to be used in the present work. Even compactly supported basis functions like Daubechies' wavelets still suffer from the same problem. In general, this problem arises whenever we try to represent a function which is defined over a finite interval by basis functions which are defined over the entire real line and have overlapping supports. Another way to avoid this difficulty is to construct multiresolution expansions which are naturally defined over a finite interval [50]. The construction of this type of expansions, however, is extremely complicated, and their use is not expected to offer any computational gain. Yet, another option is the use of periodized wavelets, which are obtained from ordinary wavelets through a periodization process [50]. The basic drawback of this latter alternative is the loss of the shift properties of multiresolution analysis. In other words, the periodized wavelets are no longer obtained by simply dilating and shifting a mother wavelet, and they are more or less elaborate entiredomain basis functions. In addition, the fast wavelet algorithm is no longer applicable to the periodized wavelets. The two-dimensional multiresolution expansions for electromagnetic fields and currents are constructed in a similar fashion like the one-dimensional case. To restore the complete spectral content of the unknown quantities, the following expansion is

adopted: f/(X,) = C o,l;,k,ko mo,k,kyy( T ) k:EZkyEZ do do 00 + E E E D Dmnx,ny m n,ny" (d- y)d m=mo i=v,h,d nEZ nyEZ do d (3.20) where mo is again the initial resolution level. From equations (3.19) and (3.20) it is apparent that a multiresolution basis can indeed contain a very large number of expansion functions. The double-index pair (m, n), identifying the resolution and shift amount, may lead to a very complicated index notation for the moment matrix elements. To facilitate this task, we adopt a simpler notation in the following manner: f (X, y) = aiFFi( (3.21) do' do with a similar counterpart for the one-dimensional case. In this notation, each single index i is identified with a dilation index mi, two shift indices n,- and n^, and the types of generating functions for x and y variables, which are either the scaling function ) or the mother wavelet b. Note that when the expanded quantity is defined over a finite domain and vanishes outside this domain, the number of shift indices will also be finite and can easily be determined from the wavelet localization lattice of Fig. 3.6. The moment integrals involve twice as many variables as the dimension of the integral equation, counting for the source and observation variables. In other words, for a one-dimensional boundary value problem, the integration is performed with respect to two variables x and x'. Similarly, in a two-dimensional problem, the moment integrals involve four variables x,y,x' and y'. For instance, a typical

7 6 two-dimensional moment integral has the following form: Kij =j j Fi( )Ge(x,Y I x', y')F( -)dxdydx'dy' (3.22) do do do do where Ge(x, y I x', y') is the dyadic Green's function of the problem, and Fi(x, y) and Fj(z', y') are the test and expansion functions, respectively. Equation (3.22) uses the shorthand notation of (3.21); otherwise, it would have been written in the following manner:,, fiyifarj,,,fy Kij K7 m,,,ny jm,. (3.23) where the superscripts f,, f xy, f, fy, denote the generating functions for the respective variables. The number of basis functions can still be reduced in an efficient way without sacrificing the accuracy by a careful exploitation of the inherent properties of wavelets. Equation (3.7) states that a wavelet has a zero average. This means that wavelets are basically suitable for the sampling of functions where they undergo abrupt discontinuities or rapid variations. In other words, the wavelet coefficients of a uniform or slowly varying function are insignificant. This property enables us to economize the moment method expansion by placing wavelet basis functions only at places where the fields or currents are expected to exhibit some sort of singularity or discontinuity, or in general where more details of the expanded function are needed. Thus, in the one-dimensional expansion of (3.19), the 1-D wavelets will be more concentrated close to the boundaries of the problem. In the two-dimensional case, i.e. equation (3.20), this stipulation takes an interesting form as follows: Recalling that there are three types of 2-D wavelets and considering their generating functions, it is concluded that the horizontal, vertical and diagonal wavelets must be primarily concentrated

77 Td Th Figure 3.9: Typical locations of 2-D multiresoluiton basis functions. around the horizontal, vertical and diagonal boundaries, respectively. For rectangular boundaries, the last type of the wavelets play a major part in the modeling of field or current distributions at the corners of the region. Fig. 3.9 shows typical locations for the placement of different types of wavelets in a rectangular region. We close this section with a comment on the implementation of the fast wavelet algorithm for the computation of moment integrals. As mentioned in earlier sections, this algorithm enables us to spare a large part of numerical integrations. It is sufficient to evaluate the moment integrals for the scaling functions at the resolution level immediately following the highest resolution involved. Let m, denote the highest resolution required for the problem under study. Then, we will compute only the Again, note that the ranges of the shift indices are finite and are determined with the aid of the wavelet localization lattice of Fig. 3.6. It is clear that these ranges vary with the size of the geometry under study. Thus, only those shift indices whose

78 i-ts effective support (width) falls within the domain of the problem or in the vicinity of its boundaries are retained for computations. Then, all other moment integrals can be computed recursively using the following relation: Kij= E E E k, EZ ky z EZ k% EZ k eZ x,* yi* Xj* yl* 7k. -2nz, 7ky -2nyi -2n ky. -2nyi mi+l,kci,kyi lmjt+l,,j,ky (3.25) where k ii,7ki,7Y,7k, are either hk or gk depending on the type of the respective generating functions. Note that the quadruple sum in equation (3.25) involves products of four hk or gk coefficients, which become negligibly small for the majority of index combinations. One can therefore disregard all index combinations which would produce very small coefficient products. 3.6 Wavelets as Electromagnetic Radiators The preceding section deliberated on the ways to economize multiresolution expansions for the method of moments. It was pointed out that some of the wavelet basis functions which do not contain significant information can be omitted from the expansion bases. Yet, the size of the bases and subsequently the number of unknowns can become very large, especially for 3-D electromagnetic problems where a 2-D multiresolution expansion will be employed. In the beginning of this chapter, we emphasized that the main motivation for the use of wavelet expansions is the generation of sparse moment matrices. This raises the question whether this sparsity can be predicted in advance, so that the computation of the negligible moment integrals can be spared. It should be noted that the sparsity resulting from the use of

i79 wavelet expansions in the method of moments is different in natre from the sparsit encountered in differential-based techniques. In the latter techniques, we have actual zero elements, while the sparsity of the moment matrices is a result, of performing a thresholding procedure, whereby all matrix entries which fall below a certain threshold level are discarded. Therefore, prediction of the sparsity pattern in advance will of course depend on the threshold level and seems to be an impossible task in general. However, it is possible to obtain a qualitative picture of the moment method interactions between wavelet basis functions. To this end, we investigate the behavior of a wavelet as an electromagnetic radiator. The moment integral of equation (3.22) can be interpreted as the reaction of the test function Fi(x, y) to the electric field radiated by the expansion function Fj(x, y). This equation can be rewritten in the following manner: Kij J, (do ) Ej(xo y) dxdy (3.26) where E(x,y) = j Ge(X,y I Xyl) F d-)dz'dy'. (3.27) It is therefore useful to study the electric field radiated by a multiresolution function, and in particular, investigate its far field. In view of equation (A.5) of appendix A and from an examination of the dyadic Green's functions given in this appendix, one can write: Ge(, y, z x', y,O) = (2 Ge ( ) j(x+y) -jr(sincos+ sinsid (.28z= ei(^x'+?1y') e-jr( sin ecos k+'i sin sin sk-j(c cos e) dCdr1, (3.28)

S() where (r, 0, <) are the usual spherical coordinates describing the observation point. and C2 = <2 + r72 - k. In the far field, the integral of (3.2S) takes an asymptotic form which can be derived using the method of stationary phase as discussed in appendix D. It can be shown that the stationary point of this integral is given by,s = ko sin 0 cos X, r7 = ko sin 0 sin. (3.29) where ko is the free-space propagation constant. For large values of kor, equation (3.28) can be approximated by the following asymptotic expression [52, 53]: Ge jkocos2 e-jkor Ge(s s) j ko(sinsinycos.+s (3.30) 27rr z/'=0 Then, the far field radiated by a multiresolution function reduces to the following form: Ejff = jkocosO k __ jk0COS~ 03kor Ejff ---- Ge(sris) 27rr Z'=o J /Si Fi(^ d d)eiko (sin 0 cosx'+sin si) (3.31) fj~. F,(, (-'Oy" dx'dy'. (3.31) But the double integral in equation (3.31) is recognized to be the two-dimensional Fourier transform of the multiresolution function Fj(x, y) evaluated at the stationary point. This function can be written explicitly in the following form: F (x, y) = f()(x). f(Y)(y) (3.32) where fig) and fY)(y) are either a dilated and shifted scaling function or a dilated and shifted wavelet depending on the type of the 2-D wavelet. Then, considering the dilation factors of multiresolution functions and the normalization parameter do, it

Sl is a straightforward task to show that the double integral of (3.31) reduces to the following expression:;/ fs ("') oc d2 f-j(2-mj kodo sin 0 cos ). fj(2-" k^odo sin 0 sin ), (3.33) Jj (...) X j (3.33) where fxj(x) and fyj(y) are either the scaling function (p or mother wavelet V/, and fxj and fyj are 4 or 4 are their Fourier transforms, respectively. The argument of these functions takes a maximum value of (2w (3.34) 2m, where we have assumed that the parameter do is an effective waelength with do < XA. Now, recall that the mother wavelet is a bandpass function with no spectral content around the origin, i.e. in(t) = 0 for small s i. Since for 2-D wavelets, at least one of fxj or fyj is always a mother wavelet, there exists a resolution level mff, depending on the type of the multiresolution analysis, such that for m > mrff the far fields radiated by wavelets are always zero anywhere in space. For the cubic spline Battle-Lemarie wavelets, this resolution level is mff = 2 when do = Ao is chosen. Returning to equation (3.26), it is seen that if the test function Fi(x, y) happens to fall in the far-field region of the expansion function Fj(x, y), then the corresponding moment matrix element will be negligible and is expected to be discarded after thresholding. Hence, this interaction may be skipped in the process of numerical integrations. The choice of the far-field range is of course a subjective one and depends on the required accuracy as well as the threshold level. In the present work, this range is determined experimentally according to the traditional criterion: 2D2 r > 2, (3.35) where D is the largest dimension of the wavelet usually taken to be 2-m I,, with hI being the effective support (width) of the wavelet.

CHAPTER IV SPECTRAL-DOMAIN WAVELET ANALYSIS OF INTEGRATED PLANAR WAVEGUIDES 4.1 Introduction The significance of rigorous full-wave analysis of integrated planar waveguides was articulated in chapter 1. It was pointed out that by increasing the complexity of the geometry under study, the numerical modeling of this type of structures can easily turn into a difficult computational task. In chapter 2, we developed a special technique which leads to a major reduction of number of unknowns through the formulation of integral equations of reduced dimensionality. We also emphasized that, though computationally very efficient, this technique practically works only for a certain class of planar dielectric structures with limited canonical shapes. The accurate modeling of general non-canonical problems inevitably requires a very fine discretization of the geometry of the structure, which in turn may amount to prohibitively large numbers of unknowns. The advantages of integral-based techniques over differential-based formulations in view of mesh sizes and ability to accurately model open-type structures were outlined in earlier chapters. Yet, we also recounted in chapter 3 the special problems encountered in the moment method treatment of large problems with big matrix sizes. It was mentioned that in the past couple of years, the theory of orthonormal 82

wavelets has been successfully exploited in the study of integral operators to generatehighly sparse moment matrices. The present work is one of the pioneering efforts to explore the application of wavelet theory in computational electromagnetics. Chronologically, we started out with incorporating the concepts of multidimensional analysis into a spectral-domain formulation of planar dielectric waveguides. The reason for this starting point was the natural continuation of the integral transform technique developed earlier during this research effort. It was pointed out the due to the special nature of that technique, the moment method solution of the reduced integral equations calls for a special entiredomain basis with good localization in both space and spatial frequency domains as well as sufficient decay and regularity. Such a basis should extend over the entire line. The options for this purpose are not many, and this search led the author to the newly developed orthonormal wavelet expansions, which were introduced in the final years of the last decade. Nonetheless, the orthogonal Hermite-Gauss functions, which share many similar features with wavelets, were eventually chosen as the expansion basis for the integral transform technique [48]. The search for a general, efficient, and accurate technique for the modeling of planar microwave structures, and in particular, open geometries, suggested the idea of incorporating wavelet concepts into the method of moments. In this way, it was hoped to produce sparse linear systems which can be handled numerically with much better ease, while maintaining the inherent advantages of integral-based approaches. The wavelet formulation was first applied to simple purely dielectric waveguides like those studied in chapter 2, and it proved very successful. Then, metallized structures and subsequently hybrid structures comprising dielectric and metallic sections were investigated, and a general two-dimensional spectral-domain formulation was

developed for the analysis of integrated planar waveguides. In this chapter, first we present the general formulation of the 2-D integral equations for an arbitrary number of dielectric and metallic sections with an arbitrary planar substrate configuration. In section 3, some special considerations for the use of multiresolution expansions in a spectral-domain formulation are discussed. The moment method implementation and related computational details are the subject of section 4. Finally, extensive numerical results are presented for a variety of integrated planar dielectric waveguides, which are of practical interest in microwave, millimeterwave and optical applications. 4.2 Formulation of Integral Equations In this section, we derive a system of spectral-domain integral equations to describe a general two-dimensional integrated planar waveguide. Before proceeding to this task, let us first introduce some definitions and conventions which will be used throughout this chapter. By Dz and Dz we denote domains in the x and z directions, respectively. The characteristic function of a domain Dt is defined as i, t E Dt, pv(t)= { (4.1) O O, otherwise. We introduce a Fourier transform along the x-axis in the following manner: f() = [f(x)] = / f(x)ejxdx. (4.2) With every x-directed domain DX, we associate a spectral-domain function S-D,(<) defined as: SI (~=) =- [pmD(x)]. (4.3)

85 A function f(x) is said to have a compact support over thie doimain Dr if f(x) = O for x f ZDx. The Fourier transform of such a function satisfies the following integral equation: f(0) = J Sv(( - (')(4) d'. (4.4) - 00 Given an arbitrary spectral-domain function 4)(<), we define its band-limited version over the domain D: by the following convolution integral: COo WvD [()(] = SvD ( -')(I(') d'. (4.5) -00 Then, the inverse Fourier transform of the band-limited version of 4)(() has a compact support over D, [54]. Fig. 4.1 shows the geometry of a general integrated planar waveguide, which consists of an arbitrary number of rectangular dielectric slabs and metallic strips, embedded in an arbitrary planar, piecewise homogeneous, multilayered, superstrate/substrate background structure. We assume that the waveguide structure supports a propagating mode of propagation constant P along the y-axis, and a time-harmonic dependence of the form ejwt is assumed and suppressed. The background layers are located parallel to the x-axis, with an infinite extent in this direction. The nth background layer is characterized by a constant relative permittivity crsn and a z-directed domain DSe, implying that z(s) < z < z() where z = z" ) and z() are the boundaries of this n - n+1 I n n+1 layer with its neighboring layers. The ith dielectric slab is characterized by its index contrast with respect to its surrounding background layer, which is denoted by 6cr~(z) and is assumed to vary along the z direction. The domains of this slab along the xand -axes are denoted by )(P) and D(P), respectively. The kth metallic strip, located at z = zc, is identified by the domain D(C). All permittivities are assumed to be, in

metallic strips a dielectric strips w /i rse i/' Figure 4.1: Geometry of an integrated planar waveguide structure.

87 general, complex to include material losses, and all metallic strips are assumed to be made of perfect conductors with zero thicknesses. It is clear that a lossy metallic strip with a nonzero thickness can be modeled as a lossy slab with a high loss tangent. The total electric field E(r) can be considered as radiated by two types of contributing currents: volume polarization currents, Jpi(r), which are defined over the cross sections of dielectric slabs as: Jpi(r) =jkoYo 6eri(z)E(r)p(p)(x)p (P)(z), -o0 < < y 00, (4.6) and surface conduction currents, Jck(), which are sustained over the surface of metallic strips. In definition (4.6), Yo = 1/Zo and ko denote the free-space characteristic admittance and propagation constant, respectively. Then, the electric field at any point can be expressed in the following form: E(r) = -jkoZo E J Ge(i i r'). Jp(x', y', z')dx'dy'dz' 3 x j3 -jkoZo E 1 GeJ(r I| r')z, = (C) Jc_(x', y)ddy', (4.7) -jkoZo>3 ~ jr(: O where the summations extend over all dielectric slabs and metallic strips, respectively, and Ge(r I r') denotes the electric field dyadic Green's function of the background structure, with its form depending on the background layers where the observation and source points are chosen. The infinite integration with respect to y' can be performed analytically. Then, according to appendix A, for the two-dimensional planar layered geometry of Fig 4.1, one can write G~(r 1 rI00 =) I Ge(r IT) 2= 1 Ge ( IP z zej (C-')df, L (x - x')(z- z') ^rsp. ^0

88 z (s) z' c S -o < y < (.S) where Ge (, 3, z | z') denotes the Fourier transform of the two-dimensional dyadic Green's function with the observation and source points located inside the,uth and vth background layers, respectively, 6,, is the Kronecker delta, and L = zz is the source dyadic which accounts for the singularity of the Green's function at the source point. In view of equations (4.7) and (4.8), one can express the transformed electric field in the,uth background layer in the following form: E(,P,z) = -jkoZo J E G ( /3 z I z'). Jpj ( z' )dz' 3 ZJ +3 1 E -L Jpj(t,/3, z) ~ i( ) ersp -jkoZo S Ge ((fltz I z(C)).JcL(,3), z E D,(s) (4.9) where the second summation involving the source dyadic extends over the dielectric slabs embedded inside the yuth layer, and the Green's functions are evaluated with the source points located in the appropriate enclosing layers. In view of definition (4.3), the Fourier transform of equation (4.6) can be written as.(fz) =jk e(z~)pp(z)() Sp(-)E(',3,z)d'. (4.10) zt -00 xi Then, from equations (4.9) and (4.10), the following set of integral equations for the transformed polarization currents are obtained: ( S+-i(z) ).A Jpi(s,z) = SE,;^iz)k2 J d'S,(SP -n) -

+ S GeQzlz~ z |c))' i ) ),, i..... (4.11) where I is the unit dyadic, and,,s. is the relative permittivity of the background layer surrounding the ith dielectric slab. To derive another set of integral equations for the transformed conduction currents, we make use of the boundary condition on the metallic strips that the tangential components of the electric field vanish on the surface of perfect conductors. This condition is implemented in the Fourier domain by way of an appropriate testing procedure. Let 4)(<) be an arbitrary spectral-domain function. Then, by testing the tangential components of equation (4.9) on the surface of the metallic strips at z = () with the band-limited version of <I>() over (C) one obtains: r00 -jkoZo 5 -'ok(p. (,, I')e Jp(, ( z')Wc[()]d[dz(C) 3 ZI oo =:k z?O) zI ~). -jkoZo e | zk I(J c[3(,) )I tX- y., x:00 xk (4.12) Note that in equation (4.12), the tangential unit vectors force the normal source dyadic terms to vanish. The integral on the left hand side of this equation involves the product of two spectral-domain functions whose inverse Fourier transforms have

90 non-overlapping supports in the space domain. In view of Parseval's theorem [54]. this integral vanishes, and equation (4.12) reduces to the following form: 0 = dJ d< S;(c)( - E')()i. 00 Or:rk o -o0 J-00 xk (E c)Ge (', Z, \Z'() (, Ji,' )dZ' ~3= 1 Jv1=1 + Ge e (^ ^ l 01 ^v) ~ / ) t = x,, k = 1,...,N, (4.13) where S*(g) denotes the complex conjugate of S(g) and use has been made of definition (4.5). A last set of boundary conditions which need a careful treatment in the spectral domain concern the compactness of polarization or conduction currents, requiring that Jpi(r) = 0 forx z D(), and JCk(x,y) = 0 for x DvkC. Comparing (4.5) and (4.11), it is seen that the latter equation indeed relates the transformed polarization currents to band-limited functions defined over the supports of these currents, thus immediately enforcing the compactness of this type of currents. Equation (4.13), however, does not guarantee the compactness of conduction currents by itself, and this condition must be implemented explicitly in the following form: Jck({1P)= J S(c)(- Jck( /)d' k= 1-.. i(C). (4.14) J-00 xk The set of equations (4.11), (4.13) and (4.14) now render a complete integral formulation of the integrated planar waveguide structure.

91 4.3 Multiresolution Expansion of Transformed Currents The construction of multiresolution expansions for the moment method treatment of electromagneitc problems was discussed in detail in chapter 3. It is important to note here that in the formulation of the preceding section, the unknown functions in the integral equations are the Fourier transforms of polarization and conduction currents. It is therefore these transform-domain quantities which are to be expanded in a proper basis to discretize the integral equations.This situation is rather different from that of the next chapter, in which the physical currents themselves are expanded in a multiresolution basis. In order to facilitate the notation and better illustrate the numerical procedure, we simplify the general geometry of Fig. 4.1 as follows. Consider the finite dielectric slab of dimensions a x b with the printed metallic strip of width w shown at the center of this figure. The dielectric slab is made up of a stack of integrated dielectric strips with different relative permittivities rgiI i = 1, 2,... and can be modeled by a piecewise constant index profile along the z-axis. Appendix A gives the expressions for the dyadic Green's function of an infinite two-layer substrate geometry with a perfectly conducting ground plane. For this simplified structure, the unknown functions include a transformed polarization current Jp(t,P,z) and a transformed conduction current J,((,P). The two-dimensional current Jp(, P, z) is expanded in a multiresolution expansion in the variable g and in a sub-domain pulse basis along the z-axis. For the latter expansion, the domain TVP) of the dielectric slab is divided into subintervals A, q = 1,..., Nb, with the associated characteristic functions pq(z), such that index profile of the slab at each sub-layer is a constant. The one-dimensional current Jc{( 3) is expanded in a multiresolution expansion in the variable (. To preserve orthogonality, an identical

92 initial resolution level is adopted for both multiresolution expansions. To siIplify index numbering, we use the one-dimensional counterpart of the shorthand notation introduced in equation (3.21) of chapter 3. In addition, we collect all the multiresolution basis functions associated with the transformed polarization current into a single set denoted by {f(P)Q()}, j = 1,..., Np,, and all the multiresolution basis functions associated with the transformed conduction current into another set denoted by {f/()()}, = 1, =,..., Nc. Then, the expansions for the two transformed currents are given by Jp(f, I, z) = (P) pq(Z), (4.15) j=1 q=l o and Jc(SP=) = al f (-), (4.16) l=1 k0 where the free-space propagation constant ko = 27r/Ao has been chosen for the normalization parameter. Note that the expansion coefficients are indeed functions of the propagation constant 3. Before implementing the method of moment, let us comment on some special features of using multiresolution expansions in a spectral-domain formulation. The multiresolution basis functions of equations (4.15) and (4.16) can be explicitly written in the following form: f i( k) ) -2'/2 ( 2 \ fi(o)= 2 i() (2m ( -ni). (4.17) Taking an inverse Fourier transform of equation (4.17) yields: Taking an inverse Fourier transform of equation (4.17) yields: ko~~X [^ )S

93 From the plots of the Fourier transforms of the scaling function and mother wavelet shown in Figs. 3.3 and 3.4, it is seen that the spectrum of +(x) has a bandwidth of 2go centered around the origin, while +(x) has a narrow passband centered around ~Aw,. Hence, in view of equation (4.18), the dilated scaling functions at the initial resolution level mo all have a bandwidth equal to 2mowsAo/7r. On the other hand, the polarization and conduction currents have compact supports in the space domain, which are equal to a and w, respectively, for the geometry under study. To get a good initial approximation of Jp and J we choose hoan initial resolution level mo such that the Fourier transform of the functions mo0,n(x) occupies a bandwidth comparable to a and w. Assuming that w < a, a rule of thumb is to have 2mown < 27r(w/Ao) < 2mo+lw,. The shift indices at each resolution level can be determined from the localization lattice of Fig. 3.5. Thus, only those wavelets whose spectra fall within the region of interest may be preserved for the moment method expansion. Finally, note that the present spectral-domain formulation calls for an interchange between the roles of the variables x and ( in this wavelet lattice. The fact that the spectral-domain unknowns are strictly band-limited due to the compactness of the space-domain currents places an upper bound on the highest resolution level required. 4.4 The Moment Method Solution and Numerical Considerations To implement the method of moments, expansions (4.15) and (4.16) are inserted into the integral equations (4.11), (4.13) and (4.14), and the resultant equations are tested by an appropriate set of functions to obtain a linear system of matrix equations. By solving this linear system, the unknown current amplitudes and subsequently the

94 field distributions are determined. The integral equation (4.11) corresponding to the polarization current is tested using the regular Galerkin procedure, in which the test functions are chosen to be {f(P)pp()}, i = 1,..., Np, p = 1..., Nb, i.e. identical to the expansion functions. Associated with the conduction current are the two independent integral equations (4.13) and (4.14), which require a split type of testing to produce an overall consistent linear system with the number of equations equal to the number of unknowns. This provision was already discussed in section 3.5 of the preceding chapter. To this end, we partition the expansion set {f(C)(g)}, k = 1,..., N, into two subsets with Na and Nc - Nc basis functions in each subset, respectively. Equation (4.13) is now tested using a modified version of Galerkin's technique, in which the test functions are taken to be the band-limited version of the first subset over the domain of the metallic strip, i.e. {W(C)[fc(c)(O)]}, k = 1,..., N.N This is equivalent to setting q(Q) = fC[)(E), k = 1,..., Nc, in equation (4.12). The first subset naturally includes those basis functions whose spectra are mainly confined inside the boundary of the problem such as dilated and shifted scaling functions at the initial resolution level. Finally, for equation (4.14), the conventional Galerkin testing is invoked again by choosing for the test functions the second subset {fkC)(f)}, k = NC+1,..., Nc. Thus, the following linear matrix system consisting of 3NpNb + 2Nc scalar algebraic equations is obtained: [R [R [a(P)] [R ] [R ] __ =0 (4.19) 0 [S] [a(C)] where [a(P) ] and [a(C)] are the unknown amplitude vectors of the transformed po

95 larization and conduction currents, respectively, [R(P] represents the interaction between the polarization current elements, [R ] represents the interaction between the conduction current elements, [Rc ] and [R( ] represent the cross-interaction between the polarization and conduction current elements, and [ S] accounts for the compactness of the conduction current. To illustrate the numerical evaluation of the moment integrals, let us consider the.= (pp) interactions [R ) ] for the case a multilayered dielectric strip. A typical integral of this type has the following form: Rjijq P - ipj + - L 1 A ij pq, (4.20) ipljq =,Pipljq ( Crc (4 where &6rp and Ap are the index contrast and thickness of the pth sub-layer (z-cell), respectively, and 0roo oo rb rb Pipljq = 6erpk J dd < dz dz' Sa(- ) fip)( )pp(z)Ge(',3, Z I zf')fj(()pq(Z') (4.21) with the band-limited kernel Sa defined in the following manner: Sa() = sn a2 (4.22) The quadruple integral of equation (4.21) can be rewritten as follows: Pipljq= I A(P)(, a) pq(,P)fjP) ()dl. (4.23) J -00 In this equation, AP) is a function of geometry defined as P)(,a)= J Sa,('- /(f' ()dJ'. (4.24) - 00

96 In view of equation (4.4). it is seen that the functions A(P) are in fact band-limited versions of the multiresolution basis functions f,(P) over the boundary of the dielectric slab. Note that the integral of equation (4.24) does not involve the Green's function and is independent of the propagation constant 3. According to appendix C, these functions can be evaluated semi-analytically for a Battle-Lemarie multiresolution basis by employing the cubic spline representations (C.9). Even the numerical integration of (4.24) is very efficient due to the fast decay of Battle-Lemarie functions. The dyadic functions 9pq in equation (4.23) are obtained by z-integration of the dyadic Green's function in the following manner: gpq(,c) = ro dz J dz'pp(z)Ge(,fz I z')pq(z'). (4.25) The above spatial integrations can be performed analytically to yield closed-form expressions for Gpq. According to appendix A, the Fourier transform of the dyadic electric field Green's function of the substrate geometry in the cover medium can be written in the following way: Ge(3,z |z') = D (3Cr0) 2 + 13) e(c+z) (4.26)( where the two terms represent the primary and secondary components, respectively, with the first term being identical to the Green's function of the unbounded space and gs depending on the substrate geometry. The dyadic functions D are defined by equation (A.11). The functions Gpq can be decomposed accordingly in the following manner: gpq(J, i) = s(, 3) + qQ(, fi) (4.27) Then, the two components are given by the following expressions: 5P,) = /2 [ (i, [l D)) )+ D (,,.) 1 ) —- - 2S L cp j

97 D (, i, ~+c) e-c(bpbq) ( e(cP -1 ) (1 - e-(cA ), ) > q, -P 6at k2 ~p(~,/3): 2C = 2 D (v, IC) e-,c(q-bp) (ercz- 1) (1P ), p < q, 5= 6 /) ec(bp+bq) ( e('cP - 1" ( c1~q- ) (4.28) where bpl and bp denote the z-coordinates of the lower and upper boundaries of the pth sub-layer of the dielectric slab (bo = 0). Note that the first equation in (4.28) corresponds to the self-cell elements. The remaining spectral-domain integrals, i.e. (4.23), are evaluated numerically using an efficient scheme such as the Gauss-Legendre quadrature. Thanks to the exponential decay of the Battle-Lemarie functions, these integrals are highly convergent, and the domain of integration is practically of finite size. Moreover, using the symmetry property of these functions, the numerical integration can be further simplified by converting the doubly infinite integrals into ones over [0, oo]. 4.5 Numerical Results By solving the eigensystem (4.19), one can determine the unknown amplitude vectors [a(P)] and [a() ] and subsequently the expansions (4.15) and (4.16) of the transformed polarization and conduction currents. Then, the electric field distribution can readily be computed at any point of the geometry by taking the inverse Fourier transform of equation (4.9). The eigenvalues of this system give the propagation constants of the propagating modes of the waveguide structure, which are found by

98 setting the determinant of the moment matrix equal to zero. It turns out that the resultant eigenmatrix is highly sparse in the sense that a very large majority of the matrix entries are extremely small compared to its largest entry. This phenomenon prompts one to delete a large portion of the moment matrix by performing a threshold procedure, whereby all those entries of the matrix whose ratios to the largest entry fall in magnitude below a certain threshold are set equal to zero. The reduced moment matrix is then very sparsely populated and can easily be handled in terms of storage and mathematical manipulations by special sparse techniques available for this kind of linear systems [55]. Pure dielectric structures with minimal metallization in the form of a ground plane are of primary interest in submillimeter-wave applications. For this type of waveguides, the formulation of section 4.2 simply reduces to integral equation (4.11) with JcI(.) = 0, and the eigensystem (4.19) consists only of the submatrix [R( ]. In this chapter, we confine ourselves to non-leaky propagating modes with real propagation constants, and search for real values of /, with 3 > pdom, where pdom is the propagation constant of the dominant mode of the infinite substrate geometry. The first example of this section consists of a single-strip dielectric waveguide resting upon a single-layer conductor-backed substrate. Fig. 4.2 shows the variation of normalized propagation constants of the first two Ex and Ez1 modes of this structure as a function of the normalized strip thickness. The waveguide parameters are erg3.8, Ers = 1.5, a = 2.25b, d = 0.5b, and w = 0. This figure also compares our results with those of the frequency-domain finite difference method [8], and excellent agreement is observed. To obtain these data, a multiresolution expansion with a total of 25-35 basis functions has been used. In all cases, only two successive resolution levels are required which depend upon the width of the dielectric strip. For example,

99 2.00.,.......... 1.75 d m 1.50 -' 1.25 E ---- This work 0 Finite difference 1.00.00.10.20.30.40.50 b/Xo Figure 4.2: Normalized propagation constants of dominant E' and Ez1 modes of a non-leaky single-strip dielectric waveguide with a single-layer substrate as a function of normalized strip thickness.

Threshold Sparsity Propagation constant Error percentage 0 - 1.32729 0 0.2% 86% 1.33372 0.48% 0.5% 91% 1.33680 0.71% 1% 95% 1.33956 0.92% Table 4.1: Effect of thresholding on the moment matrix. for a = 2.25b = 0.675Ao, the initial resolution level is mo = -1, with 9 scaling functions and 18 wavelets used at this resolution level plus 4 wavelets at m = 0. In the z direction, 1-2 sub-domain pulse functions usually give very accurate results for horizontally polarized modes, and 3-4 pulses are required for vertically polarized modes. Fig. 4.3 shows the structure of a typical moment matrix for this geometry, where a total of 31 multiresolution functions and 3 pulses have been used and a threshold of 1% has been applied. The sparsity of this moment matrix is 95%, meaning that 95% of the entire matrix entries have magnitudes less than 1% of the largest entry. Table 4.1 shows the effect of this thresholding on the computed value of the waveguide propagation constant for different threshold levels. It is seen that for a threshold of 1%, the error in the computation of 3 is less that 1% with respect to the case when no threshold is applied. The next example involves a single-strip dielectric waveguide with a two-layer grounded substrate. In this structure, large index contrasts typical of semicon

101 w,~w,,l ~ ~ ~~ir. qqqqqqe~/e~ — I ~~. ~ ~ I %..% *%%%%%% ~ %%%% *._...'q % %%%.%%..'". ~' 4%* ~~~~% %~~~~%%%...*%%% %%%. ~~~~%%% %%%*.00%%%%%% %%%. %%%~~~%..%%%%%. 4 %%%i % %%~~~. %~~%~~%%.%%%%% *%%%% * %,.. %%%%%%.%%~%%%%**%%%.,~~.~4~4%% %~l~r *%%%%%%.; %%%%% *.%%%.~~~%%%~~.* % % %%% *%:d. %%%% %%%..%%%% %%%% %%%% % *%4%....%%% ~ ~ ~~%% %% ~~irr ~ ~ ~ ~ ~ ~ ~ ~ r ~~~~ qqqq~bbb %%%* %% too %*. kq,~~HI%% %% e q~,qqqq.%, I~ %%%% %%%% %% * ~ *0 #.%%%*'% %%%/ %%~~~~~% *%%%%% ~% ZR ~~~~ %4.. %%%%%%~~~~~~% lb~~~qkqb ~ q~qbq qb~' ebbq I ~ ~ ~ ~ ~ ~ ~ qbbb~ kbbq SSllee, bbbq ~ ~ eqbbqqb~q ~ * * I~~~dA *.~~% %%%%%~~~o %*% *.'~:. *~........... *.__ ~~~~%% %%%%%* %%%%%%*"o %%% %%~~~% % ~ ~ ~ ~%** %%O %%%% %%%** %%% % %%*O%% *I% ****%%%%% ~ eq, qqq~q~ e~q qqbq %%% **%%%% ~ ~ %% %% % So.)~~~~~~% %%%%%%% %* %;4%%** %%%%% %** % %* %%%%%% %% Figure 4.3: The structure of moment matrix of Fig. 4.2 after applying a threshold of I %. rl_ ~ ~ ~ ~ ~ ~ ~ ~~'If' &~qb~' ~..I I~~~~~~ I'- it 0_'q~f" - ~q,,~qq-tq~-/.q I'A'""~L!~ ~ ~."A~IIA~. e ~ ~ *~*~~. Si~~~~~~~~~~~~~~~~~~ qbql"~fr lIq, I' =I _ ~i~1'. ~? ~ ~ ~ ~ ~~~~ ~ qbqb q bqrqbqkq b ~ ~ qbbb ql~'.q qbb bq, ~ ~ qsqbqb~b q, qbqbqsqb qbl'e eqbqbqbq ~z-~. _~]~ ~'q''~ ~~~~~~~~~~~~~~~~~ ~;'~ qb' qbqb~~~~~~qbqbqbqb.~ ~ ~ q, qkqbqbq, 2PJ"~q~qbqbqbq Figure 4.3: The struture of moment matri of Fig. 4.2 after bqlqbqb q lqbqlqb ~ ~ 1%.~~~~~~~~~~~~~~~~~~~~~~~l~bqqqq

102 ductor materials have been assumed. The waveguide parameters are 9rg =- rs2 = 12.85, Crsl = 10.0, a = 30.6pm, b = 58.24pm, di = 22.7pm, d2 = 17.1,im, and u = 0. Fig. 4.4 shows the variation of the normalized propagation constant of the dominant mode of this structure as a function of frequency. The results have been compared to those of the two-dimensional finite-difference time-domain (FDTD) method [56]. Due to the relatively large values of the strip thickness compared to the guide wavelength in this structure, more pulses (typically Nb > 5) are needed for the z-expansion of the polarization current. The formulation developed in this paper can also easily be applied to multilayered strip waveguides by assuming a piecewise-constant index profile along the z direction for the dielectric guiding region. Fig. 4.5 shows the frequency dependence of the normalized propagation constant of the dominant mode of a double-strip dielectric waveguide with a single substrate layer. The waveguide parameters in this figure are Crg1 = 10.0, Crg2 = Crs = 12.85, a = 30.6pm, bl = 22.7p/m, b = 80.9ipm, d = 17.1pm, and w = 0. This figure also compares our results to those of the mode matching technique [13]. At lower frequencies typical of microwave and millimeter-wave applications, microstripbased metallized devices play a very important role. It has been demonstrated that using an insulating dielectric strip between a microstrip and the supporting substrate structure can reduce the ohmic losses due to the conductors [57]. Metallized dielectric structures also have practical importance in waveguide transitions and feed networks. To validate the formulation developed in this paper for metallized structures, first a simple microstrip line printed on a two-layer conductor-backed substrate is considered. For this geometry, Jpj(, y) = 0, and the two integral equations (4.13) and (4.14) must be solved simultaneously, with the eigensystem (4.19) consisting only of the submatrices [R ] and [S ]. Fig.4.6 shows the variation of the normalized prop

103 2C.50aa 2.50 2.00 1.50'his work 0 FDTD 1.00...... I 250. 500. 750. 1000. 1250. f [GHz] Figure 4.4: Normalized propagation constant of the dominant mode of a single-strip dielectric waveguide with a double-layer substrate as a function of frequency.

104 3.50 a ~ 3.00 - bl t 2.50 2.00 1.50/ —- /This work 1.50 E Mode matching 1.00. I... I. I. 250. 500. 750. 1000. 1250. 1500. f [GHz] Figure 4.5: Normalized propagation constant of the dominant mode of a double-strip dielectric waveguide with a single-layer substrate as a function of frequency.

105 agation constant of the dominant mode versus frequency for a microstrip line with Crsl = 12.9, Crs2 = 8.2, w = 1.Omm, di = 0.58mm, d2 = 0.33mm, and b = 0. These results agree quite well with those based on the spectral-domain technique described in [11] using entire-domain Bessel functions for the expansion of the transformed microstrip current. Fig. 4.7 shows the structure of a typical moment matrix for this waveguide configuration with a resultant sparsity of 91% after a threshold of 1% is applied. Note that the moment matrix in this case is highly asymmetric due to the split testing procedure. The final example of this section involves a hybrid microstrip-dielectric waveguide made up of a single dielectric strip with its top surface fully metallized and resting upon a two-layer substrate structure. The waveguide parameters are Crg = 7rs2 = 8.2, ersl = 12.9, a = w = 0.2mm, b = d2 = 0.33mm, and dl = 0.58mm. The analysis of this geometry requires the simultaneous solution of the system of integral equations given by (4.11), (4.13) and (4.14), with the resulting eigensystem (4.19) including all submatrices. Fig. 4.8 shows the frequency variation of the normalized propagation constant of the dominant mode of this structure, where our results have been compared to those of the mode matching technique presented in [57]. A typical moment matrix of this geometry is shown in Fig. 4.9 after applying a threshold of 1%. At this threshold level, the sparsity of the moment matrix is 99.4%, which is phenomenal for a moment method formulation. The interaction of the conduction current with the top polarization current elements is noticeable at the lower right corner of the matrix.

106 4.00..,.... wi d2 3.50 Ze 3.00 2.50 ~ -- This work o Spectral domain 2.00..00 25.00 50.00 75.00 f [GHz] Figure 4.6: Normalized propagation constant of a microstrip line printed over a double-layer substrate as a function of frequency.

107 * ee* ~ ~~~~~* ~~~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ * ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~~ ~ ~ ~ ~ ~ ~~~~ ~ ~~~~~~~~~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~* ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~- ~~ ~ ~ ~ ~-~ I I I I I I!

108 4.00 | |. | a$ 3.50 t E 5 d 3.00 2.50 2.00 This work 1.50 oE Mode matching 1.00........... 20.00 30.00 40.00 50.00 60.00 70.00 80.00 f [GHz] Figure 4.8: Normalized propagation constant of a hybrid microstrip-dielectric waveguide with a double-layer substrate as a function of frequency.

109 ~\ * * - \ *. * * as' ^, *' * \ * u S. * l I l l Figure 4.9: The structure of moment matrix of Fig. 4.8 after applying a threshold of 1%.

110 4.6 Concluding Remarks In this chapter, we developed an accurate full-wave integral formulation for the study of integrated planar dielectric waveguide structures with printed metallized sections, which are of practical interest in millimeter-wave and submillimeter-wave applications. It was shown that by using multiresolution expansions constructed from dilated and shifted versions of the scaling function and mother wavelet in the moment method solution of the integral equations, highly sparse moment matrices can be generated. This phenomenal sparsity eradicates the traditional limitation of the method of moments imposed by the densely populated matrices of the conventional basis functions. We also underlined the merits of this new formulation by presenting numerical results for a variety of practical two-dimensional dielectric waveguide structures. Except for the spatial sub-domain pulse expansions in the normal direction for dielectric regions, the formulation of this chapter was thoroughly based on a spectraldomain approach. In other words, the unknown quantities in the integral equations were rather the Fourier transforms of volume polarization currents or surface conduction currents. This approach works quite well with significant computational efficiency for planar microwave structures with rectangular cross sections. The bottom line is that the geometry under study must be Fourier-transformable. Therefore, spectral-domain approaches have limited applicability in view of geometrical variety. This aspect of course does not suit CAD-oriented microwave engineering, which prefers general-purpose techniques with wide-range capabilities. Yet, the application of multiresolution expansions to electromagnetic problems is by no means limited to spectral-domain formulations. In the next chapter, we will develop a space-domain

111 wavelet-based formulation for three-dimensional planar structures, which has p)ractically no limits in view of geometrical shapes.

CHAPTER V SPACE-DOMAIN WAVELET ANALYSIS OF THREE-DIMENSIONAL PLANAR MICROWAVE STRUCTURES 5.1 Introduction After the successful application of multiresolution analysis theory to numerical modeling of two-dimensional planar waveguides, the development of a similar spectraldomain formulation for three-dimensional structures was investigated. However, it was found that for this type of geometries, to achieve a high degree of accuracy requires very large numbers of multiresolution basis functions. In chapter 3, we emphasized that, due to the nature of electromagnetic problems, it is essential to use a moment method expansion with significant low-frequency spectral content. The construction of a two-dimensional multiresolution expansion with this property would involve dilations and shiftings of a 2-D scaling function as well as three different types of 2-D wavelets. Although the resulting moment matrices are still highly sparse, and their numerical handling is feasible through the use of special sparse matrix techniques, the amount of computational work for evaluation of moment integrals can easily dominate the entire numerical process. On the other hand, an extensive study of three-dimensional wavelet-based space-domain formulations revealed several 112

113 possibilities to alleviate this computational load. Some of these possibilities were explored in the closing sections of chapter 3. Moreover, in contrast to spectraldomain approaches, where wavelets serve merely as pure mathematical tools, the space-domain formulations provide some physical insight in the light of approximation theory. Due to the zero-average property of wavelets, these basis functions sample the expanded function primarily in the vicinity of discontinuities and abrupt variations. For geometries whose boundaries are separable with respect to the x and y variables, this property prompts an economization of the expansion basis by placing each type of wavelets along those boundaries where the coefficients are expected to be significant. In regions where no strong field variations are suspected, the 2-D scaling functions usually provide a fair approximation of the fields by themselves. In chapter 3, we also investigated the behavior of a wavelet function as an electromagnetic radiator. It was shown that the electric field generated by a wavelet is highly localized and vanishes in the far-field region. In the context of the moment method, this property enables us to predict the weak interactions between basis function in advance. Thus, not only is the numerical evaluation of the corresponding matrix elements spared, but skipping these elements also amounts to considerable savings in the sparse storage schemes. In this chapter, we develop a wavelet-based space-domain formulation of threedimensional planar microwave structures. The primary goal of this chapter is to demonstrate the possibility and existence of a sparse moment method formulation for three-dimensional electromagnetic problems. The use of fast wavelet algorithm in the numerical implementation of this approach dramatically improves the computational efficiency. This improvement becomes much more pronounced in the treatment of 3D metallic structures like microstrip patch antennas, where convergence issues are of

114 critical importance. Not only does this algorithm minimize the amount of numerical integration procedures, but it confines these procedures to a relatively high resolution level, at which the moment integrals have better convergence behavior. Moreover, the entire numerical integration involves only scaling functions, which are more amenable in this respect than wavelets. The study of different 3-D planar structures, which will be outlined in later sections, leads to the conclusion that multiresolution expansions generate highly sparse moment matrices regardless of the dimensionality of the problem. This feature is of course contingent upon the use of a proper complete expansion basis which is in accord with the nature of electromagnetic fields or currents. Since 3-D problems often involve very large matrices whose dimensions may easily exceed the capability of available computers in view of storage and direct numerical handling, the sparsity resulting from multiresolution expansions must be exploited to the full advantage. The problems usually associated with very large sparse linear systems include efficient storage of the matrices as well as efficient numerical solution (or inversion) of the system itself. As for the former problem, numerous storage schemes have been developed which in many cases are tailored to suit specific problems [55]. In the present work, we have used a row-indexed sparse storage scheme, which is a general technique sparse systems have been a subject of intensive research in the past two decades [58]. This work utilizes a pre-conditioned bi-conjugate gradient (BiCG) method, which is usually very efficient for electromagnetic problems [59, 60], and works very effectively with the storage scheme mentioned earlier. This method is discussed in detail in appendix E. Before closing this section, it should be noted that the numerical solution of 3-D

115 eigenvalue problems, such as the determination of resonant frequencies of dielectric resonators, is an extremely difficult computational task especially when the size of the matrices involved is very big. The solution of this type of problems requires the evaluation of the determinant of a very large matrix or computation of the eigenvalues of such a matrix using an appropriate technique such as Lanczos' method [61]. For such problems, we have used an alternative approach in which, instead of solving the eigen-problem itself, we examine the response of the resonating structure to a prescribed excitation. By varying the excitation, we search for one which induces the maximum response by the structure under study. In this way, the matrix eigenvalue problem is replaced with an inhomogeneous linear system of equations, which easily renders itself to the BiCG method. In the following, first we explore some of the special numerical features of spacedomain formulations. It will be seen how the use of Battle-Lemarie functions helps facilitate numerical integrations. Then, we discuss the wavelet-based integral formulation of a dielectric resonator embedded in the free space. Numerical results will be presented to illustrate the merits of using multiresolution expansions for the moment method treatment of this type of structures. To validate the formulation for metallic structures and, in particular, for different substrate geometries (or Green's functions), next we investigate a microstrip patch antenna printed over an infinite grounded substrate. Results concerning the resonant parameters and input impedance of this planar antenna are presented and compared with other methods. Special considerations for the wavelet-based analysis of metallic structures are also emphasized.

11;( 5.2 Preliminary Numerical Considerations Before embarking on the integral formulation of three-dimensional structures, we discuss some basic numerical procedures which are frequently encountered ill spacedomain wavelet formulations. The material presented in this section are of a general nature originating from the mathematical properties of Battle-Lemarie functions and are independent of the specific problem under study. Suppose that the adopted two-dimensional multiresolution expansion consists of dilated and shifted scaling functions and wavelets at the resolution levels m = mo,...,Ima, where mo and ma are the initial and highest resolution levels, respectively. Assume further that the boundary of the problem in the x-y plane consists of a rectangle of dimensions a x b centered at the origin. Then, the fast wavelet algorithm requires the numerical computation of the moment integrals only at the resolution level ma + 1, which we call the base resolution level and denote by mb = ma + 1. According to section 3.5, these interactions have the following form: a/2 Jb/2 Ja/2 -b/2 Kij = /i dx / dy J dx' b dy' J-a/2 J-b/2 J-a/2 J-b/2 x Y xi yl 1Imbnxnny (- 0 d) G(xy | I X' y') mb.,nx(,n.(,- ) >s cl do do3Yl) do do (5.1) where <D(x,y) = 4(z)+(y), and x(zxy I x',y') is a reduced two-dimensional dyadic Green's function, which is obtained through the spatial z-integration of the threedimensional dyadic Green's function for volume polarization currents. In the case of surface conduction currents, the three-dimensional dyadic Green's function is evaluated with both observation and source points located on the plane of the currents.

117 In view of the spectral decomposition of the dyadic Green's function given by equation (A.5) of appendix A, the quadruple integrals of equation (5.1) reduce to double integrals in the spectral variables ( and 77. The resulting integrals then involve twodimensional Fourier transforms of the scaling functions in the following manner: = 1 foo r~o Kij / dl! d77 K (27 )2 dJ- J-o F* { ^4:mbnx_,ny (-T d )Pa( b(y) } g(, I) F{ mbn,. nY, (d' d )Pa(X)Pb(Y) I * do do do' Io do' (5.2) In chapter 3, we saw that the Battle-Lemarie scaling function is a lowpass function whose Fourier transform practically has a compact support of 2wo centered around the origin. We also pointed out that it is necessary to place some basis functions outside the boundaries of the problem. Note that those multiresolution functions which are located close to the boundaries are inevitably truncated. The Fourier transforms of these truncated functions no longer have the fast decay shown in Fig. 3.3, and are capable of causing serious convergence problems. In the following, we reflect upon this issue in more detail. It is convenient to normalize all the space- and spectral-domain quantities in the following manner: X do c27r- ~w (5.3) 2rt For convenience, we drop the bars hereinafter. The following parameter is also introduced:,,. do kodo f:, T> = do - 2w, (5.4) \0 27T''

11I which is in fact the normalized free-space propagation constant. Now, we define the spectral functions: _//2 Wm,n(, a)= I2 m,n(Z) ejddx (5.5) -a/2 which are the Fourier transforms of windowed, dilated and shifted scaling functions, with the window extending over [-a/2, a/2]. Other spectral function WVm,n(77, b) are also defined in a similar fashion corresponding to the y variable. The following property holds: Wm,n(-m, a) = W,-nj, a) = W n (j a). (5.6) In the terminology of Fourier analysis, equation (5.5) is the Fourier transform of a TL (time-limited) function [62]. An interesting theorem regarding this type of functions states that if a function f(x) is TL and has bounded derivatives of order up to n for every | z |< K, then its Fourier transform has the following asymptotic expansion: -f() = \ [f( )e ja/2 - f(- )e-Ja/2 1 [f,( a a e-ja/2 ( / )2 )e - f(- )-i + (-j) [f(n-1)(a)eja/2 _ f(n-1)( _a) -a/2] 0 ( 1 UO00" n 7 f 2l 2a2 J [n+l (5.7) From equation (C.9) of appendix C, it is seen that the cubic spline Battle-Lemarie functions are superpositions of cubic B-spline functions, and hence, possess continuous derivatives of up to the second order. Then, in view of theorem (5.7), it is not difficult

119 to show that the functions WVm,n of equation (5.5) call be written in the following form: k iPk where the amplitudes Ak and phases Ok depend on the indices m and n, and pk are integer powers with 1 < pk < 4 (for more details, see [63]). Equation (5.8) gives an explicit closed-form expression for the Fourier transform of a truncated scaling function. However, the determination of the amplitude and phase of various terms in this equation is relatively involved and is carried out through a computer algorithm [63]. This equation also implies that the Fourier transform of a severely truncated scaling function decays slowly as 1/i. Note that representation (5.8) is not valid for = 0, and at this point the following relation holds: {a/2 Wm,n(0, a)= J 0m,n(x)dx. (5.9) In view of definition (5.5), the moment integrals (5.2) can now be written in the following form: r0oo roo Kij = J dk J d7 Wb,n,.((, a)Wmb,fny (7, b) (., V ) W/Vmb,n, (5, a)Wmbyj(7, b), -o00 -o0 (5.10) with the necessary normalizations having been carried out. Depending on the geometry under study, these integrals may contain surface wave poles or branch point singularities. Appendix B discusses the method of extraction of singularities for the treatment of these types of irregularities. From a numerical point of view, for double integrals like (5.10), it is much easier to treat the singularities in a polar coordinate system defined by - = acosO, r7 = usinO. (5.11)

1 2() The polar forn of equation (5.10) is thenI giveI by Kij = dO adoa Wmb,nx, (a cos 0, a)WHbny, (a sin 0, b) 9(a, 0) Wmb,n, (a cos 0, a),mbny, (a sin 0., b), (5.12) for the case of a single-layer grounded substrate (see appendix A). In the above form, the branch point singularity, if any, is given by Cc = a2- CrC2 = 0, (5.13) where the parameter D defined in (5.4) has replaced the free-space propagation constant. The TM and TE surface wave poles are solutions of the following equations: DTM(a) = Crs c+ scot27r(Csd = 0, DTE(a) = (c-( tan27rCsd = 0, C = /2 - rsD2. (5.14) Note that all the singularities, regardless of poles or branch points, are located in the interval a E [0, V/c7D]. Moreover, using the symmetry properties of the Fourierdomain Green's functions and equation (5.6), the 0-integration in (5.12) can be reduced to the first quadrant of the J-r7 plane. The expansion (5.8) is particularly useful for the analytical evaluation of the tail contribution of slow ly converging moment integrals. A useful technique for the computation of infinite integrals with slow rates of convergence is to split them into different parts in the following manner: too rA Jo Jo

121 + f [f()-fa(f) ] d~ + afa(f)dg (5.15) where f(g) is a slowly decaying function, fa(,) is its asymptotic form, and A is a large number. The second integral on the right hand side of equation (5.15) converges very fast, and the third integral, which is in fact the tail contribution, is evaluated analytically. In view of representation (5.8), the tail contribution of the moment integrals (5.10) can generally be reduced to integrals of the following form: I' W=,n, ( a)Wm,nj 1 Za)d,, E Z (5.16) where A is a large number with A > 2mbwqs. Then, using expansion (5.8), one can rewrite equation (5.16) in the following manner: 00 e (Ok, )-e3 hi = E AE AkA, +P - (5.17) k1 kj AkPk,+Pk But the integrals in equation (5.17) are recognized to be generalized exponential integrals, which are formally defined in the following way: oo e-zt En(z) = dt, RZe(z) > 0. (5.18) Efficient algorithms in the form of series expansions and continued fractions are available for the evaluation of this type of integrals [47]. In particular, the following asymptotic expansion holds [64]: E(z) e-Z { + n(n 1) n(n + )(n + 2) (5.19) En,(z)I- - 23 + (5 19) Z Z Z Z For truncated 2-D scaling functions in which only one of the two 4 functions in x or y is truncated and the other remains intact, it is more convenient to evaluate the tail integrals back in the Cartesian coordinates. In this case, the tail integrals

12I are separable with respect to < and 77 variables. One of the two separated integrals converges fast, and for the other, use is made of equations (5.16)-(5.18). However, in 2-D scaling functions which are located at the corners of the geometry, both q functions in x and y are truncated. For interactions involving two basis functions of this type, the tail integrals are no longer separable, and we have to evaluate them in polar coordinates. In this case, the polar counterpart of equation (5.17) will involve quadruple sums corresponding to the product of the four Wmn functions associated with the shift indices nxi, ni, n., and ny. 5.3 A Wavelet-Based Study of Dielectric Resonators In this section, first we present a three-dimensional integral formulation of a rectangular dielectric resonator embedded in an arbitrary layered substrate configuration. Then, the case of a dielectric resonator immersed in the free space is considered, and explicit expressions for the reduced dyadic Green's function of this geometry are derived. The numerical implementation of the method of moments is described in detail, and numerical results are presented and verified by comparison with other methods. 5.3.1 Integral Formulation Consider a rectangular dielectric resonator of dimensions a x b x h with a uniform relative permittivity of Erg, which is embedded in a planar layered background structure. The geometry of the problem is shown in Fig. 5.1. The rectangular dielectric region can be replaced with a volume polarization current defined in the following

12:3 Figure 5.1: Geometry of a rectangular dielectric resonator. way: Jp(r) = jkoYo & E(r) pv(r), (5.20) where Er = Erg - Er is the index contrast, and pv(r) = pa,()pb(y)ph(Z) is the characteristic function of volume V occupied by the resonator. The cover medium is assumed to be free the space (erC = 1). The electric field radiated by this polarization current is then given by: E(r) = -jkoZo J JV Ge(r r).Jp(r')dv' + E'(r), (5.21) where Ei(r) is the incident electric field. Note that when the electric field is evaluated inside the source region, i.e. r V/, the electric field dyadic Green's function exhibits a singularity at the source point. In view of equations (5.20), (5.21) and (A.2) of appendix A, the following integral equation for the polarization current is obtained: (I+ Zr L j.4p(r) = 8ke pv(r) J JI Ge(r | r'). J (r) dv

121 + jkoYo 6,E'(r)pvJ(r), (5.22) where L = -z is the source dyadic (see appendix A). The numerical solution of equation (5.22) through the method of moments requires the expansion of the polarization current in a suitable basis. To this end, a two-dimensional multiresolution expansion in the x-y plane in conjunction with a sub-domain pulse expansion along the z-axis is utilized. The construction of 2-D multiresolution basis functions was discussed in detail in chapter 3. Using the shorthand notation introduced in this chapter, one can write: J,(r) = aip Fi( ) pp(Z), (5.23) t p Qdo' do where Fi(x,y) stands for a 2-D multiresolution basis function, and do is a characteristic length. The range of the shift indices nX and ni for each resolution level m of expansion (5.23) depends on the size of the cross section of the resonator and is determined with the aid of the wavelet localization lattice. Since in practice dielectric resonators are made of very dense dielectric materials with large permittivities, we choose for the characteristic length the material wavelength defined as do = -. (5.24) In the normal (z) direction, the dielectric region is discretized uniformly into Nh, sub-layers of thickness A = h/Nh. The integrations with respect to variables z and z' are carried out analytically in a similar fashion as in section 4.4 of the preceding chapter. We introduce a reduced planar Green's function in the following manner: -pq(x,y I x',y) = rko / dz f dz' p(z) Ge(r I') p(z'). (5.25)

1 -')) Then, the moment integrals can be written in the following form: ZiplJq Piplj - I L A) A6 pq (5.26) where ra/2 rb/2 Ta/2 b/2 Pipjq = /2 dx b/2 dx' J/ dy' J-a/2 J-b/2 J-a/2 J-b/2 Fi(- ) pq(X y I', y') Fj(' ). (5.27) do do do' do As pointed out in the preceding section, because of using fast wavelet algorithm, the actual numerical integration is performed only for 2-D scaling functions at the base resolution level. Hence, the integrals of equation (5.27) are reduced to those of equation (5.1). The related spatial quadruple integrals are then converted into spectral double integrals of the form given by equation (5.10). These latter integrals involve the Fourier transform of the reduced dyadic Green's function defined in (5.25). For the case of a rectangular dielectric resonator embedded in the free space, this Fourier transform is given by the following expressions: __ [D+(',,,rc. )+D -(C,,rc)] (- 1-CA) r=0, r(,77) 2= -2 D qrc) c (ec -) (1 -e- ), r # 0, (5.28) where the dyadic functions D are defined in equation (A.11), r = p - q, and the upper and lower signs correspond to positive and negative values of r, respectively. In the beginning of this chapter, mention was made of an alternative approach for the solution of eigenvalue problems like the one at hand, which is especially suitable when large sparse linear systems are involved. The incident field Ei(r) retained in the

126 integral equation (5.22) serves for this very purpose. The numerical solution of this integral equation yields the unknown amplitude coefficients a,p of the polarization current (or electric field) induced by this excitation inside the dielectric resonator. Applying the Galerkin technique to the discretized integral equation leads to the following system of matrix equations: [Zipljq ] = [a bi,], (5.29) where the vector on the right hand side of (5.29) is the excitation vector, with its elements given by,a/2 b/2 h o X bip= -jkoYo Ser / dx dy dzE (r)Fi(-, d-)pp(z). (5.30) J-a/2 -b/2 JO do do 5.3.2 Numerical Procedure and Results To determine the resonant frequency of the dielectric resonator, the structure is excited by a normally incident plane wave with a prescribed frequency, and the electric field induced inside the resonator is examined. When the frequency of the incident plane wave matches the resonant frequency of one of the resonating characteristic modes of the resonator, the induced response reaches a maximum. Therefore, the excitation frequency is varied over a search range, and the local maxima of the induced electric field corresponding to various resonant modes are determined. The incident plane wave has the following general form: Ei(r) = e-ko'.r, (5.31) where ko = w/c, c is the speed of light in the free space, and iu and v are two unit vectors. The divergence equation requires that u.v = 0. (5.32)

127 Note that different polarizations of the incident field excite different modes of the structure. This feature of our alternative approach to the resonator problem is very useful for several reasons. First, the study of degenerate modes in a direct solution of the eigen-problem usually poses numerical difficulties due to instability of the determinants at frequencies close to the resonant frequencies of such modes. Second, by bringing the incident field into the scene, one can also study the problem of excitation and mode coupling in dielectric resonators at no extra effort. It is also important to note that in this alternative approach, the field distribution inside the structure is determined simultaneously with the resonant frequency. As opposed to the spectral-domain formulation of the preceding chapter, in a space-domain approach, the selection of basis functions is based on physical considerations. Once given the geometry and electrical properties of the resonator, all dimensions are first normalized to the material wavelength according to equation (5.24). The structure is partitioned uniformly into smaller sub-layers along the normal direction, with the number of z-cells depending on the normalized thickness (or height) of the resonator. In most cases, 2-4 pulse functions provide very good results. For each sub-layer, a 2-D multiresolution expansion is constructed in the transverse variables x and y. For a simple geometry like the rectangular resonator at hand, which has a uniform cross section along the z-axis, the same multiresolution expansion but with different unknown coefficients can be used for all z-cells. However, for more complex structures with non-uniform cross sections along the z-axis, different multiresolution expansions for each sub-layer may be needed. The initial resolution level mo for the 2-D multiresolution expansion is chosen based on the normalized dimensions of the cross section of the resonator. It is important to remember that multiresolution expansions are very flexible in that if a

128 certain selected basis does not provide the required accuracy, one can easily enlhance it by adding wavelets of higher resolutions. In this process, the original basis functions remain intact and the related computations need not be repeated anew. This is a significant advantage over the conventional MoM schemes, in which the entire expansion basis is changed in the course of fine tuning. We usually choose the initial resolution level such that a few 2-D scaling functions are fitted inside the boundary of the structure. These scaling functions at the initial resolution level provide an crude approximation for the polarization current. Then, 2-D wavelets at this initial level and subsequent higher levels are added to the expansion basis. As described in section 3.5 of chapter 3, the placement of 2-D wavelets can be done in a selective manner to spare those basis functions which are expected to have negligible expansion coefficients. In other words, wavelets of types Tv and Ah are mostly placed in the vicinity of the vertical boundaries at x = ~a/2 and horizontal boundaries at y = ~b/2, respectively. The diagonal wavelets Ad are concentrated primarily at the corners of the rectangular cross section. The highest resolution level ma is chosen again with physical considerations in mind. This resolution must be compatible with the smallest physical scales involved in the field or current distributions. On the other hand, for efficient numerical evaluation of the moment integrals, the base resolution level, namely mb = ma + 1, must be high enough so that at least a few 2-D scaling functions stay away from truncation by the boundaries. The number of multiresolution basis functions at each resolution level (including also the base level) can be minimized outside the resonator by considering the fast decay of these functions. The wavelet localization lattice of Fig. 3.6 is especially useful for this purpose. The case of a rectangular dielectric resonator embedded in the free space has been investigated. Due to the use the free-space dyadic Green's function, the moment inte

129 grals are free from surface wave poles; however, they have a branch point singularity of order 1/2 at a = D. To remove this singularity in the course of numerical integration, the method of extraction of singularities have been employed (see appendix B). The moment integrals involving 2-D scaling functions at the base resolution level are computed numerically using a Gauss-Legendre quadrature. The integration is performed in the spectral domain and in polar coordinates a and 0 as shown in equation (5.12). Due to the spatial integration along the z direction, the integrals of (5.12) converge relatively fast even in the case of severely truncated basis functions such as those located at the corners of the rectangular domain. Therefore, the analytical evaluation of the tail contribution discussed in section 5.2 is not necessary in this problem. After all 4>-<> interactions at the base resolution level have been evaluated numerically, the elements of the moment (impedance) matrix [Z] are computed using the fast wavelet algorithm (FWA). Note that the numerical integration at the base resolution level must include all the <-4 interactions between different sub-layers (pulse basis function in z). However, for the free-space Green's function in the present problem, these interactions are a function of the relative distance of the pulse functions, i.e. r = p - q, rather than their absolute position p and q. Hence, the FWA loop is repeated for different values of r = 0,1,2,..., yielding the elements of various interpulse submatrices in the moment matrix. In the FWA computation of the moment interactions, the perfect symmetry of the Battle-Lemarie functions are fully exploited in conjunction with the symmetry of the free-space dyadic Green's function. The former symmetry with respect to the shift indices is in the following form: 4.m,n(X) = <m,-n(-X),

130 V.m,n(X) = _m-_-l(-x). (5.33) Due to these symmetries, many of the moment matrix elements are either equal or differ only in a sign. Equation (3.25) of chapter 3 shows how the various elements of the [Z] matrix are computed recursively from the base-level 4mb -<mb interactions. This procedure requires the evaluation of all possible 4)- interactions at each resolution level involved. In practice, the FWA loop starts with the highest resolution level ma = mb - 1, where first all ma,-ma interactions and then intra-level ma - ma interactions are calculated. Next, the inter-level,ma -m interactions with mo < m < ma, and the Tma4-mo interactions are computed. This completes all the interactions which involve the resolution level ma. The FWA loop then proceeds to the resolution level ma - 1, and the cycle is repeated again at this resolution. The loop continues downward to the initial resolution level mo. Note that in view of reciprocity, there is no need to compute the inter-level Tmr-m m2 interactions with mi < m2. It is important to note that, when using equation (3.25), the shift indices kx,, k,,, kxj and kyj at each pair of resolution levels mi and mj have finite ranges, which are limited by the boundaries of the geometry. Moreover, the coefficient of each term in the related quadruple sum is the product of four hk or gk coefficients, and for a vast majority of terms, this product becomes negligibly small. Such terms are therefore skipped in the course of the quadruple summation. Finally, in view of the far-field properties of wavelets as discussed in detail in section 3.6, a large number of moment matrix elements correspond to far-field interactions between 2-D wavelet basis functions. These elements are identified in advance and subsequently skipped during the FWA loop. After the entire impedance matrix [Z] has been determined, it is scanned for

131 thresholding. In fact, in addition to all far-field interactions, which are skipped in the FWA loop and replaced with zero elements, it turns out that a large number of other elements of [Z] still have magnitudes several orders smaller than its largest entry. Such elements can be discarded by establishing a threshold with respect to the largest entry of the moment matrix. The fact that the accuracy of the solution is not much affected by this threshold process is one of the most interesting features of multiresolution expansions, and it stems from the mathematical properties of the multiresolution analysis [51] (See also section 3.1). The resulting moment matrix after thresholding is very sparsely populated. This sparsity enables us to avoid the storage of the entire matrix, which can be extremely large in practice. Instead, only the non-zero elements are stored using an efficient scheme such as the row-indexed sparse storage method [47]. After computing the excitation vector, the sparse linear system (5.29) is solved iteratively using a pre-conditioned BiCG algorithm, which is described in appendix E. This algorithm is very efficient for this problem, and in most cases less than 20 iterations are sufficient for convergence. By the numerical solution of linear system (5.29), the amplitude vector [a] which contains the unknown coefficients of expansion (5.23) is determined. Then, the electric field at any point inside the dielectric resonator can be computed in the following way: E(r) = j C E aip Fi(d, ) PP(Z), r E V. (5.34) J lkoo J r, p do do To determine the resonant frequency of the structure, the frequency of the incident plane wave is varied. In general, if the field distribution of a certain mode is roughly known, then the electric field can be probed at an appropriate point inside the resonator (e.g. at its center) to find the maximum response to the excitation. If such knowledge is not available, then one can compute the stored energy inside the res

132 Geometry Mode fres (This work) fres (Ref. [65]) TMji 5.532 GHz 5.651 GHz ~....A/l l - TM6 6.012 GHz 6.104 GHz Er9 a TEz 1 6.378 GHz 6.575 GHz Table 5.1: Resonant frequencies of different modes of a rectangular dielectric resonator onator for each frequency of the excitation. In view of equation (5.34), it is not difficult to show that the stored electric field energy is given by the following expression: Fe = 2 JJv 6rg60 I E(r)|2 dv ErgoE 2 C z J-/2 dxJb/2 dyF, o)F (d, o) SC= E > I: aip. ajP dx y dy F).( -)Fj( ) 2(koYo,),, -a /2 -b/2 a0' da do' do (5.35) To validate the formulation developed above, we consider a rectangular dielectric resonator of dimensions 10mmx8mmx5mm with a relative permittivity of erg = 20, which is embedded in the free space (crc = 1). Table 5.1 shows the computed resonant frequencies of the first three modes of this resonator. The results have been compared to those based on Marcatili's approximation [65], and a good agreement is observed. The initial resolution is taken to be mo = 2. Only two resolution levels (m = 2,3) and 2-4 pulse functions along the normal direction are sufficient to provide very accurate results. To better understand how the resonant frequency is determined by

133 varying the excitation frequency, Fig. 5.2 shows the variation of tile imagnitude of electric field at the center of the resonator (x = y = 0, = //2) for the dominant TMi1 mode as a function of the frequency of the incident field. The peak of the curve identifies the resonant frequency of the dielectric resonator. Note that for the dominant mode considered above, the electric field always has a maximum at the center of the resonator. The resulting moment matrix, as expected, is highly sparse in the sense that a very large number of its entries are quite small in magnitude when compared to the largest entry. Fig. 5.3 shows the structure of the moment matrix after applying a threshold of 1%. In this case, the expansion basis consists of 2 pulse functions and a total of 231 2-D multiresolution expansion functions, and the sparsity of the moment matrix is 99.16%. Finally, Fig. 5.4 shows the variation of the resonant frequencies of the first three modes of the dielectric resonator considered above as a function of the aspect ratio h/a. 5.4 A Wavelet-Based Study of Microstrip Patch Antennas To further validate the wavelet-based space-domain MoM formulation developed in this chapter, we have applied it to 3-D metallic structures with infinite substrate geometries. This type of geometries serve a dual purpose. On the first place, the special problems associated with metallic structures such as s enforcement of boundary conditions for conduction currents and difficulties in the convergence of moment integrals are investigated. On the other hand, a different type of Green's function other than that of the free space is being examined. In this section, we study the geometry of a rectangular microstrip patch antenna. After formulating the integral equation, first we regard the patch as a planar resonator. The resonance character

134 20.0 t 10.0.000E+00.OOOE+00 * * * *1.... i.... I.... X.... i.... |....!.... 3.50 4.00 4.50 5.00 5.50 6.00 6.50 7.00 7.50 f [GHz] Figure 5.2: Variation of E-field magnitude at the center of the resonator with the frequency of excitation field.

135 1 4 I"- t- - "' 1g~\^ E m'. *'t1*'*, >. %'.9' 1.X X *I *..^ _~'I _ ~ a "Z''?',. a thresho. of ^ I'~e V. " L.. L*.9 %o' 4 * * am.. eqbe I I I I a threshold of 1%. S ~ ~.44h~~S~j' I Ib%'',*.4. *-b 4b lb ~ 1 %.4.4

136 8.00.......... 7.00 N 6.00 - 01 TEz111. 0 0~s~oo 0 0 ~ TMY1, 5.00 El TMX111 3.00 C Marcatili's method 0a This work 2.00..........30.40.50.60.70.80.90 h/a Figure 5.4: Variation of resonant frequencies of dominant modes of a dielectric resonator as a function of the aspect ratio h/a.

i:]T istics of the patch are determined in a similar manner as the dielectric resonator of the preceding section. For this purpose, a plane wave incident field is considered as the excitation. Next, we treat the patch as a radiating structure and determine its input impedance. To this end, we need to model the feed element in a realistic way. Microstrip patch antennas are usually fed by two conventional methods: (1) through a coaxial probe connected to the patch from within the dielectric substrate, or (2) via a microstrip line printed on the same substrate and connected to one edge of the patch [66]. Other more complex feed structures have also been proposed and studied [67, 68]. In the present work, we consider the first method, namely a coaxial feed. To compute the antenna input impedance, the effect of the feed element is modeled as an impressed current. The incident electric field is then assumed to be radiated by this impressed current. After the unknown patch current has been determined, by establishing proper reference voltages and currents, the antenna input impedance is computed. 5.4.1 Integral Formulation Consider a rectangular metallic patch of dimensions a x b, which has been printed over an infinite conductor-backed dielectric substrate of thickness d and relative permittivity,r. The patch is located on the x-y plane and the perfectly conducting ground plane is located at z = -d as shown in Fig. 5.5. It is postulated that the incident field excites induced current on the surface of the patch, which in turn radiate an scattered field. The total electric field can be expressed in the following form: E(r) = -jkoZo J Ge(x, y, zx',y', O). J(x',y')d'dy' + E'(r), (5.36)

Ei _ I -- a 13$ d T firs Figure 5.5: Geometry of a rectangular microstrip patch printed over a grounded substrate. where S is the surface of the patch, and Jc(x, y) is the planar conduction current sustained on this surface. The integral equation for the unknown conduction current is derived in view of the boundary condition on the surface of the metallic patch. Assuming a perfectly conducting patch, this condition requires that the tangential components of the electric field vanish everywhere on S. In the real physical problem, the patch has a finite nonzero thickness of Ah, and the planar current is mainly concentrated at the bottom surface of the patch. For numerical purposes which will be explained later, it is convenient to test the field on the top surface of the patch taking into account its nonzero thickness. Then, the boundary condition on S can be written in the following way: Et(x, y, Sh) p,(x, y) = O, t = x, y, (5.37) where ps(x,y) = pa(x)pb(y) is the characteristic function of surface S. In view of

139 equations (5.36) and (5.37), the following integral equation for the unknown patchi current is obtained: p,(x, y)t. {-jkoZo ffs Ge(, y, 6h I x', y', O). J(x', y')dx'dy' +Ei(x, y,h)} = 0, t = xy, Ah - 0. (5.38) The enforcement of compactness of polarization and conduction currents was emphasized in chapter 3 when discussing the construction of multiresolution expansions for moment method treatment of electromagnetic problems. In chapter 4, we described how to implement this condition in the spectral domain. It was seen that due to the nature of the integral equations for polarization currents, the compactness of this type of currents are automatically guaranteed. However, for conduction currents, we needed to enforce this condition explicitly. An inspection of equations (5.22) and (5.38) leads to the same conclusion for a space-domain formulation. In other words, the integral equation (5.22) implicitly forces the polarization current (and hence its multiresolution expansion) to vanish outside the dielectric region. However, the integral equation (5.38) does not impose any restriction on the compactness of the patch current. This condition is enforced in the following explicit form: [1 - ps(x, )] JC(x,y) = 0. (5.39) To determine the unknown patch current, equations (5.38) and (5.39) must be solved simultaneously. For the moment method solution of these equations, the patch current is expanded in a two-dimensional multiresolution expansion as discussed in

140 chapter 3. Using the shorthand notation of this chapter, one can write J(x, y) = aFi( - ), (5.40) do' do where the characteristic length do is chosen to be an effective wavelength defined in the following manner: Ao e~rs + 1 (1 o-=, / E, =- (5.41) io n tene e set2 equa By inserting expansion (5.40), equations (5.38) and (5.39) are discretized. The application of the Galerkin technique to these discretized equations in the usual way obviously results in an overdetermined linear system with the number of equations bigger than the number of unknowns. Another major problem is that testing equation (5.39) with those multiresolution basis functions which are well inside the boundaries and are not truncated leads to trivial or redundant linear equations. Therefore, the following testing scheme is adopted: * Test integral equation (5.38) with interior basis functions, whose centers are located on the surface of the patch. * Test the additional equation (5.39) with exterior basis functions, whose centers fall outside the surface of the patch. This testing scheme produces a consistent linear system in the following form: E [Zi, ajx + ZjYa] = bj (mi,n,ny-) S, E [Zyaj + Z^Yay] = biy, (mi,n,nbi) E S, E Sijaj~ = O, (m,,?x,,n^) ^5, =X Y 0, (mi, nx,ny,) ~ Si Sijajy = 0, (mi, ni,ny) S, (5.42)

141 where Aa/2 /b/2 fa/2 rb/2 Z0 = jkoZo /dx b/ d y / dx' I dy' J-a/2 J-b/2 J-a/2 -b/2 Fj ( Y ), - ~ (To) I x',y)F(L ) r o~= X, y, do do do' do (5.43) and fa/2 lb/2 X 7' I/ J-a/2 J-b/2 do do do' do The Green's function 5(x, y I x', y') in equation (5.43) is the dyadic Green's function of a single-layer conductor-backed substrate evaluated at z = 6h and z' = O. From section A.3 of Appendix A, the Fourier transform of the transverse components of this Green's function are given by the following expressions: 1 [ 2 DA ] e DA-_ t2 DA (4DrM = 1 [1- 1I e-c^ DTE L k DTM J (5.45) where DTE and DTM are defined by equation (A.19) and DA = Cc + ( tanh (Cd. (5.46) The excitation vector in the system of (5.42) is given by /2 dx dy y,/2 = b;, = dx d yE/2x~.,y,6 — ), = -x,y. (5.47) J-a/2 J-b/2 do do

142 5.4.2 The Microstrip Patch as a Resonating Structure Determination of the resonant frequency of a microstrip patch is based on the same methodology for the dielectric resonator of the preceding section. The patch is regarded as a 3-D resonating structure, and its response to a prescribed plane wave incident field is examined. The excitation has the same form as equation (5.31). Note that the plane wave must be polarized such that it has a tangential component in the x-y plane. Otherwise, it would give rise to a zero excitation vector, and the BiCG method would not work. The numerical procedures for the problem of the patch resonator are quite similar to those of the dielectric resonator. The major difference is that the former involves a surface integral equation, while the latter deals with a volume integral equation and hence requires more intensive computational work. Also, the Green's functions of the two problems are different, each having its own singularities and convergence properties. Apart from this dissimilarity in Green's functions, which manifests itself only in the computation of the base-level moment integrals, the implementation of the fast wavelet algorithm for the two problems is done exactly in the same manner. In fact, all the numerical aspects of the two problems as related to multiresolution expansions are identical. This includes the construction and selection of the 2-D multiresolution basis functions, all possible shortcuts in the calculation of moment interactions, and the final thresholding of moment matrices. The storage of matrices using sparse schemes and the implementation of the BiCG algorithm are also done as before. It would be appropriate here to discuss some of the special problems associated with the numerical integration of moment integrals for the patch problem. First of

143 all, the dyadic Green's function of equation (5.45) contains surface wave poles, which are indeed the eigenvalues of the infinite substrate geometry. These poles, however. are removable, and are treated using the method of extraction of singularities, which is described in detail in appendix B. The main difficulty in the numerical evaluation of the moment integrals in this problem is the slow convergence of some of these integrals which involve truncated basis functions. This problem stems from the fact that the fields are tested right at the location of the current source, or in other words. where the Green's function exhibits a singularity. The small non-zero patch thickness Sh, which we have retained in our formulation serves this very purpose. Note that this is a pure mathematical anomaly and does not emerge in the underlying physical problem. The remedy for speeding up the convergence of these ill-behaved moment integrals was prescribed in section 5.2. The example to be studied in this section consists of rectangular patch of length a = 5.95cm and varying width b which is printed over a dielectric substrate of thickness d = 0.159cm and relative permittivity crs = 2.62. For this structure, an xpolarized plane wave propagating along the negative z-axis is incident upon the plane of the patch. Fig. 5.6 shows the variation of the resonant frequency of the patch as a function of the patch width b. The results have been compared to those of the regular moment method [69], and a very good agreement is observed. Figs. 5.7 and 5.8 show the patch current distribution in the x and y directions, respectively. As expected, the current has a quasi-sinusoidal form in the direction of the incident field polarization, and is almost uniform in the orthogonal direction except at the edges where a singularity is observed. Fig. 5.9 shows the structure of the moment matrix after performing a threshold of 1%. The asymmetry of the moment matrix due to the special testing scheme used is easily noticeable. In this case, the initial resolution

144 1550................| & 1525 CZ E This work 1500 - o ) 1475. _ —-—..Conventional MoM tL rn This work.00 2.50 5.00 7.50 10.00 12.50 15.00 17.50 20.00 b [cm] Figure 5.6: Variation of resonant frequency of a microstrip patch antenna as a function of patch width.

145 is mo = 2, and two resolution levels m = 2,3 with a total of 235 2-D multiresolution expansion functions have been used. The sparsity of the moment matrix is 97.3%. 5.4.3 The Microstrip Patch as a Radiating Element The antenna characteristics of the microstrip patch generally depend on its feed mechanism. As mentioned earlier, two popular ways of feeding a patch antenna are through a microstrip line or a coaxial probe. The accurate modeling of the antenna feed is usually a very difficult task and inevitably involves some type of approximation. To this end, the effect of the feed structure is often modeled as an impressed current denoted by Jf (r). Depending on the type of modeling and the geometry under study, this current may consist of a simple known localized source [70, 71] or may involve a complex unknown distribution, which has to be determined simultaneously with the unknown patch current [67, 72, 73]. In any case, the impressed feed current is responsible for generating the incident field Ei(r) in the integral equation (5.38). The resulting incident field is then given by the following expression: E'(r) = -jkZ J Ge(T I r').Jf(r')dv', (5.48) where volume Vf is the domain of the feed current. In the implementation of the method of moments, it is seen that the feed mechanism only affects the excitation vector, and the impedance matrix of a probe-fed patch is exactly the same as that of the patch resonator with a plane wave incident field, which was studied earlier. In other words, the numerical solution of the boundary value problem again reduces to the linear system of equation (5.42). This is of course true only if the feed structure itself is not discretized. The computation of the excitation vector for a given impressed current is by far more difficult than the case of a given incident electric field. It is

146 1.25 1.00 O 7 0.7 5 0 K<.500 y=O.250.000E+00 -.50 -.40 -.30 -.20 -.10.00.10.20.30.40.50 2x/a Figure 5.7: Patch current distribution along the direction of current. 3.00...............g......... 2.00 x=o X-O 1.00.000E+00........................................ -.50 -.40 -.30 -.20 -.10.00.10.20.30.40.50 2y/b Figure 5.8: Patch current distribution normal to the direction of current.

147 IK-".':,; i...... L A4 ^.^;r.-..,~,,.......' XW';'" "*' \',4% % \,: v ' Jr 3'y- **. \*t ~. " ^' "i\ \'' ~' \ -,. % ~ ~ ^ y i I.%'"~'\ \, s,,: \ \, ~ Figure 5.9: The structure of moment matrix of a microstrip patch after performing a threshold of 1%.

148 often more convenient to make use of the reciprocity theorem [74] and interchange the role of the multiresolution test function and the feed current in equation (5.47). Assuming a test function of the form Fi(x/do, y/do), - =, y, one can write _a/2 lb/2 x yE b, a/2 J b/ F2 (, Y).E(x,y,Sh)dxdy = I V Jf,(r) Ei,(r)dv /* /* /* r Ja/2 rb/2 =-jko~ J J Jvd J-a/2 J-b/2 Jf (r). Ge(, y, z I', y', Ah). ftFi(d, o-)dx'dy' = x,y, (5.49) where Ei,(r) is the electric field radiated by the test function. The main reason for this manipulation is that when the current source is embedded inside the substrate (as is the case for a probe-fed patch), the dyadic Green's function may become quite complicated, while the Green's function of the last integral in equation (5.49) is due to a source located on the interface (or in the cover medium). As mentioned above, various simple or complicated models for a coaxial probe feed have been developed. Since the primary objective of the present work is rather to validate the wavelet formulation for 3-D metallic structures, we adopt the simple model of reference [70]. This model usually provides satisfactory results for thin substrates with thicknesses limited to a few hundredths of the free-space wavelength. It assumes a uniform current distribution over the coaxial probe and ignores the effect of the coaxial aperture. The accurate analysis of thicker substrates with big coaxial

149 microstrip patch y x coaxial feed Figure 5.10: Geometry of a coaxial-probe-fed microstrip patch antenna. apertures of course require complicated models such as those proposed by [72, 73]. Fig. 5.10 shows the geometry of a coaxial probe-fed patch antenna. The simple model of [70] replaces the coaxial probe with a z-directed current filament in the following form: Jf(r) = zS(x - f)S(y - yf), -d < z < 0, (5.50) which is connected to the patch at point (xf, yf). Inserting (5.50) in equation (5.49) leads to the following expression: /a/2 rb/2 rO _ j Zl y..= k -j coZo j/dxz dy' dzG (xf,yf, z IxyS h)Fi(-,) -a/2 J-b/2 J-d ezpy do' do) 1, = x,y, (5.51) where the dyadic Green's functions is evaluated with the observation point inside the substrate (See section A.3 of appendix A). In view of numerical integration, the integrals of equation (5.51) are very similar to those of the impedance matrix given by equation (5.43), and one can therefore use the fast wavelet algorithm to facilitate the computation of the bi coefficients. The numerical evaluation of these coefficients would therefore be limited to the base resolution level mb only. These latter integrals

150 can be written in the following form: jkoZo doc, 0 bi = - I dJ J d77 dz. A(27) Jo J-oo J-d Goezi (, 1, z I a) e-J(+Y) W blX (, a) W'^n b (,7 & (), - = xy, (5.52) Having found the patch current distribution, one can compute the electric field everywhere in the structure. By defining an input voltage and current, the input impedance of the patch antenna can be determined in the following manner: Vi n f~d z. E(xf,yf, z) dz Zin = - i - (5 53) Iin Io where the input current is assumed to be Io = 1A, and the definition of the input voltage is based on the total electric field inside the substrate, which is the sum of the fields produced by both patch and feed currents. In view of equations (5.36),(5.40) and (5.51), it is not difficult to show that J Ez(xf, f, z)dz = ai.bi, (554) where E8(r) denotes the scattered field radiated by the patch current. On the other hand, the contribution of the feed current to the input voltage can be expressed in the following form: E'(xf,yf,z)dz = -jkoZo / Gz(xf, yf, z xf.yf. z)dzdz. -d J-d J -d (5.55).=III where Ge (r I r') is the dyadic Green's function due to a source embedded inside the substrate, with the observation point being also located inside the substrate. The

151 Fourier transform of the zz component of this Green's function is given by cosh +d)[ C cosh s sinh'] - cosh,d Erc cosh (z - Ers ( sih z> z', (5.56) where the various quantities are defined in section A.3 of appendix A. The contribution of the excitation field to the input impedance is usually very small around the resonance of the patch antenna. However, it becomes important at low frequencies [75]. Note that references [70, 71] have ignored this effect in the computation of the input impedance, and to account for the self-inductance of the probe, they add an additional reactance term jXp given by (see [66]) Xp = Zo tan(V/kod). (5.57) The first numerical example of this section consists of a rectangular patch antenna with a = 13.97cm, b = 20.45cm, ~rc = 1, Crs = 2.59, and d = 0.1588cm. The location of the coaxial feed is at xf = -6.35cm and yf = 0. The resonant frequency of the patch is fr = 658MHz. Fig. 5.11 shows the variation of the, patch input impedance (normalized to 509) as a function of frequency. This figure also compares our results with those of [70], and very good agreement is observed. The next example is that of a rectangular patch antenna with a = 7.62cm, b = 11.43cm, r = 1,,rs = 2.64, and d = 0.159cm. The location of the coaxial feed in this case is at zf = -3.05cm and yf = 0.385cm, and the antenna resonates at fr = 1189MHz. Fig. 5.12 shows the variation of the normalized input impedance

152 - 5Conventional MoM * This work Figure 5.11: Variation of patch input impedance as a function of frequency. Increment: 5MHz (increasing clockwise).

153 of this antenna as a function of frequency, where the results have been compared to those of [71]. 5.5 Concluding Remarks In this chapter, we developed a general wavelet-based space-domain formulation for the analysis of three-dimensional planar microwave structures. To validate the formulation, it was applied to two different 3-D geometries: (1) a rectangular dielectric resonator embedded in the free space, and (2) a rectangular microstrip patch antenna printed over an infinite grounded substrate. The study of these two structures enables us to investigate various computational aspects of our wavelet-based formulation in view of accuracy, efficiency and convergence issues. In both cases, it was demonstrated that the use of 2-D multiresolution expansions with the method of moments leads to the generation of very sparsely populated moment matrices. Although the sizes of the matrices involved are quite big, due to the resulting sparsities, one can take advantage of highly efficient numerical tools such as the bi-conjugate gradient method. Moreover, the use of special sparse storage schemes makes the numerical treatment of large-scale problems involving thousands of unknowns practically feasible. In addition to the sparsity of moment matrices, the combination of fast wavelet algorithm and the special features of the Battle-Lemarie multiresolution analysis make a dramatic improvement in the computational efficiency of our wavelet analysis. The overall outcome is a very accurate, efficient and versatile, rigorous full-wave technique, which can easily compete with heavily intensive differential-based approaches while maintaining all the advantages of integral-based formulations. Before closing this chapter, it is worth noting that the formulations presented in

154 ---- Conventional MoM * This work Figure 5.12: Variation of patch input impedance as a function of frequency. Increment: 10MHz (increasing clockwise).

155 this work are the first steps in the application of wavelet theory to the numerical modeling of 3-D electromagnetic problems. This is of course a starting point for the development of a general-purpose wavelet-based MoM formulation for more complex geometries of arbitrary shapes and with possible material inhomogeneities. The methodology described in this chapter is very general in nature, and easily renders itself to geometrical extensions. Moreover, although the Battle-Lemarie expansions were used throughout this work, the basic formulation is independent of the type of the multiresolution analysis used. The only limiting aspect of the procedures discussed in this chapter is due to the use of semi-closed-form expressions for the Fourier transforms of windowed basis functions, which drastically speeded up the numerical evaluation of base-level moment integrals. These expressions are based on the special mathematical properties of the Battle-Lemarie functions, and do not apply to other types of multiresolution expansions.

CHAPTER VI CONCLUSION 6.1 Summary of Achievements This dissertation is concerned with the numerical modeling of planar passive microwave structures using integral-based techniques. Although for many electromagnetic problems, especially those involving open geometries, the integral equation approach provides the most accurate formulation of the underlying boundary value problem, its practical use has been traditionally limited to small-scale or mediumscale problems. This limitation is due to the fact that all conventional types of basis functions regularly used in the method of moments generate full moment matrices. The analysis of large and complex structures is usually accompanied by matrices of prohibitively large sizes, whose numerical handling may easily exceed the capability of the available computing resources. The present work proposes two different approaches to eliminate this major limitation. The first approach involves the derivation of a new integral equation for planar dielectric structures. This technique, which has been named the integral transform technique, employs higher-order boundary conditions in conjunction with Taylor expansions of the fields in an appropriate integral transform domain. In this way, an equivalent system of integral equations with reduced dimensionality are obtained. By 156

157 using a proper expansion basis such as the set of orthogonal Hermite-Gauss functions, it is possible to achieve very accurate results with only a few expansion terms. This technique was validated for a variety of planar dielectric waveguides. The primary limitation of this highly efficient approach is its critical dependence on the geometry of the structure as well as availability of amicable Green's functions with certain properties. The second approach to improve the computational efficiency of the method of moments exploits the newly developed orthonormal wavelet expansions. In contrast to the first technique, this approach is very general in nature and can be generalized to any arbitrary geometry or arbitrary material inhomogeneity. In this work, first we presented a spectral-domain wavelet-based formulation for integrated planar waveguides with arbitrary numbers of dielectric and metallic sections. Extensive numerical results were given for various two-dimensional structures, and the accuracy of the results was verified by comparison with other techniques. Then, acknowledging the limitations of spectral-domain formulations in view of geometrical variety, we developed a wavelet-based space-domain formulation for three-dimensional planar microwave structures. This latter formulation, of course, can also be specialized to two-dimensional geometries. The 3-D space-domain MoM formulation is based on the use of 2-D multiresolution expansions for approximation of the unknown fields and currents. Two different 3-D open geometries, namely a rectangular dielectric resonator and a microstrip patch antenna, were studied in detail to examine the accuracy and efficiency of our novel formulation. In both cases, highly sparse moment matrices were generated, which easily render themselves to efficient iterative techniques such as the bi-conjugate gradient method. It was found out that using fast wavelet algorithm in the numerical implementation of the method of moments dramatically improves

158 the computational efficiency of the overall procedure. 6.2 Recommendations for Future Work Both of the two approaches proposed in this dissertation for the improvement of the computational efficiency of the method of moments are quite recent and have been developed by the author during the past four years. Although the scope of the integral transform technique is practically limited, it is highly effective within its range of applicability. This technique can easily be extended to three-dimensional dielectric structures, for which it is expected to yield a very high degree of computational efficiency. It can also be combined with regular planar integral equations for metallic structures to treat hybrid configurations involving both dielectric and metallic sections. Another possibility is the application of the integral transform technique to three-dimensional planar geometries with circular cross sections. In this case, the appropriate integral transform is of course a generalized Hankel transform. As for the wavelet-based MoM formulations, the present work is one of the pioneering research efforts in this field. In particular, the study of the three-dimensional planar structures are being reported for the first time. The entire field of wavelet theory itself is very young, and in the near future many new developments will certainly be witnessed. The present dissertation lays a firm foundation for the application of multiresolution analysis theory to electromagnetic problems. Although we treated only planar structures, the scope of wavelet-based formulations is obviously much wider than this type of microwave circuits as evidenced by other recent works reported in the literature. In fact, the methodology developed in this thesis can directly be applied to many other disciplines of engineering and applied sciences, where the

159 numerical modeling of various phenomena lies upon the integral equation technique. With regard to electromagnetic problems, the application of 3-D multiresolution expansions for the study of non-planar structures may open further new possibilities for the application of integral-based techniques to highly complex large-scale problems. The development of a general-purpose wavelet-based numerical code for the analysis of electromagnetic problems can be a very important milestone in the field of computational electromagnetics. From a numerical point of view, the implementation of the fast wavelet algorithm renders itself to code parallelization. A parallelized version of the present formulation will definitely surpass new limits of computational speed and efficiency. Finally, various interesting mathematical properties of wavelet expansions give an incentive to investigate the potential application of multiresolution analysis theory to other types of numerical techniques, especially differential-based methods.

APPENDICES 160

161 APPENDIX A DYADIC GREEN'S FUNCTIONS FOR PLANAR LAYERED GEOMETRIES In this appendix, the dyadic Green's functions of some planar layered background geometries are discussed. Although there are general recursive algorithms to derive the dyadic Green's functions of multilayered structures [41, 52], in the following we present explicit expressions for the Green's functions of some simple background structures, which have been considered throughout the present work. The electricfield Green's functions can be obtained by introducing different types of potentials [11]. We adopt the Sommerfeld representation based on the magnetic vector potential A. The dyadic Green's function corresponding to this potential is denoted simply by G(r I r'). Then, the electric field dyadic Green's function of the structure is given by Ge(r I') ( - I +C ) VV.G(r I r'). (A.1) Note that G(r I r') has a singularity at r = r', where the observation and source points coincide. Then, the second derivatives involved in equation (A.1) become undefined at the sined at the singular point. However, this singularity is removable and contributes a residual term after integration of the Green's function. It has been shown that

162 this integration can be performed in the Cauchy sense and the contribution of the singularity can be embodied in the form of a source dyadic L as follows [76, 77]: G,( I ) = P.V. (I + VV). G(r I r) - S(r - r), (A.2) where P.V. stands for the Cauchy principal value. The source dyadic L depends on the infinitesimal volume around the singular point, which is excluded in the course of Cauchy integration. However, the final result of the integration, i.e. the sum of the principal value and the residual term, is of course independent of the shape of the excluded volume. For the Sommerfeld representations used in this work, the source dyadic is given by [77] L = Zi. (A.3) It is well known that in an unbounded medium, an infinitesimal dipole produces a vector magnetic potential which is parallel to itself. It has been shown by Sommerfeld that in the presence of an infinite substrate, an infinitesimal vertical dipole produces a vertical magnetic potential, while an infinitesimal horizontal dipole produces a vector magnetic potential which has both horizontal and vertical components [78]. For an arbitrarily oriented dipole, the potential can be expressed as the superposition of the two horizontal and vertical cases. It is advantageous to decompose the total vector magnetic potential into two components: a primary component which is equal to the potential of the same dipole in an unbounded medium, and a secondary component which accounts for the presence of the substrate. This is written in the following manner: G(r I r') = Gp(r I r') + Gs(r I r'). (A.4) In this form, the primary component contains the source singularity, while the secondary component contains substrate poles.

163 The dyadic Green's function of layered structures are usually expressed as a spectral decomposition into elementary plane waves as follows: G(r r) (2 | ) G (,77, I z') e-(x' )e-37(y-')d dq, (A.5) (27) J-OO 2 -o00 where the substrate layers have been assumed to be parallel to the x-y plane and normal to the z-axis. Equation (A.5) is indeed a two-dimensional inverse Fourier transform, with G(, r, z I z') being the Fourier transform of the Green's function. It is seen that the x- and y-dependences in this equation are of convolutional type. The Fourier transform of the electric field Green's function is then given by 1(L 2 1L )(2 z ( ) (A.6) Ge(,,z I z') =(+ -- 7). G(7,v z I') - ~ 8z —z (A.6) where V denotes the Fourier transform of the gradient operator. Before proceeding to examine specific background geometries, note that for a twodimensional geometry of infinite extent along the y-axis with a propagating mode of propagation constant /3 along this direction, the y-dependence of the fields will be in the form of e-jPY. In this case, equation (A.5) can be integrated analytically with respect to variable y, and using the identity: L ej('1-)Ydy = 27r (r7-), (A.7) — 00 one can obtain the two-dimensional dyadic Green's function in the following manner: 1 Joo 2 G(r I r') - j G(, 3, z I z') eC-J(x-')d~. (A.8) A.1 The Unbounded Space For an unbounded medium filled with a dielectric material of relative permittivity rc), the secondary component is equal to zero. The dyadic Green's function for this

164 medium is given by the well-known formula: - _ - Jkcr-r'l G(r r') = I 47r I r - r" 1 0 f00 p-^C\Z-Z'\ -I __ [ [0 __ ^_(X-X')e3(YY')dd, (A.9) (27r)2 -oo -oo 2 e- -')- )d (A.9) where kc = V/k0o, and + -jJv/rko-2 _ _, 2 + 2 < ercko2 (A.10) It is useful to introduce the following spectral matrices: 1 -'2/(eko2) -r(/(Erko2) ~jC/(Erko2) D (E,, 1r) =- -,7/(Erko) 1 - 2/(rko2) ~j774/l(ko) | ~jfC/(f ko2) ~+j77C/(erk) (2 + 72)/(Erko2) (A.11) Then, the Fourier transform of the electric field Green's function of the unbounded medium can be written in the following form: =f e-CcIZ-Z'l P.V. Ge(~, 7, I z') = D (X, rc) 2-, (A.12) 2G, where the + and - signs correspond to the cases z > z' and z < z', respectively. A.2 Unbounded Half-space For an unbounded half-space with an infinite perfectly conducting ground plane located at z = 0, the dyadic Green's function can easily be found by employing

165 the image theory. In this case, the dyadic Green's function for the vector magnetic potential in the upper half-space is given by G= r l _"I e*-jkcR.so - kcRiO G(r r') = I — R- + I -J47o (A.13) 47rRso 47Rio where Rso = J - x'S2 + (y - y)2 + (Z - z')2 is the distance from the source point, Ro = I(x 2 + (y - y')2 + (z + z')2 is the distance from the image point, and the dyadic Ii is defined in the following way: -1 0 0 Ii= 0 -1 0 0 0 1 (A.14) The Fourier transform of equation (A.13) is given by G( z j z') = - [I eIz-z'l + I e-C(z+z')] (A.15) A.3 Single-layer Conductor-backed Substrate Consider a single-layer dielectric substrate of thickness d, with a perfectly conducting ground plane located at z = -d. The relative permittivities of the substrate and the unbounded cover medium are denoted by c,, and fC, respectively. The dyadic Green's function for the vector magnetic potential has the following general form: G~ 0 0 G(rl r') O Gyy 0 Gz Gzy Gzz (A.16)

166 In the cover medium (region I), i.e. z > 0, the Fourier transform of components of the dyadic Green's function are given by G = 2 [e-(z-' + RTE e- (z+z')] =,y, VV 1 [_c rlr, e I'>,yI,,G GZy - j(Ers - Crc) e(c(z+Z) r D? DTMDTE (A.17) where the source point is assumed to be in the upper half-space, and the various quantities in the above equations are defined as follows: DTE = Cc+ CcothC d, DTM = ErCc + Ct tanh Cd, NTE = Cc - coth Cd, NTM = r(Cc -s C tanh sd, RTE = NTE/DTE, RTM = NTM/DTM. (A.18) Note that the Green's function has first-order poles which are the roots of the following characteristic equations: DTE(, r/) = 0, DTM(~, r) = 0. (A.19) These roots are identical to the propagation constants of the TE and TM surface wave modes of the substrate [80]. Since the dominant TMo mode of an infinite conductor

167 backed dielectric layer has zero cut-off, the Green's function of this structure always has at least one surface wave pole. From equations (A.17), it is seen that the Green's function in the cover medium is in the form of the superposition of primary and secondary components, and the functions RTE and RTM may be interpreted as the reflection coefficients of the interface between the cover medium and the dielectric substrate. Then, in view of definition (A.11), the Fourier transform of the electric field Green's function in this region can be written in the following form: Ge( 77, TZ z') = D (, 5, 6rc).Gp([ i, Z I z') + D (^ ( 7, Ec). Gs(, r, Z Z'). (A.20) Inside the substrate (region II), i.e. -d < z < 0, the components of the dyadic Green's function are given by G II 1 sinh(z + d):,,e- (CZ- V = x, y, DTE sinhd' -= "ll = Crs cosh ~,(z + d) -eCZ' DTM cosh d d G - Gi -j(Ers, -rc) cosh ((z + d) e_ 7r DTMDTE cosh 4,d (A.21) A.4 Double-layer Conductor-backed Substrate Fig. A.1 shows the geometry of an infinite two-layer substrate with a perfectly conducting ground plane. The thickness and relative permittivity of the two layers are denoted by d1, Ersi, d2, and ~rs2, respectively, and a cover medium of relative permittivity Arc is assumed to fill the upper half-space. Here again, we assume that

168 Erc X dl ~rs1 d2 Figure A.1: Geometry of a double-layer grounded substrate. the source point is located in the upper half-space. Then, in the cover medium (region I), i.e. z > 0, the Fourier transform of components of the dyadic Green's function are given by G 1, l [e —'l+R EC(z+z')]] = y, G111- 2([ e +RT, v,y -z 2Cc I = -1 [ e-(Iz-''+RTM'+ e-C(z+z')] Gz~ GI zyI -jrc ____ r rsl =1 F4F5 + ( _ 1 F3 e7 T7 DTMDTE FE rs2 (A.22) and in the first substrate layer (region II), i.e. -d1 < z < 0, we have the following expressions: GV = sinh ( ) 1 - [ coth s(z + di) + s2 coth s2d2 ] e-Z', v = x,y, DTE sinh ~ssi dJ ~ II - ri, cosh 8 1(z + d1) [ 1 -} ~tanh,1(: - + d1)tanh(s2d2l e-cz' ZZ DTM cosh Fldl [,s2 /'l

169 G_ Gz _ -jErsi coshsl( (+d1) [F + tanh ( (z + dl) eC, < r7 DTMDTE cosh sldl tanh (sldl (A.23) where the various quantities are defined as follows: DTE = (ci+CS21F2 DTM = ers lcF4 + ErcF3, NTE = CFI-(2iF2, NTM = erCcF4 - ErcF3, RTE = NTE/DTE, RTM = NTDTM, F1 = 4i coth (51dl + (s2 coth (2d2, F2 1 + -2 coth (ldl coth Cs2d2, ~s1 qs ah(d Ersl F3 =,l tanh lSdl +-C2 tanh (s2d2, Ers2 F4 = 1+ rsl 2 tanh ldl tanh s2d2, Ers2 (sl F5 = Ersl (1 - ) (l coth Cld + (El 1 (s2 coth (2d2, E rc frs2 Erc F6 = -F5 _-(_-)(c, Crsl Crs2 F7 = rc (s2 tanh Cldl tanh (2d2 Fs + (-rsl (+ -rc F3 Ers2 Csl Ors2 Ersl (A.24) and ~;2 = 2 + r2 _ Erik2 i C, S1, S2. (A.25)

170 APPENDIX B METHOD OF EXTRACTION OF SINGULARITIES The Green's function encountered in the present work have weak singularities in the form of surface wave poles or branch points. The evaluation of the elements of moment matrices involves numerical integration of these Green's functions. Although this type of singularities are removable in the course of integration, ignoring them may lead to numerical instability and serious convergence problems. There are several methods for the treatment of removable singularities. In this work, we have used the method of extraction of singularities, which is based on the idea of decomposing the improper integral into a well-behaved part plus another part with an explicit singularity and evaluating the latter integral analytically [79]. In the following, note that the integration is carried out on the real axis, but the singularities may be complex in general, and necessary contour deformation must be performed to capture the singularities.

171 B.1 Extraction of Branch Point Singularities This type of singularities often appear in the form of the following improper integral: LZ2 f(z) I1 = (2 f) 2dz a> 1 a EZ, Jz (2 - 2)o/2 Imzl =Im z2 = 0, z1 < Re zo < z2, (B.1) where f(z) is an analytic function free from any singularities. Then, letting F(z) f(z)(z + Zo)-a/2, one can write I [2 -F(z) F(z) (dzB.2) I, = z + ( z- zo)~/2 (B.2) Z 2 (z- zo)a/2 The first integral of equation (B.2) is well-behaved and can be evaluated numerically without any difficulty. In fact, at the singular point z = zo, the integrand of this integral vanishes provided that F(z) has an analytic derivative at this point. The second integral of (B.2) is evaluated analytically to yield the following expression: -F(zo) - (z-1o)/-2 -_( - a/o)1-/2 ].(B.3) 1 - 0/2 (Z2 - Z0) /] It is important to choose the right branches of the square roots in equation (B.3). B.2 Extraction of Surface Wave Poles As discussed in appendix A, the Sommerfeld integrals contain first-order poles, which are identical to the TE and TM eigenvalues of the infinite substrate geometry. Consider the following integral: rI = f Z N(z) 12 I f(z)dz, f(z) - (z) (B.4) Jz d D(z)

172 Im(z) 01 02 Z1 Z2 Re(z) log branch cut Figure B.1: Branch cut for the logarithmic function in the extraction of surface wave poles. where D(zo) =, Im Im m z =, Z1 < Rezo < z2. (B.5) This integral can be rewritten in the following manner: Z2 = R [()- Resf(z)lz= + Res f(z)l=zo., (B.6) Jc z\l [Z - ZO J z - Zo where N(zo) Res f(z)\l=o = lim(z - zo) ) = (z) (B.7) z-.*rzo DI(zo) The first integral of equation (B.6) is again well-behaved and is computed numerically. The second integral is evaluated analytically to yield: Res f(z)l,=z. [ln(z2 - Z) - ln(zi - Zo)]. (B.8) The logarithmic function ln(z - zo) is a multi-valued function, and the right branch must be selected. Fig. B.1 shows the branch cut for this function at z = zo. Setting Z - Z - I Z - Zo edl0, Z2-ZO = Z2-Zo | e'j2, (B.9)

173 where Im Zo 91 = tan-' + zr, I z - Re zo I tan' Im zo 02 = tan-1 - ZRo (B.10) I z2 - Re zo | equation (B.8) reduces to the following form: D ~( \\Im zo Im) } Res f(z)l=Z. { n Z2Z - -j ltan + tan + t[ z\ - Z - - | Rezo Z2 - Re t-Reo (B.11) Note that when the pole is located on the real axis, the two tan-1 terms of equation (B.11) vanish and we have the usual jir residual contribution.

174 APPENDIX C THE BATTLE-LEMARIE MULTIRESOLUTION ANALYSIS One of the most interesting properties of multiresolution expansions is that the entire basis functions are generated by simple dilation and shifting of the two characteristic functions k and b. The construction of these characteristic functions, however, is itself a complicated task in general. Reference [50] presents a general procedure for the construction of the scaling function and mother wavelet of a multiresolution analysis (MRA) based on the properties listed in the beginning of chapter 3. In signal processing applications, wavelets are usually employed in the context of discrete transforms, and the expansion coefficients are obtained through convolution with the discrete filters {hn} and {gn} as discussed in section 3.3 of chapter 3. A procedure, called the cascade algorithm, which is based on the two-scale properties (3.8), is often used to construct the scaling function and mother wavelet given the sequences {hn} and {gn} [50]. This algorithm provides an efficient way of constructing complicated wavelets such as Daubechies' compactly supported wavelets. In the present work, we have used cubic spline Battle-Lemarie wavelets. These wavelets have good regularity, perfect symmetry and sufficient decay. Most important of all, there exist analytical expressions for the Fourier transforms of the Battle-Lemarie scaling function and

175 mother wavelet. This appendix briefly describes the construction of the Battle-Lemarie family of orthonormal wavelets. For the mathematical proofs, the reader is referred to the monograph [50]. The polynomial spline functions are the building blocks for the construction of this type of multiresolution analysis. We denote by BN(z) the Nthorder B-spline function. The Fourier transform of this function is given by N (S) —sin /2 N+1 (C. BN(() = ) I (C.1) Note that the zeroth- and first-order B-splines are in fact the usual pulse and triangular basis functions, respectively, which are traditionally used in the method of moments. Based on an orthonormalization process which is described in [50], the Fourier transform of the Nth-order Battle-Lemarie scaling function can be expressed in terms of the Fourier transform of the Nth-order B-spline function in the following way: )(N = BN(7) 2 (C.2) [El_- I BN( + 2rl) 1]1/2' Note that the denominator on the right hand side of equation (C.2) is a 27r-periodic function of ~, and can be expressed in the form of an Nth-order polynomial in sin2(~/2). Now define MN(() = (N)() (C.3) This function is also 27r-periodic and has the following Fourier series expansion: MN() = - E hE N ) e-jn, (C.4) X nEZ where Z is the set of integers. The Fourier coefficients of (C.4) are indeed identical to the discrete sequence {hn} of equation (3.8). It can be shown, albeit through

176 a complicated proof, that the Fourier transform of the Nth-order Battle-Lemarie mother wavelet is then given by the following equation: N)(N) -= e/2 M ( +,r) +(N)((), (C.5) where MA (,) denotes the complex conjugate of MN((). Throughout this thesis, the cubic spline Battle-Lemarie multiresolution analysis (N = 3) has been used for the moment method expansions. The Fourier transforms of the cubic spline Battle-Lemarie scaling function and mother wavelet are given by the following closed-form expressions: -(3) 1/2,(C.6) P (sin2 ~)]/2 and.(e)()= e/( sin ( s2 ) 1/2, e4 sin C P (sin2 ]) P (sin2 ) (C.7) where P(x) is a cubic polynomial given by 4 22 4 3 P(x)=l —+ --. (C.8) P() = 1 3 5 315 (C+8) In view of equations (C.6) and (C.7), it is not difficult to show that the cubic spline Battle-Lemarie scaling function and mother wavelet can be expressed in terms of cubic B-spline functions in the following way: +(3)(X) = b3)B3(x- k), kEZ =+(3)(x) CE C) B3(2 - k). kEZ (C.9)

177 k b(3) c3) h?3) k k k 0 1.96976160 -2.00521039 0.76613005 1 -0.67243039 2.89173391 0.43392265 2 0.26870415 -2.00521039 -0.05020169 3 -0.11851986 0.54227884 -0.11003702 4 0.05519138 -0.01207123 0.03208086 5 -0.02652026 0.14408849 0.04206833 6 0.01299809 -0.14591233 -0.01717627 7 -0.00645742 0.00301822 -0.01798229 8 0.00323979 0.02834407 0.00868520 9 -0.00163770 0.01914911 0.00820140 10 0.00083276 -0.02246183 -0.00435366 Table C.1: Characteristic coefficients of the Battle-Lemarie multiresolution analysis. with b3) = b(3 and ck) = c2 Table C.1 gives the values of some of these expansion coefficients. Equations (C.9) imply that the Battle-Lemarie functions are essentially superpositions of simple polynomial functions (see [63]).

178 APPENDIX D METHOD OF STATIONARY PHASE The numerical integration of complex functions with very rapid oscillations is a difficult computational task and often leads to erroneous results due to the cancelation effect. This type of integrals are frequently encountered in certain problems such as the evaluation of far fields, which involve complex exponentials of very large arguments. In such cases, it is possible to derive asymptotic forms of the relevant integrals in an analytical fashion. To this end, we use the method of stationary phase which is described as follows [52, 53]. First, we consider the following one-dimensional integral: I(Q) = f(x) ej9g(x) dx, Q > 0, (D.1) 1 where f(x) is a slowly varying, well-behaved function, S is a large parameter, and g(x) is a real function with a continuous derivative in xz <~: < x2, and a stationary point at x = x,, i.e., g'(x) = 0, x1 < xs, < x2. (D.2) In the vicinity of the stationary point, the Taylor expansion of the function g(x) is given by g(x) = g(x,) + i g"(x,)(x - x8)2 +... (D.3)

179 By extracting the neighborhood of the stationary point from integral (D.1) and applying integration by parts to the two resulting integrals, it can be shown that the following asymptotic expression for equation (D.1) is obtained: f7rx3) 2wlg(x,)+ Oa I (Q). fx) I g"(xS) e ] + 2 [ f(X2) ejg(X2) _ f(Xl) ejg(x1) 9g'(2)'(x1) J + O(n-2), (D.4) where = sgn g"(Xs). (D.5) Note that the contributions from the stationary point and the two endpoints are of orders 1/v/7 and 1/fI, respectively. The above result can be extended to the two-dimensional case in a straight-forward manner. In particular consider the double integral I()= / f(x, y) eig("Y) dxdy, Q > 0, (D.6) J -00 -00 where f(x, y) is a slowly varying function, and g(x, y) is a peal function with a stationary point at (x8, y,) given by ag(xs,Ys) ag(x,,ys) = 0 (D.7) ax ay Then, the asymptotic form of integral (D.6) is given by I2(n). Rf(xy) I det(D) 11/2] e, (D.8)

180 where the matrix D is defined as / a2g(x2YS) ag( x.y,) D= ax xay, (D.9) a2g(rx,,y) a2g(x3,ys) ay\ x y2 and a = sgn dl + sgn d2, (D.10) with dl and d2 being the eigenvalues of D. In this case, the endpoint contributions vanish.

181 APPENDIX E THE BI-CONJUGATE GRADIENT METHOD The solution of linear systems of equations is one of the oldest problems of numerical analysis. There are a large number of technique to solve linear system of the form A.x = b, (E.1) where A is the coefficient matrix, and x and b are the unknown and constant vectors, respectively. In the context of the present work, these quantities correspond to the moment (impedance) matrix, and amplitude and excitation vectors, respectively. The methods of solution of equation (E.1) are traditionally divided into two categories: direct methods and iterative methods. Based on these two general categories, specialized techniques have been developed for solving very large sparse linear systems like those encountered in chapter 5 [55]. In this work, we have used a pre-conditioned Bi-Conjugate Gradient (BiCG) method, which belongs to the second category and offers a very high degree of efficiency for the type of matrices encountered in electromagnetic problems [59, 60]. In the following, the adjoint of a matrix is defined as A = A*T, (E.2)

182 where * denotes the complex conjugate, and the superscript T denotes the transpose of a matrix. Moreover, the following inner product is defined: < xl, X2 >= -xT. X2. (E.3) The iteration process starts with an initial solution xo (guess vector). An intuitive choice is xo = 0. Let A denote a pre-conditioner matrix, which must be fairly close in form to A-'. Then, we have A.A.x = A.b. (E.4) If A happens to be identical to the inverse of A, then the solution is complete. Since the moment matrices in wavelet formulations are highly sparse with the biggest elements on their diagonals, a good choice for the pre-conditioner matrix is the inverse of the diagonal part of A. We introduce the vector pairs (rk, fk), (Zk, Zk), and (pk, pk). The initial choices for these vectors are as follows: o = b - A. ro = ro A. zo to Po = Zo Po Zo (E.5) Then the iteration process proceeds in the following order: < Tk, Zk > < k, A. P >'

183 Xk+l = Xk + cXkpk, rk+l - rk - CkA.pk, rk+l - rck — kA'. Pk A. Zk+l = rk+l, Aa. Zk+l = rk+l, P < rf+l Zk+l > Pk = < Tk, Zk > Pk+1 Zk+1 + Pk Pk, pfk+i = Zkc+l + iPk. (E.6) A relative error of the following typical form is defined and tested in each cycle of iteration: I b-A. xk+l Ek= Ibi (E.7) The mathematical proof of convergence of this scheme can be found in [81].

BIBLIOGRAPHY 184

185 BIBLIOGRAPHY [1] E.A.J. Marcatili,"Dielectric rectangular waveguide and directional coupler for integrated optics," Bell Syst. Tech. J., vol.48, pp. 2071-2102, Sept. 1969. [2] J.E. Goell,"A circular-harmonic computer analysis of rectangular dielectric waveguides," Bell Syst. Tech. J., vol.48, pp. 2133-2160, Sept. 1969. [3] R. Mittra, Y.L. Hou and V. Jamnejad,"Analysis of open dielectric waveguides using mode-matching technique and variational methods," IEEE Trans. Microwave Theory Tech., vol. MTT-28, pp. 36-43, Jan. 1980. [4] E. Schweig and W.B. Bridges, "Computer analysis of dielectric waveguides: a finite difference method," IEEE Trans. Microwave Theory Tech., vol. MTT-32, pp. 531-541, May 1984. [5] J.S. Bagby, D.P. Nyquist and B.C. Drachman," Integral formulation for analysis of integrated dielectric waveguides," IEEE Trans. Microwave Theory Tech., vol. MTT-33, pp. 906-915, Oct. 1985. [6] R.B. Wu and C.H. Chen,"On the variational reaction theory for dielectric waveguides," IEEE Trans. Microwave Theory Tech., vol. MTT-33, pp. 477-483, June 1985. [7] K. Hayata, M. Koshiba, M. Eguchi and M. Suzuki, "Vectorial finite-element method without any spurious solutions for dielectric waveguiding problems using transverse magnetic-field component," IEEE Trans. Microwave Theory Tech., vol. MTT-34, pp. 1120-1124, Nov. 1986. [8] K. Bierwirth, N. Schulz and F. Arndt, "Finite-difference analysis of rectangular dielectric waveguide structures," IEEE Trans. Microwave Theory Tech., vol. MTT-34, pp. 1104-1113, Nov. 1986. [9] D. Yevick and B. Hermansson,"New formulations of the matrix beam propagation method: application to rib waveguides," IEEE J. Quantum Electron., vol.25, pp. 221-229, Feb. 1989. [10] K. Wu and R. Vahldieck, "Comprehensive MoL analysis of a class of semiconductor-based transmission lines suitable for microwave and optoelectronic application," Int. J. Num. Modelling: Electronic Networks, Devices and Fields, Vol.4, pp. 45-62, 1991.

186 [11] T. Itoh, Ed., Numerical Techniques for Microwave and Millimeter- Wave Passive Structures. New York, NY: John Wiley & Sons, 1989. [12] S.T.Chu and S.K. Chaudhuri, "Combining modal analysis and the finitedifference time-domain method in the study of dielectric waveguide problems," IEEE Trans. Microwave Theory Tech., vol. MTT-38, pp. 1755-1760, Nov. 1990. [13] A.G. Engel, N.I. Dib and Linda P.B. Katehi," Characterization of a shielded transition to a dielectric waveguide," IEEE Trans. Microwave Theory Tech., vol. MTT-42, pp. 847-854, Sept. 1994. [14] J.A. Pereda, L.A. Vielva and A. Perieto, "Computation of resonant frequencies and Quality factors of open dielectric resonators by a combination of the finite-Difference time-domain (FDTD) and Prony's methods," IEEE Microwave Guided wave Lett., vol.2, pp. 431-433, Nov. 1992. [15] A. Bayliss and E. Turkel,"Radiation boundary conditions for wave-like equations," Comm. Pure Appl. Math., vol. 33, pp. 707-725, 1980. [16] B. Engquist and A. Majda,"Absorbing boundary conditions for wave-like equations," Math. Comp., vol. 31, pp. 629-651, 1977. [17] J. Berenger,"A perfectly matched layer for the absorption of electromagnetic waves," J. Comp. Phys, vol. 114, No. 2, pp. 185,200, Oct. 1994. [18] K.K. Mei,"On the integral equations of thin-wire antennas," IEEE Trans. Antennas Propagat., vol. AP-13, pp. 374-378, May 1965. [19] E.W. Kolk, N.H.G. Baken and H. Block," Domain integral equation analysis of integrated optical channel and ridge waveguides in strtified media," IEEE Trans. Microwave Theory Tech., vol. MTT-38, pp. 78-84, Jan. 1990. [20] J.F. Kiang, S.M. Ali and J.A. Kong," Integral equation solution to the guidance and leakage properties of coupled dielectric strip waveguides," IEEE Trans. Microwave Theory Tech., vol. MTT-38, pp. 193-203, Feb. 1990. [21] H.Y. Yang, J.A. Castaneda and N.G. Alexopoulos," An integral equation analysis of an infinite array of rectangular dielectric waveguides," IEEE Trans. Microwave Theory Tech., vol. MTT-38, pp. 873-880, July 1990. [22] K. Sabetfakhri and P.B. Katehi,"An integral transform technique for analysis of planar dielectric structures," IEEE Trans. Microwave Theory Tech., vol. MTT42, pp. 1052-1062, June 1994. [23] K. Sabetfakhri and P.B. Katehi,"Analysis of integrated millimeter-wave and submillimeter-wave waveguides using orthonormal wavelets expansions," IEEE Trans. Microwave Theory Tech., vol. MTT-42, pp. 2412-2422, Dec. 1994.

187 [24] R.F. Harrington, Field Computation by Moment Methods. New York: Macmillan, 1968. [25] G. Beylkin, R. Coifman and V. Rokhlin,"Fast wavelet transforms and numerical algorithms I," Commun. Pure Appl. Math., vol. 44, pp 141-183, 1991. [26] I. Daubechies,"Orthonormal bases of compactly supported wavelets," Comm. Pure Appl. Math., vol. XLI, pp. 909-996, 1988. [27] S.G. Mallat,"A theory for multiresolution signal decomposition: The wavelet representation," IEEE Trans. Pattern Anal. Machine Intell., vol. PAMI-7, pp. 674-693, July 1989. [28] S.G. Mallat,"Multifrequency channel decompositions of images and wavelet models," IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-37, pp. 20912110, Dec. 1989. [29] S.G. Mallat,"Multiresolution approximation and wavelet orthonormal bases of L2," Trans. Amer. Math. Soc., vol. 3-15, pp. 69-87, Sept. 1989. [30] A.E. Yagle,"Image reconstruction from projections under wavelet constraints," IEEE Trans. Signal Proc., vol.41, pp. 3579-3584, Dec. 1993. [31] A. Grossman and J. Morlet, "Decomposition of Hardy functions into square integrable wavelets of constant shape," SIAM J. math. anal., vol. 15, no. 4, pp. 723-736, July 1984. [32] B.Z. Steinberg and Y. Leviatan,"On the use of wavelet expansions in the method of moments,"IEEE Trans. Antennas Propagat., vol. 41, pp. 610-619, May 1993. [33] R.L. Wagner, G.P. Otto and W.C. Chew,"Fast waveguide mode computation using using wavelet-like basis functions," IEEE Microwave Guided Wave Lett., vol. 3, pp. 208-210, July 1993. [34] G. Wang and G.-W. Pan, "Full wave analysis of microstrip floating line structures by wavelet expansion method," IEEE Trans. Microwave Theory Tech., vol. MTT-43, pp. 131-142, Jan. 1995. [35] M. Swaminathan, T.K. Sarkar and A.T. Adams,"Computation of TE and TM modes in waveguides based on a surface integral equation," IEEE Trans. Microwave Theory Tech., vol. MTT-40, pp. 285-297, Feb. 1992. [36] T.E. van Deventer and P.B. Katehi,"A planar integral equation method for the analysis of dielectric ridge structures using generalized boundary conditions," in IEEE MTT-S Dig., 1991, pp. 647- 650. [37] T.E. van Deventer,"Characterization of two-dimensional high frequency microstrip and dielectric interconnects," Ph.D. dissertation, University of Michigan, Dec. 1992.

188 [38] K. Sabetfakhri and P.B. Katehi," A study of open dielectric waveguide problems using the generalized integral equation method," in URSI Radio Sciences Meeting, Chicago, July 1992. [39] N. Bleistein and R.A. Handelman, Asymptotic Expansions of Integrals. New York, NY: Holt, Rinehart and Winston, 1975. [40] C.T. Tai,"Some essential formulas in dyadic analysis and their applications," Radio Science, vol. 22, pp. 1283-1288, Dec. 1987. [41] S. Barkeshli and P.H. Pathak,"On the dyadic Green's function for a planar multilayered dielectric/magnetic media," IEEE Trans. Microwave Theory Tech., vol. MTT-40, pp. 128-142, Jan. 1992. [42] S.T. Peng and A.A. Oliner,"Guidance and leakage properties of a class of open dielectric waveguides: Part I-Mathematical Formulations," IEEE Trans. Microwave Theory Tech., vol. MTT-29, pp. 843-855, Sept. 1981. [43] N.K. Das and D.M. pozar,"Full-wave spectral-domain computation of material, radiation and guided wave losses in infinite multilayared printed transmission lines," IEEE Trans. Microwave Theory Tech., vol. MTT-39, pp. 54-63, Jan. 1991. [44] A.A. Oliner,"Leakage from higher modes on microstrip line with application to antennas," Radio Science, vol. 22, pp. 907-912, Nov. 1987. [45] G. Sansone, A.H. Diamond and E. Hille, Orthogonal Functions. New York, NY: Interscience Publishers,Inc., 1958. [46] T. Tamir, Ed., Guided-Wave Optoelectronics. Berlin: Springer-Verlag, 1988. [47] W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery, Numerical Recipes in FORTRAN: The Art of Scientific Computing. Cambridge: Cambridge University Press, 1992. [48] K. Sabetfakhri and P.B. Katehi," On the application of quasi-wavelet expansions to open dielectric waveguide Problems," in Proc. Antennas Propagat. Soc. Int. Symp., Ann Arbor, July 1993. [49] P.J. Davis and P. Rabinowitz, Numerical Integration. Waltham, MA: Blaisdell Publishing Co., 1967. [50] I. Daubechies, Ten Lectures on Wavelets, Philadelphia, PA: Society for Industrial and Applied Mathematics, 1992. [51] C.K. Chui, ed., Wavelets: A Tutorial in Theory and Applications. New York: Academic Press, 1992. [52] J.A.Kong, Theory of Electromagnetic Waves. New York, NY: John Wiley & Sons, 1975.

189 [53] L.B. Felsen and N. Marcuwitz, Radiation and Scattering of Waves. Englewood Cliffs, NJ: Prentice-Hall, 1973. [54] A. Papoulis, The Fourier Integral and its Applications. New York, NY: McGrawHill, 1962. [55] D.J. Evans, Ed., Sparsity and its Applications. Cambridge, UK: Cambridge University Press, 1985. [56] N.I. Dib, Personal communication. [57] b. Young and T. Itoh, "Analysis and design of microslab waveguide," IEEE Trans. Microwave Theory Tech., vol. MTT-35, pp. 850-857, Sept. 1987. [58] I.S. Duff,"A survey of sparse matrix research," Proc. IEEE, vol. 65, No. 4, Apr. 1977. [59] T.K. Sarkar,"The application of the conjugate gradient method for the solution of operator equations arising in electromagnetic scattering from wire antennas," Radio Science, vol. 19, pp. 1156-1172, Sept.-Oct. 1984. [60] C.F. Smith, A.F. Peterson and R. Mittra,"The biconjugate gradient method for electromagnetic scattering," IEEE Trans. Antennas Propagat., vol. AP-38, pp. 938-940, June 1990. [61] C. Lanczos,"An iteration method for the solution of the eigenvalue problem of linear differential and integral operators," J. Research Nat'l. Bureau Standards, vol. 45, No. 4, pp. 255-282, Oct. 1950. [62] A. Papoulis, Signal Analysis. New York, NY: McGraw-Hill, 1977. [63] K. Sabetfakhri, "On the numerical aspects of Battle-Lemarie multiresolution analysis," Michigan Radiation Lab. Rep., April 1995. [64] M. Abramowitz and I.A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables. New York, NY: Dover, 1972. [65] R.K. Mongia, "Theoretical and experimental resonant frequencies of rectangular dielectric resonators," IEE Proc.-H, Vol.139, pp. 98-104, Feb. 1992. [66] K.R. Carver and J.W. Mink, "Microstrip antenna technology," IEEE Trans. Antennas Propagat., vol. AP-29, pp. 2-24, Jan. 1981. [67] D.M. Pozar and S.M. Voda,"A rigorous analysis of microstripline fed patch antenna," IEEE Trans. Antennas Propagat., vol. AP-35, pp. 1343-1350, Dec. 1987. [68] D.H. Schaubert,"Microstrip antennas," Electromagnetics, vo. 12, pp. 381-401, 1992.

190 [69] M.C. Bailey and M.D. Deshpande,"Integral equation formulation of microstrip antennas," IEEE Trans. Antennas Propagat., vol. AP-30, pp. 651-656, July 1982. [70] D.M. Pozar,"Input impedance and mutual coupling of rectangular microstrip antennas," IEEE Trans. Antennas Propagat., vol. AP-30, pp. 1191-1196, Nov. 1982. [71] M.C. Bailey and M.D. Deshpande," Input impedance of microstrip antennas," IEEE Trans. Antennas Propagat., vol. AP-30, pp. 645-650, July 1982. [72] R.C. Hall and J.R. Mosig,"The analysis of coaxially fed microstrip antennas with electrically thick substrates," Electromagnetics, vo. 9, pp. 367-384, 1989. [73] J.T. Aberle and D.M. Pozar,"Accurate and versatile solutions for probe-fed microstrip patch antennas and arrays," Electromagnetics, vo. 11, pp. 1-19, 1991. [74] R.F. Harrington, Time-Harmonic Electromagnetic Fields. New York: McGrawHill, 1961. [75] J.R. Mosig and F.E. Gardiol,"General integral equation formulation for microstrip antennas and scatterers," Proc. IEE, Pt. H, vol. 132, pp. 424-432, 1985. [76] A.D. Yaghjian,"Electric dyadic Green's functions in the source region," Proc. IEEE, vol. 68, pp. 248-263, Feb. 1980. [77] M.S. Viola and D.P. Nyquist," An observation on the Sommerfeld-integral representation of the electric field dyadic Green's function for layered media," IEEE Trans. Microwave Theory Tech., vol. MTT-36, pp. 1289-1292, Aug. 1988. [78] A. Sommerfeld, Partial Differential Equations in Physics. New York, NY: Academic Press, 1964. [79] P.B. Katehi and N.G. Alexopoulos,"Real axis integration of Sommerfeld integrals with applications to printed circuit antennas," J. Math Phys., vol. 14, pp. 527533, 1983. [80] R.E. Collin, Field Theory of Guided Waves. New York, NY: IEEE Press, 1991. [81] L.A. Hageman and D.M. Young, Applied Iterative Methods. New York, NY: Academic Press, 1981.

UNIVERSITY OF MICHIGAN 3 9015 03526 8823