RL-944 ELECTROMAGNETIC MODELING OF HIGH-SPEED HIGH-FREQUENCY INTERCONNECTS by Jong-Gwan Yook A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy (Electrical Engineering) in The University of Michigan 1996 RL-944 = RL-944

Jong-Gwan Yook 1996 All Rights Reserved

This thesis is dedicated to my parents in memory of their devotion and selfless love ii

ACKNOWLEDGEMENTS This piece of work presented as a form of Ph.D. dissertation would not have been possible without help of many individuals, of whom it is almost impossible to mention all, encountered in my life. Herein, I would like to take a chance to thank those people who have encouraged and cultivated my intellectual growth and enriched my daily life in Ann Arbor for the last five and half years, and hopefully forever. My first sincere gratitude goes to my advisor, Professor Linda Katehi, since without her support, understanding, and technical expertise, this research would have been a much more difficult and painful journey. Her enthusiasm for research and teaching has made her a symbol for me of a devoted scientist and educator. I would also like to thank all of my committee members: Professor Karem Sakallah for his guidance on the digital circuit aspect of my research, Professor Kamal Sarabandi for his sincere interest in and suggestions on my dissertation, Dr. Rainee Simons in NASA Lewis Research Center for his expertise on high-frequency packages and transitions, and Dr. Hal Marshall for his endless support on parallel computing. The diversity and abundance of technical facilities and programs at the Radiation Laboratory in the University of Michigan have led me to taste and realize the importance of interdisciplinary educations. In particular, my friends, Dr. Paul Siqueira, Dr. Kazem Sabetfakhri, Messrs. Tayfun Ozdemir, Stephane Legault, Sunil Bindiganavale, and many more, too numerous to mention, deserve my deepest thanks for iii 111

their close friendship and technical discussions immersed in this thesis. Drs. Yisok Oh and Howard Yoon have to be recognized for their encouragement and support throughout my graduate career. As a successful graduate student studying far from home, I felt always comfortable and enjoyable with many friends in the church, Ann Arbor Korean Catholic Student Association, and they deserve my warm appreciation. Especially, I wish to thank loving God for letting me know such good friends: Mr. Sang Yeol Park, Dr. Young Ro Byun, Messrs. Jong Yul Choi, Dong Wook Lee, Eun Suk Choi, and their families. I have received financial support from Intel Co. for the last four and half years, and I would like to thank them for their interest in the research I have done. I always found it helpful and informative to talk to people at Intel, especially Drs. Tawfik Arabi and Tim Schreyer. They have my sincere appreciation. I also wish to express my appreciation to Dr. Rod Johnson for his time and effort to review and correct my dissertation from the acknowledgement to abstract. Most importantly, I wish to express my deep gratitude to my family members: my father for his intelligence, endurance and patience, my mother for her wisdom and selfless devotion, my elder brother for his humanity, and my younger sister for her understanding and support. iv

ABSTRACT ELECTROMAGNETIC MODELING OF HIGH-SPEED HIGH-FREQUENCY INTERCONNECTS by Jong-Gwan Yook Chair: Linda P. B. Katehi This dissertation develops in detail an approach to a three-dimensional full-wave elec — trornagnetic field simulator, namely the finite element method (FEM), and applies it to various high-speed high-frequency interconnects, ranging from a simple via hole to an entire K/Ka-band MMIC phase shifter package. This approach has several essential features: the utilization of tetrahedron-based edge basis functions rendering spurious-free solutions, a non-uniform structured mesh generator giving flexible modeling capabilities, and implementation of artificial absorbing layers helping to simulate open boundary problems. In addition, the FEM is fully parallelized on a distributed memory machine (the IBM SP2). Two different parallelization strategies are implemented, among which the task parallelization strategy provides linearly scalable performance improvement due to the minimal communication overhead among the processors.

1 \V'ith all the above features, the parallelized FEM has been successfully applied to the characterization of planar/non-planar high frequency interconnects, hermetic transitions, and K/Ka-band MMIC packages. In particular, the undesirable internal package resonances and energy leakages in the IMMIC package are identified, and furthermore, a few mechanisms for the suppression of the above phenomena are suggested, and these are supported by rigorous numerical data. The lumped equivalent circuits for the vertical via holes are utilized in the system level EM modeling to provide appropriate inductances and capacitances. In this dissertation, characterization of various high-frequency interconnects has been stressed as well as parametric study. Special attention is devoted to development of the system level electromagnetic modeling for high-speed digital circuits and packages. This system level modeling is achieved by combining the FEM and the well known circuit simulator, HSPICE. It is based on the so-called tiling processor, which generates equivalent circuits for PCBs and ICs, and then the circuit simulator is employed for time or frequency domain characterizations. The hybrid approach implemented in this study renders effective simulations in frequency as well as time domain for a given geometry. The detailed derivation and applications of the system level electromagnetic modeling tool are described. Its validity and accuracy are proved by modeling the INTEL P6 board and package containing eight layers of signal, power, and ground planes and hundreds of signal traces.

TABLE OF CONTENTS DEDICATION...................................................... ii ACKNOWLEDGEMENTS.............................. iii LIST OF FIGURES..................................ix LIST OF TABLES................................... xix LIST OF APPENDICES................................. xx CHAPTER I. INTRODUCTION.............................. 1 1.1 Motivation and Background......................... 1 1.2 Objectives................................. 4 1.3 O verview.................................. 6 II. THEORETICAL METHODOLOGY: 3D VECTOR FINITE ELEMENT METHOD.............................. 9 2.1 Introduction............................... 9 2.2 Vector Wave Equation......................... 10 2.3 Weak Form of Wave Equation and Its Discretization.......... 12 2.4 Basis and Weighting Functions..................... 14 2.5 Assembly of Matrix Equation and Its Solution Methods........ 18 2.6 Treatment of Boundary and Symmetry Conditions.......... 21 2.6.1 Boundary Conditions..................... 21 2.6.2 Symmetry Conditions.................... 24 2.7 The Question of Spurious Solutions: A Rationale for Using the Edge Basis Function..................................... 28 2.7.1 The Origin of Spurious Solutions.............. 28 2.7.2 Why the Whitney Edge Element?.............. 29 2.7.3 Eigenvalue Analysis with Edge Element............. 31 2.8 Excitation Mechanism............................... 33 2.9 Customized Generation of Tetrahedral Mesh............. 35 2.10 Treatment of Open Boundary Problems....................... 39 2.10.1 Isotropic Absorbers for MMIC Applications........ 43 2.10.2 Isotropic Absorbers for Waveguide Applications....... 45 v

2.11 Passive Lumped Elements............................. 46 2.11.1 Volume Current Method....................... 48 2.11.2 Surface Impedance Method....................... 49 III. VALIDATION AND PARALLELIZATION.................. 52 3.1 Introduction............................... 52 3.2 Eigenvalue Problems.......................... 53 3.2.1 Lossless Cavities....................... 55 3.2.2 Lossy Cavities......................... 60 3.2.3 Sum m ary........................... 62 3.3 Parallelization.............................. 62 3.3.1 M otivations.......................... 62 3.3.2 Parallelization Strategies................... 63 3.3.3 Discussion........................... 71 3.4 Performance of the Artificial Absorbers................... 72 3.4.1 Microstrip Thru Line..................... 73 3.4.2 Microstrip Junction Terminated with Absorber.......... 74 3.4.3 CPW Series Stub....................... 75 3.4.4 CPW Series Stub: Open Structure................ 77 3.4.5 Waveguide Absorber.......................... 79 3.4.6 Summary................................. 80 3.5 Conclusions................................... 83 IV. PLANAR/NON-PLANAR HIGH FREQUENCY INTERCONNECTS.................................................... 84 4.1 Introduction............................... 84 4.2 Microstrip Short Circuit........................ 87 4.3 Microstrip Ground Pad......................... 93 4.4 CPW-to-Microstrip Through-Via Transition................ 98 4.4.1 W ithout Via Hole....................... 98 4.4.2 W ith Via Hole........................ 101 4.5 Channelized CPW-to-Microstrip Transition................ 103 4.6 Asymmetric 2-port Via with Plated-Through-Hole.......... 107 4.7 Double Via Hole.............................. 111 4.8 CPW-to-Slotline Transitions...................... 118 4.8.1 Planar Transition............................ 118 4.8.2 Non-Planar Transition........................ 121 4.9 CPW-to-Waveguide Transition.............................. 123 4.10 Numerical Convergence................................... 128 4.11 Conclusions............................................ 130 V. HERMETIC TRANSITIONS................................. 132 5.1 Introduction............................... 133 5.2 Intra-chip Transition: Microstrip-through-CPW Hermetic Transition 135 vi

5.2.1 Open Transition Through a CPW Section......... 5.2.2 Closed Transition Through a CPW Section........ 5.2.3 Vertical Wire Grounding................... 5.2.4 Rectangular Via Hole Grounding.............. 5.2.5 Extended Rectangular Via Hole Grounding........ 5.2.6 PEC Wall Grounding..................... 5.2.7 Microstrip-through-CPW Transition With Hermetic Wall 5.2.8 W-Band Application of the Hermetic Transition...... 5.3 Inter-Chip Transition: Hermetic Bead Transition........... 5.3.1 Squire Dielectric Ring.................... 5.3.2 Circular Dielectric Ring................... 5.4 Conclusions............................... 137 13S 140 143 146 148 149 152 154 154 156 156 VI. K/Ka-BAND MMIC PACKAGE...................... 159 6.1 Introduction............................... 6.2 Modeling of the Package........................ 6.3 Numerical and Experimental Results and Discussions........ 6.3.1 Modeling of Isolated Symmetric Package.......... 6.3.2 Modeling of Isolated Asymmetric Package......... 6.3.3 Modeling of Asymmetric Package Placed in a Test Fixture 6.4 Features Which Can Enhance Package Performance......... 6.4.1 Asymmetric Package with Additional Metal Filled Vias 6.4.2 Asymmetric Package with Resistive Coated Lid...... 6.4.3 Asymmetric Package with Lossy Dielectric Frame..... 6.5 Conclusions............................... 159 162 164 164 166 169 170 170 175 175 177 VII. SYSTEM LEVEL EM MODELING OF DIGITAL ICs AND PACKAGES........................................ 179 7.1 Introduction........................... 7.2 Problem Statement and Previous Work............ 7.3 PCB Tiles and their Lumped Equivalent Circuits...... 7.3.1 Tiling Procedure................... 7.3.2 Equivalent Circuits.................. 7.4 Simulation Results........................ 7.4.1 Convergence Test: Tiling Optimization....... 7.4.2 Seven Drivers with Decoupling Capacitors..... 7.4.3 Modeling of Signal Traces.............. 7.4.4 Split PCBs...................... 7.4.5 Intel P6 Gadget.................... 7.4.6 Summary....................... 7.5 Conclusions........................... VIII. CONCLUSIONS........................... 8.1 Summary of Achievements..................... 179.... 182.... 187.... 187... 190.... 195... 196.... 200.... 203... 205... 209.... 211.... 213 215 215 vii

8.2 Suggested Future Work....................... 217 A PPEN D ICES.......................................219 BIBLIOGRAPHY....................................... 233 o** Vlll

LIST OF FIGURES Figure 2.1 Geometry of three dimensional volume which has current sources Ji and/or Mi with possible material inhomogeneities and geometrical discontinuities. The boundary of the volume V is closed by some combination of PEC, PMC, and ABC......................................... 11 2.2 Tetrahedral element and associated nodes and edges definitions. The circled numbers 1 to 4 represent four node labels and El to E6 six edges associated with those nodes. ri1 and ri2 are node vectors associated with the ith edge and r is an arbitrary vector which lies in the tetrahedron...15 2.3 Behavior of edge element Wp for edge p {i, j}............... 17 2.4 PEC or PMC surface S with outward unit normal vector i......... 22 2.5 Material interface S with unit normal vector i................. 23 2.6 An example of a circuit configuration (microstrip line) bearing a symmetrical plane at the center and its electric and magnetic field distributions. As presented on the right, the original problem can be replaced with half of the structure, with PMC symmetry condition on one side............ 24 2.7 An example of a circuit configuration (double strip line) containing symmetry at the center and its electric and magnetic field distributions. The original problem can be replaced with the half of the structure with PEC symmetry condition at the center as shown on the right............ 25 2.8 Reduced equivalent circuits for symmetric 2-port network depending on excitation configurations: (a) original network, (b) reduced equivalent network for even mode excitation, and (c) reduced equivalent network for odd mode excitation.......................................... 26 2.9 Two types of ideal excitation configurations: (a) horizontal current source on the transmission line plane. (b) vertical current source............ 35 2.10 Two types of rectangular bricks which containing five tetrahedral elements: (a) Type I, (b) Type II....................................... 36 ix

2.11 Treatment of open boundary problem: (a) Original open boundary problem, (b) Equivalent closed/finite boundary problem............... 41 2.12 Impedance load......................................... 48 2.13 Resistive sheet and equivalent thin dielectric layer................. 50 3.1 Rectangular cavity having PEC walls and dimensions a, b, and c in x-, y-, and z-direction, respectively................................ 56 3.2 Partially-filled rectangular cavity having dielectric material with er and thickness c'................................... 57 3.3 Rectangular cavity having dielectric cube at the center of the cavity bottom having dielectric constant Er and a', b', and c' in x-, y-, and z-dimension, respectively.................................. 59 3.4 Variation of the resonant frequency of the dominant mode for various size and er of the dielectric material................................. 60 3.5 Typical resonance curve of a lossy cavity. The quality factor Q is defined based on the center frequency, fo, and two adjacent frequency points, fi and f2............................................... 61 3.6 Parallelization of the linear equation solver, BiCG method, in the FEM scheme..................................................... 64 3.7 Parallelization of the matrix-vector multiplication. Step 1: broadcast the partitioned FEM matrix [Ai]. Step 2: broadcast the vector [x] in every iteration. Step 3: perform partial matrix-vector multiplication [AJ][x] on each machine, where [A] = U [Ai] for i = 1,..., N. N is the number of slave nodes. Step 4: gather the results to form [b].................... 66 3.8 Speedup curves of the parallel matrix-vector multiplication for two different problem sizes. M indicates the number of unknowns............... 68 3.9 The concept of the task parallelization................... 70 3.10 Speedup curve for the task parallelization...................... 71 3.11 Computed return loss for the microstrip line terminated with artificial absorbers. Geometrical and material factors: L = 6.38 mm, W = 1.3 mm, D = 5.0 mm, H1 = 0.635 mm, H2 = 3.365 mm, cr = 3.4, ur = 1.0, el = 1.0 + j3.0, It, = 1.0 + j3.0, ~2 = 3.4 + j10.0, and [2 = 1.0 + j2.9412 74 x

3.12 Computed return loss for the microstrip junction terminated with artificial absorbers. Absorbers are placed at the side and top walls also in addition to the termination end. Geometrical factors: W1 = 1.3mm, W2 = 0.6mm, and T = 3.0mm. The thickness and the dielectric constant of the substrate are H = 0.635mm and Er = 3.4. The absorbers at the end are chosen as: E1 = 1.0 + j15.0 and p1 = 1.0 + j15.0 for the air side and 62 = 3.4 + j15.0 and P/2 = 1.0 + j4.412 for the dielectric side...................... 76 3.13 Computed return loss for the CPW series stub terminated with artificial absorbers (closed structure) at 4 different positions. Geometrical factors: G = 0.225mm, S = 0.45mm, L = 1.35mm, W = 0.225mm, and T = 2.0mm. The thickness and the dielectric constant of the substrate are H = 0.635mm and ~r = 9.9. The absorbers at the end are chosen as: 1e = 1.0 -+ j15.0 and p1 = 1.0 + j15.0 for the air side and 62 = 9.9 + j20.0 and Pu2 = 1.0 + j2.02 for the dielectric side....................... 77 3.14 Computed return loss for the CPW series stub terminated with artificial absorbers (open structure). Geometrical factors: G = 0.225mm, S = 0.45mm, L = 1.35mm, W = 0.225mm, D = 1.775mm and T = 2.0mm. The thickness and the dielectric constant of the substrate are H = 0.635mm and ~r = 9.9. The absorbers at the end are chosen as: E1 = 1.0 + j15.0 and i1 = 1.0 + j15.0 for the air side and C2 = 9.9 + j20.0 and I2 = 1.0 + j2.02 for the dielectric side....................... 78 3.15 Computed reflection coefficient of the artificial absorber in a waveguide designed for absorption of the dominant TE1o mode as a function of mesh size (6x, by, and 6z) in the absorber region. Geometrical and material parameters of the waveguide: a = 47.6 mm, b = 22.0 mm, e61 = 1.0, and /,r1 = 1.0. The thickness of the absorber, 6t, is 0.041Ag (5.0 mm) at 4.0 GHz. fTE1o, fTE20, and fTEol are the cutoff frequency of each mode.... 81 3.16 Computed transmission coefficient for a conducting cube of dimensions dl = d2 = 9.0 mm and d3 = 6.0 mm. The width and height of the waveguide: a = 22.86 mm and b = 10.16 mm. The thickness of the absorber at the end of the waveguide is fixed to 5.0 mm which is 0.06Xg at 8.0 GHz..... 82 4.1 Geometry of a microstrip short circuit via hole at the center of the line. The dashed lines represent PEC walls to terminate FEM meshes. W1 = 85am, W2 = 600tim, W3 = 1785pum, H1 = 250pm, H2 = 600kim, 1e = 12.9, and e2 = 1. The dotted lines around the geometry represent the outer boundaries for FEM mesh termination.................... 88 4.2 Scattering parameters for the via.................... 89 xi

4.3 Inductance of the microstrip short circuit with via hole. The inductance is calculated using two different methods. In FEM, the flux density around the via is used (f = 10 GHz). In FDTD, however, the via inductance is extracted by fitting the derived S-parameters to a T-equivalent circuit.. 90 4.4 Scattering parameters of the microstrip short circuit for two different substrate heights (H1 = 0.2, 0.3 mm) with via hole grounding at the center of the line.................................. 91 4.5 Magnitude of the electric field distribution at f = 20 GHz under the microstrip short circuit (z = 0.19 mm plane). The edge singularity and the standing wave patterns of the electric field distribution are clearly revealed. 92 4.6 Geometry of a 1-port microstrip ground pad. This is a board level via. Wo = 2.30 mm, W1 = 1.15 mm, W2 = 2.88 mm, W3 = 10.0 mm, H1 = 1.0 mm, H2 = 2.0 mm, e1 = 3.4, e2 = 1. The dotted lines around the geometry represent the outer boundaries for FEM mesh termination... 93 4.7 The via inductance calculated from the flux/current definition at 10 GHz as a function of via height.................................. 94 4.8 Phase of the S11 of the microstrip ground pad computed using the FEM at several different reference planes (Dy = +0.41, 0.0, -0.44, -1.44, -2.88, and -3.24 mm from the top to the bottom of the graphs). Dy is 0.0 mm at the geometrical center of the via. H1 = 1.0 mm................ 95 4.9 The via inductance of the microstrip ground pad calculated from the flux/current definition for three different via heights (H1 = 0.5, 1.0, 1.5 mm) as a function of frequency........................ 96 4.10 Magnitude of the total electric field distribution at (a) f = 10 GHz, (b) f =15 GHz, and (c) f = 20 GHz under the microstrip line at z = 0.9 mm plane. The edge singularity and the standing wave patterns of the electric field distribution are clearly revealed..................... 97 4.11 Magnitude of the magnetic energy distribution on the ground plane at f = 10 GHz. The via inductance is mainly determined by the magnetic energy surrounding the via............................ 98 4.12 Geometry of the transition from CPW to microstrip on different layers without via hole. Wg = 50/m, WI = 75/n, Wo = 75/m, W, = 200um, W = 1000,m, Hi = 100/lm, H2 = H3 = 400[/m, and er = 12.9..... 99 4.13 Scattering parameters for the transition from CPW-to-microstrip on different layers without via hole connection........................ 100 Xli

4.14 Transition from CPW to microstrip on different layers through via hole. "I = 50/pm, W' = 75pm, 'Wo = 75pm, Ie = 200pm, WH = 1000tnm H1 = 100pm, H2 = H3 = 400ypm and cr = 12.9. The vertical via is an empty rectangular cavity with PEC walls.......................... 102 4.15 Scattering parameters of the CPW-to-microstrip transition with vertical via hole connection............................... 103 4.16 Geometry of a channelized CPW to microstrip transition. W' = 50[im, 14c = 75,m, W = 1000/pm, H1 = 400pm, H2 = 500pm and r = 12.9... 104 4.17 Scattering parameters for the channelized CPW-to-microstrip transition. 105 4.18 Magnitude of transverse component of the electric field at z = 390pm plane...................................... 106 4.19 Geometry of asymmetric two port via hole passing through a platedthrough-hole (PTH) in the middle of the substrate. H1 = H2 = 0.2 mm, W1 = 0.4 mm, W2 = 0.8 mm, r1 = -r2 = 4.5, D = 0.2, 0.4, 0.7 mm (3 cases)...................................... 108 4.20 Comparison of the S-parameters for 3 different PTH sizes........ 109 4.21 Via inductances of two different via height: H1 = H2 = 0.1 and 0.15mm and D = 0.4 mm. Inductances are evaluated from the magnetic flux density surrounding the via.............................. 110 4.22 Type I equivalent circuit used for the asymmetric 2-port via hole geometry. Note that if the PTH opening becomes smaller, C increases and finally becomes a short circuit............................ 111 4.23 Matching the S parameters using the type I equivalent circuit (Circle: FEM data, Line: Equivalent circuit). H1 = H2 = 0.1 mm and D = 0.2 mm. Lvia = 25 pH, C = 0.28 pF with type I equivalent circuit......... 112 4.24 Geometry of double via hole. H = 1.6 mm, W1 = 1.4 mm, W2 = 1.0 mm (85 F), W3 = 0.2 mm, W4 = 1.0 mm, W5 = 1.4 mm, 1rl = 4.5, D = 0.2, 0.4, 0.6, 0.8 mm (4 cases)........................... 113 4.25 Comparison of the S-parameters for 4 different separation distances.... 114 4.26 Magnitude of the total electric field distribution at (a) f = 4 GHz and (b) f = 10 GHz under the microstrip line at z = 1.5 mm plane. The edge singularity and the standing wave patterns of the electric field distribution are clearly revealed............................ 115 xiii

4.27 Type II equivalent circuit used for the double via hole geometry. The L12 and C12 represent the mutual coupling between the via holes. Note that if the separation distance D increases, the L becomes smaller and finally becomes ideal short circuit............................... 116 4.28 Matching the S parameters with the type II equivalent circuit. Lvia = 380 pH, L12 = 80 pH, and C12 = 0.1 pF. H = 1.6 mm, 14' = 1.4 mm, WI2 = 1.0 mm, W3 = 0.2 mm, 44 = 1.0 mm, 14W = 1.4 mm, Crl = 4.5, D = 0.8 mm. Circle: FEM data, Line: Equivalent circuit.............. 1.17 4.29 Geometry of the CPW-to-slotline planar transition with air bridge: L, = 4.751 mm, S1 = 0.114 mm, S2 = 0.129 mm, W1 = 0.305 mm, W2 = 0.254 mm, W3 = 0.178 mm, H = 0.635 mm, and Cr = 10.5............... 119 4.30 Scattering parameters for the coplanar waveguide-to-slotline transition with air bridge............................... 120 4.31 Geometry of the CPW-to-slotline non-planar transition: L1 = 2.54 mm, L2 = 2.54 mm, L3 = 5.0 mm, L4 = 6.229 mm, W1 = 2.54 mm, W2 = 0.152 mm, W3 = 0.318 mm, W4 = 2.54 mm, W1 = 0.178 mm, H = 0.635 mm, and Cr = 10.5.................................... 122 4.32 Scattering parameters for the coplanar waveguide-to-slotline transition under shielding environment. The design frequency is 6 GHz.............. 123 4.33 An example of a CPW-to-waveguide probe station............. 125 4.34 CPW-to-waveguide transition: a = 47.6 mm, b = 22.0 mm, L1 = 10.8 mm, L2 = 0.8 mm, W1 = 0.5 mm, W2 = 0.8 mm, W3 = 0.5 mm, W4 = 11.8 mm, Ws = 3.8 mm, h1 = h2 = 2.0 mm, h3 = 1.2 mm, h4 = 13.2 mm, and Cr = 12.0................................... 126 4.35 Scattering parameters of the CPW-to-waveguide transition.......... 127 4.36 An example of the convergence test of the FEM. The inductance of the microstrip ground pad shown in Fig.4.6 as a function of the number of cells per guided wavelength.......................... 129 5.1 Geometry of the Microstrip-to-CPW transition without hermetic wall. W1 = 0.48 mm, W2 = 0.14 mm, W3 = 0.2 mm, W4 = 2.0 mm, Ws = 3.0 mm, H = 0.635 mm, Cr = 12.5. Case 1: D = 2.26 mm. Case 2: D = 1.76 mm........................................ 136 5.2 S-parameters for the open transition shown in Fig. 5.1. In this case, duroid substrate, which has a relative dielectric constant C. of 10.8, is used. The CPW ground planes and the PEC underneath the substrate are connected in the measurements.............................. 138 xiv

5.3 Comparison between FEM and FDTD data for (a) case 1 (D = 2.26 mm) and (b) case 2 (D = 1.76 mm)............................ 139 5.4 S-parameters for two different size of CPW ground plane (14 in Fig. 5.1). Other geometrical factors are remain the same.................... 141 5.5 Geometry of the microstrip-through-CPW transition without hermetic wall: vertical wire grounding. W1 = 0.466 mm, W2 = 2.034 mm, D1 = 0.40 mm, and D2 = 0.40 mm. Other geometrical factors are the same as those in the previous example........................... 142 5.6 S-parameters of the vertical wire grounding case............... 143 5.7 Electric field (E,) distribution around the transition region at f=10 GHz. (a) Microstrip-through-CPW transition, (b) Microstrip-through-CPW transition with vertical grounding wires, and (c) Microstrip-through-CPW transition with extended rectangular vias.................... 144 5.8 Geometry of the microstrip-through-CPW transition without hermetic wall: Rectangular via hole grounding with finite size CPW ground plane. The size of the via hole is 0.452 x 0.6 mm. WI = 0.226 mm, W2 = 0.452 mm, W3= 0.452 mm, W4 = 1.13 mm. Other geometrical factors are identical to the first example.............................. 145 5.9 S-parameters for the rectangular via hole grounding case........... 146 5.10 Geometry of the microstrip-through-CPW transition without hermetic wall: extended rectangular via hole grounding. The top view shows one quarter of the whole geometry. W1 = 0.1 mm, W2 = 0.14 mm, W3 = 0.226 mm, W4 = 0.2 mm, Ws = 0.6 mm, W6 = 0.2 mm, W7 = 1.5 mm, Ws = 0.5 mm, W9 == 2.034 mm............................... 147 5.11 S-parameters of the extended rectangular via hole grounding case..... 148 5.12 Geometry of the microstrip-through-CPW transition without hermetic wall: PEC wall grounding. W1 = 2.034 mm and other geometrical factors are identical to the first example............................... 149 5.13 S-parameters of the PEC wall grounding case................. 150 5.14 Geometry of the microstrip-through-CPW transition with hermetic wall on top of the extended PEC via grounding case. The width of the hermetic wall is 1.6 mm. W1 = 0.48 mm, W2 = 0.14 mm, W3 = 0.2 mm, W4 = 1.6 mm, H1 = H2 = 0.635 mm, Er1 = 12.5, and Er2 = 2.3............... 151 xv

5.15 Comparison of the S-parameters for the transition with hermetic wall on top and extended PEC via under the circuit plane...............152 5.16 S-parameters for the transition with and without hermetic wall on top of the CPW section. The geometrical factors: 14i = 70 /tm, W2 = 30 /m, W3= 40 iim, W4 = 300 tm, and I5 = 350 um (refer to Figure 2). The width and height of the dielectric brick under the PEC wall are 250 /m and 100 um, respectively........................... 153 5.17 Geometry of the hermetic bead transition. W1 = 0.55 mm, W2 = 0.21 mm, W3 = 1.27 mm, 14 = 2.225 mm, 1Vs = 5.0 mm, H1 = 0.635 mm, H2 = 4.0 mm, L1 = 1.50 mm, L2 = 0.40 mm, and c, = 10.8............ 155 5.18 S-parameters of the hermetic bead transition with square dielectric ring. 156 5.19 Comparison of the S-parameters for the hermetic bead transition of three different air gap width L2. (a) IS111 and (b) IS21.............. 157 5.20 Comparison of the S-parameters of the hermetic bead transition with square and circular dielectric ring of radius equals H1............ 158 6.1 Top view of the K/Ka-band MMIC package.................. 161 6.2 Schematic diagram of the K/Ka-band hermetic package designed and manufactured for a MMIC phase shifter chip. The diameter of the vertical metal filled vias are D, = 0.203 mm and the distance between these vias are D, = 1.016 mm. The upper and lower alumina layers are each 0.381 mm thick (Dt)............................... 163 6.3 Schematic showing the location and extent of PECs in addition to the artificial absorbers which simulate the test fixture and open environment. Two different types of absorbers (i = a = 1.0 + j10.0, del = 9 5 + j15.0, and djit = 1.0 + jl.5789) are designed and placed in either side of the package facing the DC bias lines (not fully shown in the figure for simplicity)..................................... 165 6.4 Computed S-parameters for the isolated symmetric package......... 166 6.5 Computed S-parameters for the isolated asymmetric package......... 167 6.6 Computed vertical electric field distribution (dB scale) at (a) f = 20.5 GHz and (b) f = 29.5 GHz in the isolated asymmetric hermetic package.... 168 6.7 Comparison between the computed S-parameters for the isolated asymmetric package with (original) and without (simplified) DC bias lines and bond wires.................................. 169 xvi

6.8 Measured and computed S-parameters ( (a) |IS11 and (b) IS211 ) for the asymmetric package residing in the test fixture...................... 171 6.9 Computed vertical electric field distribution (dB scale) in the asymmetric hermetic package placed in the test fixture at (a) f = 20.5 GHz and (b) f = 29.5 GHz in the isolated asymmetric hermetic package............ 172 6.10 Schematic diagram of the asymmetric hermetic package showing four additional vias at the input and output ends of the package......... 173 6.11 Computed S-parameters for the asymmetric package with four additional vias at the input and output ends of the package................. 174 6.12 Computed S-parameters ( (a) IS1l and (b) 1S211 ) for the asymmetric package residing in the test fixture with resistive coating on the inside surface of the lid. The resistivity and the thickness of the coating is indicated as R in [Q/O] and t in [mm], respectively............................ 176 6.13 Computed S-parameters for the asymmetric package with frame constructed from dielectric material with high loss tangent. Case 1: er = 9.5(1.0+j0.1) and r- = 1.0 + jO.0. Case 2: cr = 9.5(1.0 + jO.5) and /, = 1.0 + j0.0... 177 7.1 Traditional way of modeling the inductances of the PCBs and package bonding wires................................. 183 7.2 Modeling of simultaneous switching noise (SSN) and simulation methodology...................................... 186 7.3 Equivalent circuits for three different tiles, such as power/ground tile, power/ground pin tile, and power/signal/ground tile............... 189 7.4 Example PCB used in the validation experiments................. 195 7.5 (a) Board layout and current waveform for validation experiments. Noise maps on the (b) power plane and (c) ground plane for three different number of tiles................................ 197 7.6 CPU and memory usage versus the number of tiles................ 198 7.7 Variation in the position of the connector................. 198 7.8 Noise map (peak-to-peak voltage variations) on the power and ground planes with connector in different place: (a) left, (b) top, (c) right, (d) bottom............................................... 199 xvii

7.9 Noise maps (minimum on power and maximum on ground plane) (a) without and (b) with decoupling capacitors. Note that the scale of the vertical axis is not identical for both cases........................ 201 7.10 Configuration of the decoupling capacitance......................202 7.11 Noise variation on power and ground planes with various decoupling capacitances.........................................202 7.12 Top and cross sectional views (not to scale), and source/stripline configuration. The parenthesized numbers are coordinates in [mm] relative to lower left corner of the PCB. V denotes connection to ideal ground, i.e., node 0 of HSPICE.................................... 204 7.13 Noise maps on the power and ground planes (a) with signal trace and (b) without signal trace............................ 205 7.14 Split PCB and short circuit pins....................... 206 7.15 Gap in a microstrip line (a) and its lumped equivalent circuit (b)...... 206 7.16 Configuration of the split PCB having 25 short circuit pins ('@') and 11 current sources (Imax = 100 mA and 0.5 nsec rise/fall time). The dots ('.') represent electrical nodes on the planes. Note that the connector is placed at the left hand side of the board supplying 0.0 V only for the lower plane. 208 7.17 Noise maps on the (a) upper (Vmax = 420 mV) and (b) lower planes (Vmax = 260 mV).............................. 208 7.18 Schematics of the Intel P6 gadget........................ 210 7.19 Input current waveform (0.8 nsec rise and fall time) and its frequency component............................................ 211 7.20 Output current waveforms at 4 different points with 64 input current sources switching simultaneously: magnitude Imax = 15 mA and rise/fall time tr = tf = 0.8 nsec............................ 212 C.1 Vertical via hole as a current path........................ 227 C.2 Equivalent circuits for one- and two-port network.................. 228 C.3 Symmetric two port network with the even and odd mode excitations corresponding to placing magnetic and electric walls at the center of the circuit.230 xviii

LIST OF TABLES Table 2.1 Summary of FLoating point OPeration (FLOP) counts for reduced equivalent problems. In the table, the number of iterations for convergence are assumed to be the maximum which is the matrix size............... 28 2.2 Number of zero eigenvalues for simply connected cavity with E-field formulation............................................. 33 2.3 Number of zero eigenvalues for simply connected cavity with H-field formulation................................... 33 2.4 Two types of tetrahedral labeling systems.................. 36 2.5 Estimations of some of the important quantities in the structural tetrahedral mesh. It is assumed that the domain is simply connected and there are no interior PEC surfaces which will reduce the total edge counts. N, Ny and Nz are the number of nodes in x, y and z direction, respectively.. 38 3.1 First few eigenvalues of an empty rectangular cavity having a = b = c = 10 mm............................................... 56 3.2 First few eigenvalues of an empty rectangular cavity having a = 40, b = 30, and c = 20 mm.................................... 56 3.3 First few eigenvalues of a partially-filled rectangular cavity having a = 40, b = 30, c:= 20, c' = 5 mm, and c,, = 2.5......................... 58 3.4 First few eigenvalues of a partially-filled rectangular cavity having a = 40, b = 30, c = 20, c' = 5 mm, and e, = 12.9.................... 58 3.5 Resonant frequencies and Q-factors for lossy dielectric loaded cavities having various amount of losses............................... 62 4.1 The equivalent circuit parameters of the double via hole structure and corresponding resonance frequency. The via inductance Lvia is extracted from the flux model, while the shunt terms, L12 and C12, are determined to match the FEM data with the equivalent circuit................. 118 xix

LIST OF APPENDICES Appendix A. Calculation of the Element Matrix Entries........................... 220 B. Pre-conditioned Conjugate Orthogonal Conjugate Gradient Method...... 223 C. Network Parameters............................................ 226 C.1 Inductance of Vertical Via Hole............................. 226 C.2 One-Port Network....................................... 228 C.3 Symmetric Two-Port Network.............................. 229 C.4 Asymmetric Two-Port Network......................... 231 xx

CHAPTER I INTRODUCTION 1.1 Motivation and Background In the recent advancement of satellite, mobile and personal communication systems, military and commercial radar systems, remote sensing technologies, and numerous other areas, high-speed high-frequency MIMICs (monolithic microwave integrated circuits) play an essential role in terms of small volume, low weight, low power consumption, and low cost. In spite of these advantages of using MMICs, their design and optimization are not an easy task due to their complexity and sophisticated functionalities, especially as the operating frequency increases up to the X-band. In many cases, both the microwave and millimeter-wave subsystem with MMICs chips and the overall system are very complex and interrelated, and the functions of the circuits and interconnects are very sensitive to their mechanical and electrical environment. Depending on the system, there exist many different interconnection structures, ranging from simple through lines to whole packages, which need to be characterized and fully understood. Furthermore. the advent of high-speed computers and their subsidiary systems require high-speed operations reaching Gigabits per second switching times and data transfer, and the resulting electromagnetic interference and spurious radiation arising from the circuit discontinuities, in many cases, 1

2 cause serious malfunction of the system if not carefully designed. For these reasons, the successful development of efficient design and analysis tools has been an essential step toward advances in these technologies. For low frequency lumped electric circuits and solid state devices, there are numerous circuit simulators providing accurate time domain characterization of digital and analog circuits. Similarly, design and optimization software is also available for canonical planar circuits in the microwave frequency region. Existing electric circuit simulators and microwave circuit design tools, however, have their own limitations, including lack of modeling capability for electromagnetic interactions between closely or sometimes remotely placed circuit components. In the case of microwave and millimeter-wave circuit design and characterization tools, in particular, when there are material inhomogeneities and three-dimensional structural discontinuities, the limitations of those tools become apparent. Recently, to overcome the shortcomings of the available simulation tools a number of full-wave analysis methods have been developed based on the integral representation of Maxwell's equations [1] or on partial differential operators [2, 3, 4]. The method of moment (MoM) and mode matching techniques are examples of the former approach, while the finite element method (FEM) and the finite difference tinme domain (FDTD) method are most frequently encountered as examples of the latter approach. The aforementioned full-wave techniques are considered to provide complete solutions to complicated structures by solving Maxwell's equations rigorously under proper boundary conditions. The MoM, in general, provides accurate and efficient solutions, and, thus, is recommended for use whenever applicable. In practical applications, however, the MoM is mainly employed for planar two-dimensional and fairly

3 simple three-dimensional circuit characterizations due to the difficulties in obtaining relevant Green's functions for non-canonical geometries. In contrast to the MoM, the partial differential equation (PDE) techniques, such as the FEM and FDTD methods, do not require the knowledge of Green's functions for a given problem, thus providing flexible modeling capability regardless of two and three-dimensional discontinuities. In principle, the FEM and FDTD methods are similar to the MoM in terms of using basis functions on discretization space and performing weighting operations. In the FDTD method [5, 6], however, the pulse functions are used in both the basis and weighting function spaces, and as a result the FDTD method requires discretization of the whole three-dimensional space with a dense sampling rate (at least 15 samples/Ag) due to its simple basis and weighting functions. On the other hand, FEM employs similar discretization and weighting operations to find weak solutions of Maxwell's equations, except for the fact that any arbitrary functions of higher order can be used as basis and weighting functions [2, 3, 7]. With higher order basis and weighting functions, the same accuracy can be achieved with less discretization, because of their modeling efficiency for highly varying electromagnetic fields. Another important feature of the FEM is the usage of non-uniform discretizations for adaptive modeling of rapidly changing local electromagnetic field distributions, thus allowing optimum mesh or sub-domain sizes. The meshing flexibility of the FEM also contributes to reduced overall problem size and faster solution time, and poses clear advantages over the regular FDTD method. In addition, recent developments in high-speed parallel computing facilities including shared and distributed memory machines have an impact on practical application and real time computation of numerical techniques.

4 In view of the above discussion, the frequency domain finite element method is developed and parallelized in this study for the characterization of arbitrary-shaped three-dimensional high-speed, high-frequency interconnects. In the following, the objectives of the study and the overview of this thesis will be described in detail. 1.2 Objectives The goal of the present dissertation has been to develop and apply an accurate and fast modeling method for the characterization of high-speed, high-frequency circuits and interconnects, and to explore the three-dimensional finite element method using two fundamental features, tetrahedral sub-domain elements and their corresponding edge-based vector basis functions. While the conventional node-based finite element method has been plagued by spurious modes and requires special treatment to obtain physically correct solutions [8]-[12], the edge-based finite element method provides spurious-free solutions by enforcing correct boundary conditions at the material interface and sharp metal discontinuities [13]-[19], making the approach attractive. As a result, the edge-based finite element solution scheme has been developed to provide reliable solutions for very complicated structures without using any further post-processing or sorting of the solution. Furthermore, to overcome the difficulties in applying time-consuming full-wave techniques to the real-time characterization of a given circuit or interconnect, the development of efficient parallelization strategies is desirable. Therefore, the frequency domain FEM has been parallelized on distributed memory parallel computers, such as the IBM SP2, using the message passing paradigm [20]. The linear scalability of the FEM code is rooted in the frequency independent nature of the method, which requires minimal inter-processor communication overhead. The parallelized FEM

5 can provide near-real-time characterization of a given circuit and in turn opens the door to the use of the full-wave techniques for circuit design, as well as tools for optimization. This developed tool has been applied for the characterization of various highspeed high-frequency interconnects and the theoretical results have been compared with measurements or other theoretical techniques whenever possible for validation. Accurate characterizations of the aforementioned individual interconnects are of critical importance for successful design of more complicated circuits. Furthermore, to tackle open boundary problems often encountered in real situations, design equations for isotropic absorbing layers have also been sought, and their validity and effects have been fully tested under several different circumstances. The characterization of high-speed digital circuits and packages in the frequency as well as the time domain is an important but difficult task, especially with full-wave techniques, due to their heavy computational burdens. To overcome this limitation, a hybrid approach has been explored in this study which combines the frequency domain full-wave FEM and time domain circuit simulator, HSPICE [21]-[25]. In this approach, various digital ICs and packages are discretized into a number of smaller pieces, called tiles, and for each tile an equivalent circuit is derived using FEM and a simpler microwave network theory. After all the equivalent circuit extraction procedures have been completed, time domain circuit simulations are performed to model the whole structure, and develop noise maps and time domain potential fluctuations plots. The information obtained from the system level modeling is essential for diagnosis of simultaneous switching noise (SSN) and for effective suppression of the noise arising from the high-speed digital circuits.

6 1.3 Overview This thesis is focused on the development of a rigorous theoretical technique, namely the FEM, for the characterization of three-dimensional high-speed highfrequency interconnects and its applications to various geometries. Each geometry has its own unique characteristics and complexities which have prevented them from being fully characterized with conventional approach. In Chapter 2, the theoretical framework of the three-dimensional finite element method is briefly summarized. Unlike the conventional approach to FEM theory, the Hilbert space representation of the method is adopted starting from the weak form of the wave equations. Since appropriate boundary conditions should be specified to obtain a unique solution of the wave equations, the most fundamental boundary and interface conditions often encountered in the practical problems are described. Furthermore, in view of the importance of using the edge basis functions, as described in the previous section, for obtaining spurious-free solutions, the rationale of using the Whitney edge basis functions is fully discussed. In addition, excitation mechanisms in the FEM system and customized tetrahedral mesh generation schemes are presented. For the simulation of unbounded problems which are in many cases unavoidable, isotropic lossy materials are utilized to simulate the open boundary conditions. The last part of this chapter contains the theoretical background of the idea of incorporating passive lumped elements into the FEM framework. Even though an extensive study of this method is not presented in this dissertation, a preliminary study shows its effectiveness and accuracy. In Chapter 3, several different eigenvalue problems are tackled to validate the developed FEM and to assess its accuracy. For a rectangular cavity having lossy

( dielectric material, the linear system of equations incorporating an excitation mechanism is solved, and the resonant frequency and the quality factor are computed. From the above information, the real and imaginary parts of the dielectric material can be extracted easily. After the FEM is validated, it is also parallelized on a distributed memory parallel computer (the IBM SP2) and two different parallelization strategies are examined. With the "task parallelization" strategy, in particular, almost linearly scalable performance improvement has been obtained, which is quite remarkable. Furthermore, the design concept of the artificial absorbers proposed in Chapter 2 is examined through several different circuit configurations supporting quasi-TEM as well as non-TEM waves. Starting from Chapter 4, applications of the developed FEM to various practical high-speed, high-frequency interconnects are presented. In Chapter 4, eight different types of basic building blocks of the MMICs are characterized [26, 27, 28] and in some cases the equivalent circuits are obtained [29, 30]. In particular, the equivalent circuits for via hole geometries are used in Chapter 7 to provide appropriate equivalent circuits for the power and ground pin tiles. In Chapter 5, the performance of intra- and inter-chip hermetic transition structures are fully characterized to provide design guidelines to MMIC designers [31, 32]. Specifically, the microstrip-through-CPW configuration is employed for an intra-chip transition, while hermetic bead structure is utilized for inter-chip transition. The effect of various grounding via holes, cavity size, and the shape of the hermetic beads are carefully examined. In Chapter 6, the K/Ka-band (18 to 40 GHz) MMIC package fabricated for a phase shifter is characterized and the effects of the structural symmetry/asymmetry, bonding wires, gaps, side walls formed by 12 via holes, and the DC bias lines are

8 also studied [33, 34]. Furthermore, the spurious package resonances are identified and several techniques for suppressing them are suggested with rigorous numerical data. To the best of the author's knowledge this is the first comprehensive study of RF energy leakage and internal resonances in a practical millimeter-wave package. In Chapter 7, the detailed procedure for the system level electromagnetic modeling of high-speed digital ICs and packages is presented. In the tiling procedure, three types of primitive tiles are identified and their equivalent circuits are devised, including power/ground plane tile, power/signal/ground plane tile, and power/ground pin tile. Upon building the equivalent circuit library, several different printed circuit boards and source configurations are studied along with an Intel test circuit for validation of the proposed approach. The simulation and measurement results show very good agreement. Considering that there are 8 layers in the PCB and hundreds of vertical pins and traces, the remarkably good comparison between the measurement and the simulation proves the accuracy and effectiveness of the proposed hybrid approach. Finally, in Chapter 8, a brief summary of the achievements and original contributions of the present research are presented. Also, concluding remarks and directions for future studies are provided for later study.

CHAPTER II THEORETICAL METHODOLOGY: 3D VECTOR FINITE ELEMENT METHOD 2.1 Introduction The finite element method (FEM) has been extensively applied to scattering, radiation and propagation problems [35]-[40]. Recently, it has found application to circuit problems as well and has led to successful treatment of microwave and millimeter wave circuits with complex geometrical shapes and inhomogeneous material combinations [33]-1[30]. In this chapter, the basic principles and general procedures of the three-dimensional (3D) vector finite element method are presented starting from Maxwell's equations and subsequent vector wave equations. The 3D vector FEM is based on the discretization of the weak form of the wave equation in frequency domain, and provides accurate full wave solutions to the most general electromagnetic problems. In this study, even though both electric and magnetic field formulations are given for the sake of completeness, only the electric field formulation is implemented for several practical reasons, such as ease of application of boundary conditions and reduction of overall problem size. The current formulation presented is valid for geometries having isotropic dielectric materials whether lossy, lossless or a combination and perfect conductors. However, it can be easily extended to incorporate anisotropic 9

10 materials and lossy conducting surfaces. The FEMN developed in the study is capable of handling closed as well as open boundary problems. In the following sections, the essential formulation of the finite element method and its discretization procedure will be discussed in detail. 2.2 Vector Wave Equation Consider an arbitrary shaped three dimensional geometry which contains material inhomogeneities and geometrical discontinuities, as shown in Figure 2.1. The electromagnetic fields in the domain Q are characterized by Maxwell's equations aB V x E = -t -M (2.1) V x H = t +J (2.2) V D = Pe (2.3) V B = pm (2.4) and the appropriate constitutive relations D = cE (2.5) B = puH (2.6) where E and H are electric and magnetic field intensity in [Volt/m] and [Amp/m], respectively and similarly, D and B represent electric and magnetic flux density in [Coul/m2] and [Weber/m2], and the sources and material parameters are defined as: Ji: impressed electric current in [Amp/m2] Mi: impressed magnetic current in [Volt/m2] p,: electric charge density in [Coul/m3] pm ' magnetic charge density in [Weber/m3]

11 (E, p) Ji (E, H) - ml l. "s' — / PEC, PMC or/and ABC Figure 2.1: Geometry of three dimensional volume which has current sources J; and/or M; with possible material inhomogeneities and geometrical discontinuities. The boundary of the volume V is closed by some combination of PEC, PMC, and ABC. e:= CoCr ~ permittivity in [F/m] -' =: olr: permeability in [H/m] Without loss of generality, it is assumed that fields have time harmonic variations (ewt ) and the impressed sources are not present in the computational domain. Then the source-free vector wave equations can be derived from Maxwell's equations: 1(Vx E) E re V x -E) k 0 (2.7) L61(V x H)!rH where ko (= wVVJ0) is the wave number in the medium. Later sections will address the fact that the circuits can be excited by using the surface integral term which is equivalent to the impressed sources. For a numerical solution of the electromagnetic fields in a given domain, the wave equations, which are second order partial differential equations, need to be solved with appropriate boundary conditions for a unique solution. However, since finding

12 a strong solution of the wave equations is not always possible, weak solutions will be sought in this study with appropriate weighting and testing operations. Among the two equivalent formulations, E- and H-field formulation, the E-field formulation will be presented and tackled. As mentioned before, the E-field formulation has clear advantages over the H-field formulation when there are many perfect electrically conducting (PEC) surfaces or bodies. Moreover, with the E-field formulation, the implementation of the boundary condition for a perfect magnetically conducting surface, which is the natural boundary condition, is also straightforward and will be discussed in the next section. 2.3 Weak Form of Wave Equation and Its Discretization As a first step towards finding the finite element solution, a given solution space Q is subdivided into number of smaller sub-domains Q,, and in each e, the electric field is expressed as a linear combination of basis functions P = {P1, P2,,Pm,}T C } C1, a Hilbert space, 1 as follows; m E = >xiPi (2.8) i=1 where x = {xi, x2,...,xm} C Cm, m-vectors of complex coefficients, is the unknown coefficient vector to be determined. Note that due to the nature of approximation of the discretization procedure, the original problem is modified to become a new problem which has a different solution space [43]. After substituting Eq. (2.8) into 'Let (7i, (,.)) be an inner-product space. If the associated normed linear space (}l, 1 11II) is complete, then (7-, (,.)) is called a Hilbert space. The inner product (.,.) on a linear space Ri satisfies (a) (v,v) > O Vv E X (b) (v,v) = 0 v = 0.

13 Eq. (2.7), we obtain discretized wave equation: m x, ( X i( V XV Pi - k2erPi) ) = 0. (2.9) i=1 where the inner product (A,B) is defined as follows: Th aid of some ofnext stepctor is the discretizatihe range odifferential operatorem3, Eq. (2.10) V X (UrVx) - k 2Er, and defining the weighting functions Q = {Q1, Q2,... I Qn.}T C 7i 2, a Hilbert space. Now, the weighting operation on the discretized wave equation is performed as m E X-(Q, (XV X (rlV x P i, - k x Q)) = for j=l,...,n (2.10) i=1 (A, B) = lit A B dv. Now, with aid of some of vector identities2 and the divergence theorem3, Eq. (2.10) can be written as m E ^{(- 1V x P., V x Q) k 2o(Pi,, Qi)} i=1 m + X I j n. -(VXPi ) X Q]ds=O for j=l,...,n. (2.11) where the last term can be modified to accommodate an equivalent surface current as -ja"Io j (f x H). Q ds. (2.12) In view of the continuity of the magnetic field across any material interface bearing no surface current, it is recognized that the contribution from the last term of Eq. (2.11) for interior Qe is zero except the cases when there are surface current sources. 2B*V xA= V (A x B) +A V xB 3 V. AdQ = A. dS

14 The linear system of equation (2.11) can be represented in a matrix form for each subregion Qe as follows; [Ae]T [Xe] = [be] (2.13) where [Ae] = [Ue]-k2 [Ve] MUe = (1V x P,, V x Q j) Vi = (C rPir~ j) Xi [X,...,xam e E [bt...b]-=-jpo (n x H).Qds JAe and 1 < i < m, 1 < j < n. The column vector [be] in the right hand side of Eq. (2.13) can be used as a forcing function to impress an arbitrary surface current. 2.4 Basis and Weighting Functions Given the limitation of the node-based FEM, such as incorrect modeling of boundary conditions, and solutions contaminated by unwanted spurious modes, we applied in this study the e(lge-based vector basis functions We with tetrahedral sub-domain elements [13]-[19], which are known as Whitney 1-form,4 as both basis and weighting functions (P- = W~ and Qj = Wj) leading to Galerkin's method.5 In terms of global coordinate systems the edge basis function We is defined as follows (refer to Fig. 2.2) for j = 1,...,6, f + gj x r, r E, WJ - (2.14) 0, r e 4See section 2.7.2. 5Galerkin's method is known to provide second-order accuracy [44].

15 i~ YE3 E2 x Figure 2.2: Tetrahedral element and associated nodes and edges definitions. The circled numbers 1 to 4 represent four node labels and El to E6 six edges associated with those nodes. ri1 and ri2 are node vectors associated with the ith edge and r is an arbitrary vector which lies in the tetrahedron.

16 where bj fj = 6 r(7-)l x r(7-j)2 bjb7-j 93 6Qe gj = - e7-j = (r(7-j )2- r(7-/)l)/bj b7-3 = I r(7-j - r(7-)l and e, is the unit vector for the jth edge and bj is the length of the edge. The r(7-j) and r(7-j)2 are the vectors from the origin to the two nodes consisting of the (7 - j)th edge. Alternatively, the definition of the edge basis function with respect to local coordinates has a quite different form. For an edge p = {i, j}, We = A-VA- - AJVA, (2.15) where An(r) (n = i or j node) is the barycentric weight of r in its tetrahedron with respect to vertex n and forms a continuous piecewise-linear function.6 The behavior of the basis function Wp can be visualized as shown in Fig. 2.3, that is, the fields turn around the axis k -1 ("central axis") remaining normal to planes containing the central axis. Also, the field has a magnitude proportional to the distance to the axis. Note that VA3 and VAs are orthogonal to facets {i, k, l} and {j, k, }, respectively. 6For example, Ai is defined in a tetrahedron {i,j,k,l} in Fig. 2.3 as Ai = volume{ir,'k.} and can be evaluated as A 1 if x coincides with the ith node i 0 if x coincides with the jth(j $ i) nodes. For an arbitrary point r = (x, y, z) e, the relationship between the global and local coordinate systems can be found as X X1 X2 X3 X4 A1 y = yl Y2 y3 y4 A2 Z ZI Z2 Z3 Z4 A3 1 I 1 1 1 A4 where (xi, yi, zi) denotes the global coordinates of the 4 vertices of a tetrahedron.

17 1 i,> -Vki kiVxj- jV^i < 1 Figure 2.3: Behavior of edge element We for edge p = {i,j}. This basis function has several other important properties as follows: V-Wp =0 (2.16) V x W; = 2gp = 2VAi x VAj (2.17) W. eq Ir on the pth edge= 6Pq (2.18) where Jpq is the Kronecker delta function. The basis functions satisfy "the divergencefree condition" since the E field obtained from a linear combination of the basis functions exactly satisfy V-E = 0. More importantly, the edge basis function enforces the tangential continuity of the E fields across the material interfaces and allows discontinuity in the normal component across the interface, which render correct modeling of the physical boundary conditions. In addition, the modeling capability of the edge basis functions for the null-space of the curl operator contribute to the fact that the solutions are not contaminated with spurious modes [43]-[47], a serious drawback of the node based FEM. Equation (2.17) is interpreted as meaning that the curl of the basis function is a constant vector and provides a convenient tool for

18 calculation of the H field from a given E field as 6 H= - Xg, (2.19) However, the H field computed above has a constant value in each tetrahedral element and thus provides limited accuracy in modeling the magnetic fields. The third property, Eq. (2.18), implies that We has no tangential components along any of the edges except pth one. Therefore, coefficient xz can be considered as a projection of the electric field onto the ith edge. This property can be easily verified by considering the facts that A = 0 on the edges {j,k},, {j, l}, and {1, k} and Aj = 0 on the edges {i, k}, {i, 1}, and {1, k}, and the fact that We has no tangential component on the facets {j, k, } and {i, k, }. 2.5 Assembly of Matrix Equation and Its Solution Methods From the sub-domain matrix equation, [Ae]T[x] = [be], the global matrix equation is constructed by assembling the contributions from all the tetrahedrons. In practice, a unique global edge numbering system, global-to-local indexing, is introduced to every edge and the augmenting procedure is performed based upon this edge numbering system. In other words, for every tetrahedron a 6 x 6 element matrix is created according to the local numbering system and a lookup table, I(e,i), for global edge numbering system. The lookup table relates the local edge number, i, on a particular element, e, to its position in the global data structure. To assemble the global matrix equation, each row and column in [Ae] is added to the global matrix [A] by referring to the lookup table, and a similar procedure is applied to [xe] and [b6] to form [x] and [b], respectively. As a result, the final matrix equation is assembled as: [A] [x] = [b] (2.20)

19 where M [A] = U{[U6] - [V] e=l M [XI = IX] e=1l M [b] = U[b]. e=l The M denotes the total number of tetrahedrons in the domain Q, which is, in general, different from the total number of edges, say L. The set operation UeM represents assembly operations based upon the global edge numbering system for all M tetrahedrons. Therefore, the matrices [A], [x], and [b] become L x L, L x 1, and L x 1 matrices, respectively. With the tetrahedral edge basis functions, Eq. (2.14) or (2.15), the matrix elements for [Ue] and [Ve] are evaluated as (see Appendix A for details): Uij = ( V x W, V xWj) (2.21) Vi. = (CrWi, Wj). (2.22) It is crucial to be aware of the properties of the final FEM matrix [A]. For lossy or lossless materials and even for uniaxial anisotropic materials, the ith and jth edges in a tetrahedron have identical contributions to [Ae] and eventually to [A] rendering the matrix [A] complex symmetric, but not Hermitian. Also, it is observed that the set of eigenvalues of the source-free equation7 contains multiple zero eigenvalues corresponding to the solution of the DC gradient electromagnetic fields. As a result, the matrix [A] becomes positive semi-definite system. Another very important and fundamental property of the matrix [A] is its sparsity, which is a natural consequence of the compact supported expansion and weighting functions. 7The eigenvalues are the solution to the generalized eigenvalue problem: [U][x] = A[V][x]

20 The degree of the sparsity varies case by case, but for structured mesh it can be predicted, as will be presented in later section. In this work, the conjugate (CG) and bi-conjugate gradient (BiCG) method [48] — [57] for real symmetric or Hermitian matrices and the conjugate orthogonal conjugate gradient (COCG) method [58] for complex symmetric matrix have been implemented with diagonal preconditioning to speed up the convergence rate and to improve the solution accuracy (refer to Appendix B for detail). In general, the CG type linear equation solver requires 0(N2) floating point operations in each iterattion for the full N x N matrix [A].. As a result, if it takes n iterations to solve a full matrix equation, [A][x] = [b], the required operation count is 0(nN2), which could be 0(N3) for large n comparable to N. In this case, the CG type solver loses its attractiveness and feasibility for solving huge matrix equations. For solving hugethe matrix equation generated by the FEM is very sparse and its sparsity 8 reaches higher than 99% in the usual practical problems. This sparsity gives a real attraction to the CG type solver, since the required operation counts become 0(nmN) < 0(nN2), where m is the average number of non-zero entries in a sparse matrix 9. As described in detail in section 2.9, m becomes at most 20 for structured mesh and N is several orders of magnitude greater than m. In particular, the BiCG and COCG methods are computationally efficient for positive semi-definite matrices since these methods require only 1 matrix-vector multiplication in each iteration, compared to 2 in the CG case. Even though the CG method guarantees its monotonic convergence for non-singular matrix regardless of whether 8The definition of the sparsity is given: N2 - M sparsity:= N 2-: x 100 [%], for N x N matrix where Mo is the total number of non-zero matrix entries. 9The definition of the average number of non-zero entries in a matrix is given as m:= Mo/N

21 the matrix is positive definite or not, in this study, the BiCG and COCG methods with diagonal pre-conditioning are applied in view of their performance which is superior to the CG method. These methods reveal faster convergence (fewer iterations) and require a reduced computational burden in almost all cases. 2.6 Treatment of Boundary and Symmetry Conditions 2.6.1 Boundary Conditions To obtain a unique solution of the wave equation in a given space (Q), appropriate boundary conditions are required on the surface (90Q) enclosing the volume Q. There are two important boundary conditions encountered in the electromagnetics applications, namely perfect electric conducting (PEC) and perfect magnetic conducting (PMC) boundaries. On the PEC (PMC) surface, the tangential electric (magnetic) field should be zero. The mathematical form of these conditions can be written as: n xE = 0 (PEC) (2.23) nx H = 0 (PMC) (2.24) where ni is the unit normal vector pointing outward from the surface (refer to Fig. 2.4). The implementation of the PEC or PMC boundary conditions with the edge basis functions in the context of the FEM can be done in a straightforward manner. For example, in the E field formulation, the PMC surface becomes the natural boundary having zero normal electric field component on the surface, and is implemented naturally due to the divergence-free property of the basis function, as explained in section 2.3. However, the PEC condition can be enforced by setting to zero the coefficients of the edges comprising the PEC surface, because the tangential field on a triangular surface of a tetrahedron is generated only by the three edge vectors which make up the triangle. It should be noted that when there are PEC

22 E, H O 0 E EH= OPEC/PMC surface l Figure 2.4: PEC or PMC surface S with outward unit normal vector n. surfaces, the total number of unknowns is reduced accordingly since the unknown coefficients for the edges on the surfaces are removed, producing a reduced system of equations. Another important outer boundary condition is the radiation boundary condition (RBC) describing the field behavior at infinity, as written below: lim R[V x F +jkir xF, F] = F = {E or H} (RBC) (2.25) R- oo where r is the radial unit vector and R = x2 + y2 + z2 is the distance from the origin to the observation point. The RBC demands that the electric and magnetic fields decrease as 1/fR in 3-dimensional space when R -+ oo. At the interface between two different materials, more general boundary conditions are required since the fields on both sides have non-zero values. When there are surface currents (J,) and charges (Ps) at the interface, the tangential electric field should be continuous and the normal electric flux density is allowed to be discontinuous due to ps. Similarly, at an interface bearing a surface current source the normal magnetic flux density is continuous and the tangential magnetic field should be discontinuous with the amount of J8. The boundary conditions are summarized

23 (medium 1) n El. HI + E2, H2 material (medium 2) interface Figure 2.5: Material interface S with unit normal vector n. as follows: n x (E - E2) = 0 (2.26) n (B1 -B2) = 0 (2.27) n x (H1 - H2) = Js (2.28) ni.(D1-D2) - ps (2.29) where the ni is a unit normal vector pointing from medium 2 to 1 as defined in Fig. 2.5. The surface charge density, ps and the surface current density, J,, are related to each other through the continuity equation: VS J,= (2.30) where the V.- denotes the surface divergence operator. It is worth mentioning that the equations (2.26) and (2.27) are equivalent and similarly the equation (2.28) and (2.29) are equivalent to each other. If one of the media is a perfect electric or perfect magnetic conductor, the electric and magnetic fields are zero inside and thus equations (2.23) and (2.24) are derived from equations (2.26)-(2.29). Note that no

24 Symmetry Plane Symmetry Plane I..- " 1 "- " ---- —. -^ H / / / I j\ \ \ ( ----' ---^ -------- / / / --- I PMC wall Figure 2.6: An example of a circuit configuration (microstrip line) bearing a symmetrical plane at the center and its electric and magnetic field distributions. As presented on the right, the original problem can be replaced with half of the structure, with PMC symmetry condition on one side. surface current or charges can exist on the surface of perfect conducting bodies. 2.6.2 Symmetry Conditions In general, many microwave and millimeter-wave circuits and interconnects contain geometrical and electrical symmetries, which, if handled properly, can confer computational advantages and simplicity. In this study, two types of symmetry boundary conditions (refer to Fig. 2.6 and Fig. 2.7) are exploited based on the electric and magnetic fields distributions. As illustrated in Fig. 2.6, when there are only tangential electric and normal magnetic fields components on an arbitrary plane or surface, the plane can be modeled as a PMC surface. By applying the PMC boundary condition, the original problem is

25 Symmetry Plane E Symmetry Plane H I —. --- -..- - -- -- I ____ - - - - \,, - ^! -.t^ ^ '\1^, — >j PEC wall Figure 2.7: An example of a circuit configuration (double strip line) containing symmetry at the center and its electric and magnetic field distributions. The original problem can be replaced with the half of the structure with PEC symmetry condition at the center as shown on the right. replaced with the reduced one as shown on the right hand side of the figure without changing the original field distributions and thus the solution of the original electromagnetic problem can be obtained by solving the half structure. Needless to say, it is much easier and more efficient to handle the reduced equivalent problem than the original one. Analogously, as shown in Fig. 2.7, a fictitious surface bearing only a tangential magnetic and normal electric field components can be interchanged with the PEC boundary, resulting in a reduced equivalent problem. Even though Fig. 2.6 and 2.7 illustrate the symmetries in the microstrip and double strip lines, similar symmetries can be found in the coplanar waveguide (CPW) and in the slot-line with their dominant mode distributions. In many case, the symmetry planes are generated due to specific source con

26 (a) - ------ ---- --. IIgi ------ ----—..PMC wall (b) ----------------- - I 91 --- -- PEC wall Figure 2.8: Reduced equivalent circuits for symmetric 2-port network depending on excitation configurations: (a) original network, (b) reduced equivalent network for even mode excitation, and (c) reduced equivalent network for odd mode excitation.

27 figurations. For example, a symmetric 2-port network with even and odd mode excitations 10 can be transformed into either of two simplified equivalent networks, depending on source arrangements, as shown in Fig. 2.8. In this case, to characterize the original network, only half of the problem can be solved with PEC or PMC condition according to a given source configuration. Obviously, with the CG type linear equation solver the reduced equivalent problem requires much less computational costs, since the BiCG and COCG routines used in this study demand 0(nmN) floating point operations for an N x N matrix ' and O((nmN) for N x TN matrix, where n < n for NT < N. For example, if an N x N matrix takes n = N iterations for convergence and an 1N x N (N = N/2) matrix needs n = N/2 iterations, then O((nmN) = 0(nmN)/4. Hence, the half-reduced equivalent problem takes only 1/4 of the computational cost of the initial problem. For structures having two symmetry planes, when the symmetry is fully accounted, the required operation counts for the equivalent problem are reduced to 1/16 of the operation counts for the original problem. As mentioned in the previous section, the PEC and PMC boundary conditions are implemented without any further theoretical and practical complications. On the contrary, those will provide great advantages on reducing the problem size to a half, a quarter, or even to a 1/8 of the original problem. When the number of FLoating point OPerations (FLOPs) is considered, as summarized in Table 2.1, it becomes clear what is the real advantage gained due to the reduction in problem size. This is one of the most important aspects of the symmetry boundary conditions and the reason why the symmetry of a given problem should be fully exploited whenever possible. 10See appendix C. "l See section 2.5.

28 Symmetry Planes Problem Size FLOP Count 0 mN,0(mN2) 1 mN/2 0(mN 2)/4 2 mN/4 ((mN2)/16 3 mN/8 0(mN2)/64 -k mN/(2k) 0(mN2)/(22k) Table 2.1: Summary of FLoating point OPeration (FLOP) counts for reduced equivalent problems. In the table, the number of iterations for convergence are assumed to be the maximum which is the matrix size. 2.7 The Question of Spurious Solutions: A Rationale for Using the Edge Basis Function 2.7.1 The Origin of Spurious Solutions The occurrence of physically incorrect spurious solutions in the various numerical techniques, such as the FEM, the method of moment (MoM), etc, has been noted, and many attempts have been made to explain and cure spurious modes [40], [43]-[17], [45], [46], [59]-[71]12. It has been found that these spurious solutions satisfy the finite element equations, but not the physical boundary conditions, and occur regardless what combination of vector components are involved. Thus, the spurious modes and the physically correct solutions may be discriminated by applying the exact boundary conditions on the inter-element and inter-material interfaces. On the interface of two different materials or on the boundaries of interior sub-domain elements, the tangential electric (ni x E) and magnetic (nf x H) fields should be continuous if there are no current sources, while the normal components ni. E and n - H are allowed to be discontinuous across the material interfaces. It has also been observed that the spurious solutions exhibit non-zero divergence fields even without containing any sources. Hence, the spurious solutions can be distinguished a posteriori by monitoring the divergence of the computed solutions or the behavior of eigenvectors. 12From the abundance of literature dealing with the spurious modes phenomenon, one can imagine how difficult and how important this task is.

29 This observation has led to the speculation that the enforcement of the divergencefree condition through the penalty function method would cure the occurrence of the non-physical solutions [8]-[12]. However, it has been proved [43, 46, 72] that the enforcement of the divergence-free condition cannot insure spurious-free solutions. Furthermore, even though this penalty function approach has gained partial success by removing or relocating the spurious solutions from the spectral region of interests, it is considered to be "sweeping dust under the rug" [16]. Recently, a rigorous theoretical analysis has been presented [43] on the cause of the spurious modes in numerical electromagnetic field eigenvalue problems. In view of the vastness of the spurious modes phenomenon across many different numerical methods, a source of the phenomenon has been identified as a misconception of the discretization of a parameter dependent indefinite or semi-definite linear operator. When such an operator family is applied to an n-dimensional subspace of its domain, the family image in general may contain a subspace of dimension higher than n. As a result the discretization of the parameter independent set of expansion and weighting functions of equal cardinality has the potential to fail. This reasoning is not incompatible with the conjecture in [73] which describes the root cause of the spurious solutions as the improper approximation of the null-space of the curl operator. 2.7.2 Why the Whitney Edge Element? Based upon the above observations, the Whitney edge element [13]-[18], [74] is employed in this study for the implementation of a robust finite element method known a priori to produce spurious-free solutions. As is already discussed in section 2.4, the edge basis function has many desirable properties, including its divergence

30 free condition 6 V-W =0 = V-E=V.(ZxiW )=0 j=l and accurate modeling of boundary conditions for tangential as well as normal electric field components. Furthermore, the discrete algebraic-geometric-differential [75]-[77] structure of the Whitney element closely matches the continuous one, providing a good approximation method. The Whitney complex can be represented as +gradl w' 2 2 I W3 (ascending sequence) (2.31) W gads I f~ W2 - W3 (descending sequence) (2.32) where WP is the finite dimensional subspace generated by taking linear combinations, with complex coefficients, of p-Whitneys. The adjoint operators grad*, curl*, and div* coincide with -div, curl, and -grad, respectively, on its complement domain and the Whitney elements of order p = 0,1,2,3 correspond to node, edge, facet, and volume elements, respectively. Considering the inherent mathematical symmetry of structure in electromagnetics -grad - div - -s E -4 B — 0 (2.33) curl -grad div Q - D -- pe, (2.34) where 1 and Q denote electric and magnetic scalar potentials, respectively, the remarkable similarity of the Whitney complex to Maxwell's equations leads to its excellent modeling capability for electromagnetic fields. Further, the conformity of the Whitney complex W~ C domain(grad) (2.35) W1 C domain(curl) (2.36) W2 C domain(div) (2.37)

31 reinforces the similarity between them and gives some of the basic properties of the Whitney elements: * Function W~ (node element) is continuous across facets * The tangential component of W1 (edge element) is continuous across facets * The normal component of W2 (facet element) is continuous across facets. In this study, the Whitney 1-form, W1, the so-called edge element, is employed to model the electric field distributions in each sub-domain. As mentioned above, WI1 shares many basic properties with the electric field, E, and thus renders spurious-free solutions. 2.7.3 Eigenvalue Analysis with Edge Element From the source-free vector wave equation (2.7), the system of linear equations (2.20) obtained after applying Galerkin's procedure is reduced to the generalized eigenvalue problem. The generalized eigenvalue problem can be stated as: ([U] - A[V])[x]= 0 (2.38) where M M [U] = U [V [] U [ e=l e=l and the matrix entries of [U6] and [V6] are defined as in equations (2.21) and (2.22). The eigenvalue A is the square of a propagation constant, ko = v/A, and has real non-negative values for a lossless isotropic medium. A set of eigenvalue and eigenvector, {A, x}, represents the electromagnetic field distribution at a specific frequency, fo = (c/27r)ko. The spurious eigenvalues are identified by examining the corresponding eigenvectors. In general, the eigenvectors of the spurious modes exhibit nonphysical field distributions, such as node-to-node variations, or have significantly

32 higher divergence values: jII IV E 2 d >. With the edge element, it has been found that no spurious modes are interspersed among the physical solutions. However, the computed set of eigenvalues of the generalized eigenvalue equation (2.38) contains multiple zero eigenvalues (A = 0) which correspond to the DC gradient electric fields, 4). Note that any gradient fields satisfy the vector wave equation with V x (V4)) = 0. These zero eigenvalues can be interpreted as a separate class of spurious modes since their existence does not interfere with finding physical solutions, except when the discrete numerical scheme does not approximate them properly. This inaccurate approximation has been considered as a possible source of the non-physical spurious modes [66, 73]. The number of zero eigenvalues, N,,o,, can be accurately predicted for structured tetrahedral mesh, in which the domain is divided into bricks and then each brick is subdivided into 5 tetrahedrons, as follows: (Nx - 2)(Ny - 2)(N - 2), for E-field formulation Nzero = (2.39) N:NyNz - 1, for H-field formulation where N1, Ny and Nz are the number of nodes in x, y, and z directions, respectively. For the E-field formulation, the Nzero is coincident with the number of interior nodes in a simply connected domain surrounded with the PEC boundary conditions. For the H-field formulation, the number of zero eigenvalues is identical to the number of edge-tree branches in the domain with PEC outer surface, which will be the natural boundary condition. The number of edge-tree branches can be computed from the difference between the total number of edges and the number of independent algebraic equations [78]. Tables 2.2 and 2.3 show examples of the number of zero eigenvalues

33 (Nr N, N ) Total # of # of surface Matrix Size # of interior AN ero edges edges M nodes (3,33) 90 72 18 1 1 (4,4,4) 252 162 90 8 8 (5,5,5) 540 288 252 27 27 (6,6,6) 990 450 540 64 64 (7,7,7) 1638 648 990 125 125 Table 2.2: Number of zero eigenvalues for simply connected cavity with E-field formulation. (Nt, Ny, Nz) Total # of # of surface Matrix Size Total # of Nzero edges edges M nodes (2,2,2) 18 8 7 (3,3,3) 90 0 90 27 26 (4,4,4) 252 0 252 64 63 (5,5,5) 540 0 540 125 124 (6,6,6) 990 0 990 216 215 (7,7,7) 1638 0 1638 343 342 Table 2.3: Number tion. of zero eigenvalues for simply connected cavity with H-field formula and other relevant quantities for the E- and H-field formulations in various mesh configurations. As shown in the tables, the E-field formulation requires a lower computational cost, i.e., a smaller matrix size M, than that of the H-field, due to the PEC boundary condition, as described in section 2.6.1. 2.8 Excitation Mechanism In the application of the FEM for the characterization of high frequency and high-speed interconnects, many different types of excitation mechanisms could be employed. Since the extracted circuit parameters should be independent of the excitations, one can choose a convenient excitation method suited to specific applications. For faster convergence, two dimensional (2D) field distributions on each face of the port can be used: however, this approach requires another 2D electromagnetic field simulator. In this study, ideal impulse current sources are employed and imposed on

34 the corresponding edges of the tetrahedrons. Even though the ideal current sources generate impulses at the sources, the electromagnetic fields quickly conform to the appropriate field distribution according to the boundary conditions, even in the close vicinity of the sources. In the finite element formulation, the ideal current source could be implemented either in conjunction with the inhomogeneous wave equation, or using the appropriate boundary conditions with the homogeneous wave equation (see Eq. (2.7)). Consider the right-hand side of the weak form of the homogeneous wave equation (2.13): b = -jkoZo j W. (n x H) dS. (2.40) In the kernel of the integral, the n x H can be interpreted as a current source J and the dot product between J and We can be considered as a line current Jl on the edge of the surface triangle aQe. As a result, the excitation vector [be] with an appropriate scaling factor is given as X W. J dS = i-j for a source at jth edge. (2.41) Jabe When there are several sources, as in the case of a multi-port network, the excitation vector [be] has non-zero elements and the positions of those non-zeros correspond to the global indices of the excitation edges, as shown below: [be]T = [0 0... 0 bi 0. 0 bj 0... 0] (2.42) where the bp (1 < p < N,p = {i,j}) are an appropriate non-zero pair of numbers for various sets of excitation mechanisms such as the even, odd, or sum excitation discussed in Appendix C. In Figure 2.9, two different configurations of the current source excitations used in this study are illustrated. As shown in the figure, horizontal or vertical excitations produce identical results.

35 (a) (b) Figure 2.9: Two types of ideal excitation configurations: (a) horizontal current source on the transmission line plane. (b) vertical current source. 2.9 Customized Generation of Tetrahedral Mesh The first step of the finite element analysis is a discretization of the domain of interest into appropriate sub-domain elements. In this study, for the discretization of a given three-dimensional space Q into a number of tetrahedral sub-domain elements Qe, a simple and efficient structural mesh generation algorithm is developed. In particular, to accommodate and model effectively the large contrast in geometrical details commonly found in high-speed high-frequency interconnects, both uniform and non-uniform meshing schemes are implemented. In non-uniform mesh, variable size rectangular bricks are created first and then each brick is subdivided into five tetrahedrons in two different ways as shown in Figure 2.10 (see Table 2.4 for labeling system). Note that those two different ways of subdividing a brick are mirror images of each other. For continuity of the edges across the brick faces, all of the neighboring bricks sharing a face with each other are assigned to have different types of dispositions. One advantage of using the structural mesh is its regularity, which allows accurate prediction of the total number of edges (i.e., the total number of unknowns) and tetrahedrons, the edge connectivities, and the resulting memory requirement for a given problem. For an edge which is parallel to the Cartesian axis, there are 4

36 (a) (b) Figure 2.10: Two types of rectangular bricks which containing five tetrahedral elements: (a) Type I, (b) Type II Tetrahedron Type I Type II {1,2,4, 6} {1,2, 3, 5} {1, 5, 6, 7} {2, 3, 4, 8} {ij,k,l} {1, 3, 4, 7} {3, 5, 7, 8} {4, 6, 7, 8} {2, 5, 6, 8} _{1,4,6, 7} {2,3,5,8} Table 2.4: Two types of tetrahedral labeling systems.

37 tetrahedrons surrounding the edge and 13 edges are associated with that. However. for a typical diagonal edge, 6 tetrahedrons are connected to that edge and therefore 19 edges are associated with that. As a result, the total number of non-zero entries in a row or column of an arbitrary FEM matrix is pre-determined to be 13 or 19 with the structural meshing scheme. Further, the average edge connectivity would be determined to be 15 from the fact that there are 12 edges parallel to the Cartesian axis and 6 edges diagonal in a rectangular brick 13. The pre-determined edge connectivity allows efficient memory management and planning and building effective data structures for the best possible use of computing facilities, since the total number of edges in Q determines the overall problem size to be solved. The relationship between the number of tetrahedrons, nodes, and edges for a simply connected region, whether the structural meshing scheme is used or not, can be estimated as follows [79] Nedge = Ntet + Nnode + Nsurf - 3 (2.43) where Ntet, Nedge, Nnode and NSurf are the total number of tetrahedrons, edges, nodes, and surface nodes, respectively. This formula is very useful for the prediction of a problem size in either case, whether structured or unstructured mesh is used. However, for an unstructured mesh, it is not always possible to predict accurately the number of tetrahedrons, nodes and surface nodes, and even further, the edge connectivity varies over a large range. Roughly, the total number of edges can be estimated as six times the total number of nodes and as a result the estimated number of total non-zero entries in a FEM matrix is 90 times of the total number of nodes. The average matrix density also can be estimated from the average number of non-zero entries in a row (or column) 13Average edge connectivity = (12 edges x 13 + 6 edges x 19)/18 edges = 15

38 Total # of nodes: Anode ANxAy N Total # of tetrahedrons: Ntet 5(N. - 1)(N - 1)(AN - 1) Total # of surface nodes: 2(NxNy + NyAN + AT A) Nsurf -4(/' + Ny + N,) + 8 Total # of edges: Nedge 6N:NyN, - 3(N:Ny + NyN, + NNrx) (# of unknowns) +5(Nx ~ Ny + Ny + Nz) + 2 # of non-zero entries 13 or 19 in a row or column Average # of non-zero 15 entries in a row or column Total # of non-zero 90NxNyNz -75(NxNy + NyNz + NZNx) entries in FEM matrix: Mo -75(NA + Ny + Nz) + 30 Average matrix density 2.5/(NxNyNz) Table 2.5: Estimations of some of the important quantities in the structural tetrahedral mesh. It is assumed that the domain is simply connected and there are no interior PEC surfaces which will reduce the total edge counts. N~, Ny and Nz are the number of nodes in x, y and z direction, respectively. and the size of the square matrix. Since the average edge connectivity has been determined to be 15 for any arbitrary problem, the overall problem size is directly proportional to the total edge count. Therefore, the sparsity of the FEM matrix can be estimated as given below sparsity = (1 - 2 ) x 100 [%] (2.44) Nedge 2.5 N N x 100 [%]. (2.45) Note that the sparsity of the FEM matrix is approaching 100 % as the size of the problem, Nedge,, is becoming larger. These relations are summarized in Table 2.5. It has been found numerically that the iterative linear equation solvers mentioned in section 2.5 converge much faster with the uniform mesh than the non-uniform scheme, implying that the uniform mesh leads to a better matrix condition number. However, even though uniform mesh has the definite advantage of fast solution time for a fixed number of unknowns, it is not always advantageous to use the uniform mesh when there are large variations in detailed geometries as found in usual high

39 speed high-frequency interconnects. In general, the uniform mesh generates a much larger matrix system than the non-uniform case when there are very small details in a structure and eventually causes a longer solution time. With these observations, the non-uniform structural meshing scheme is developed in this study to provide effective mesh data for circuits with large variations in detailed dimensions. Even in the non-uniform meshing scheme, the quality of each sub-domain element is also very important, since one badly distorted tetrahedral element could leads to a poor matrix condition number. 2.10 Treatment of Open Boundary Problems Accurate modeling of high-speed high-frequency interconnects often requires treatment of closed as well as open boundary problems. However, due to the fact that one of the fundamental ideas based on partial differential equation (PDE) techniques is the restriction of the computational domain to a finite region, the open domain problem requires special techniques for appropriate approximation of the outer boundaries, as shown in Figure 2.11. In general, the difficulties in the treatment of open domain problems with appropriate outer termination boundaries have become a hindrance in the efficient PDE solution of practical problems. In the context of the finite element method, the open domain problems can be simulated in several different ways, such as the FEM-boundary integral (FEM-BI) type method [80]-[87], the application of higher-order absorbing boundary conditions (ABCs) [38, 39], [88][117] or the use of the perfectly matched layer (PML) [118]-[131]. The FEM-BI type method is considered more rigorous than the other methods since it applies exact boundary conditions through Green's function. However, the FEM-BI type method produces a fully populated matrix subsystem, which causes difficulties in storing and

40 solving large systems of equations. Compared to the global boundary condition of the FEM-BI method, the ABC and PML type approaches produce local boundary conditions. The locality of the boundary conditions contributes to the generation of a sparse matrix system and offers great advantages in storing and solving the matrix system. In higher order ABCs, accurate prediction of the local curvature of the wave front is crucial to the overall accuracy as well as convergence of numerical methods. Unfortunately, however, the curvature of the wave front is, in general, unknown and therefore the accuracy of ABCs is limited by its poor approximation of the tangential electromagnetic fields at the outer surface. Also, ABCs reveal difficulties in the termination of the inhomogeneous boundary value problems since they were developed for homogeneous open spaces with plane wave incidence. Moreover, these ABCs do not provide optimum design procedures for a given error tolerance except by increasing the separation distance between the object and the boundary surface, which will eventually intensify the computational burden and require the use of higher order boundary conditions. As a results, the ABC type approaches lead to inferior numerical performance such as slow convergence and spurious reflections at the outer boundary. With these observations, the design of absorbing material with the perfectly matched layer (PML) concept becomes an alternate way of terminating outer boundaries and thus simulating infinite space. The PML technique is based on two essential ideas, namely impedance matching at the interface and damping of the waves in the material for least amount of reflections coming from the finite thickness of the material. These types of absorbing layers, in general, produce a symmetric FEM matrix and provide flexible design parameters in order to satisfy a given error bound. Moreover, the PML type absorber can be placed within fractions of a wavelength

41 Open Space Es HSfields (a) Closed/Finite Space Es HSfields Boundary (b) Figure 2.11: Treatment of open boundary problem: (a) Original open boundary problem, (b) Equivalent closed/finite boundary problem.

42 without sacrificing accuracy. Among the known artificial absorbers, Berenger's PML has demonstrated excellent performance for 2D and 3D problems in the context of a finite difference time domain technique. However, Berenger's PML equations are non-Maxwellian and suffer from instabilities when applied to layered inhomogeneous problems. Furthermore, the implementation of the absorber in the context of FEM is not straightforward [131]. Another similar type of PML has been presented recently, introducing uniaxial anisotropic materials to implement a non-reflecting interface for an arbitrary angle of incidence. Such an interface condition with a material having broadband characteristics, satisfies Maxwell's equations and has been implemented using anisotropic material. For example, the characteristics of the anisotropic absorbing material facing z = constant plane are diagonal tensors: p 0 0 []PML = 0 p 0 c, (2.46) 0 0 Il/p p0 0 []PML = 0 p 0 P (2.47) o o i/p where p = a - /3 and the values of a and (3 have to be carefully chosen for best numerical performance. In theory, the above PML material yields a reflectionless interface for an arbitrary angle of incidence, polarization, and frequency. However, the performance of the PML concept in the FEM context reveals several practical difficulties, including the difficulty of choosing appropriate a and / for various structures, slow convergence in the iterative matrix solver, and direction dependency of the material.

43 In contrast to the above anisotropic absorber, an isotropic artificial absorber has also been devised to simulate a free space buffer layer. This type of material is not independent of the angle of incidence. However, its loss mechanism is clearly understood and thus provides a well defined design procedure for optimum performance. Also, it can be designed to be a PML absorber for well defined field distributions as inside a waveguide. Furthermore, the absorber can be easily incorporated into an existing FEM code and renders much faster convergence compared to the anisotropic one. In view of the advantages of using the isotropic absorber, in the following subsections, design procedures are presented the absorber, for MMIC and waveguide applications. 2.10.1 Isotropic Absorbers for MMIC Applications Transmission line structures, as commonly encountered in usual MIC/MMIC applications such as microstrip, CPW and slot line, support quasi-TEM or hybrid type modes, as opposed to the TEM nature of plane waves excited in free space. Moreover, when there are circuit discontinuities, evanescent or higher order modes are excited which leak power into the substrate. As a result, design of an appropriate layer capable of absorbing all of these waves, including the dominant mode, is a rather challenging task. In our application of the absorber to MMICs, we find three different situations. The first is encountered when the absorber is used as a termination boundary in transmission lines to simulate a matched load, the second occurs when we need to simulate free space, and the third corresponds to the termination of a substrate. For the design of effective absorbing layers for the above three cases, matching the wave

44 impedances at the absorber interface will result in a non-reflecting surface: Zr = ZAA where Zr and ZAA are the wave impedances in the dielectric and the artificial absorber, respectively. Since the wave impedance for a TEM wave is identical to the intrinsic impedance of the medium in which it propagates, matching the intrinsic impedances across the material-absorber interface is equivalent to matching the wave impedances. Similarly, in a homogeneous waveguide terminated with an artificial absorber, TE and TM mode wave impedances in both media are identical, if the intrinsic impedances are the same. Therefore for effective absorption and minimal reflection of the surface waves and higher order evanescent modes created by a circuit discontinuity, the intrinsic impedances at the absorber interface need to be matched. In addition, to avoid any secondary reflections at the end of the absorbing layer, the material parameters need to be optimized to maximize the attenuation of the waves transmitted into the absorber. This leads to complex values for the electric permittivity (EAA = CAA - JCAA) and magnetic permeability (AA LAA - JAA UAA) of the absorber with large imaginary parts. Based on the above observations, the design procedure of the artificial absorbing layer interfacing a dielectric layer with E:r and jur can be summarized as follows: 1. Choose the operating frequency (fO = c/A\) and an arbitrary thickness of the absorber (St). 2. Compute relative intrinsic impedance of the dielectric layer (i7r = a/r/Cr). 3. Choose Re(EAA) = Cr and Im(CAA) such that Im(CAA) > -~ ()O) l ln(a), where

45 a is the maximum magnitude of the reflected wave which can be tolerated at the interface. 4. Choose P/AA = r2CAA (impedance matching condition). It is noted that even though the absorbing layer is designed under the assumption of normal incidence, the design procedure can be modified for the elimination of a wave incident from a certain angle. 2.10.2 Isotropic Absorbers for Waveguide Applications In homogeneous waveguides where mode patterns are well defined, an isotropic absorber can be optimized to simulate an infinite waveguide. Since for each propagating mode at a certain frequency one can define an equivalent incident angle using ray optics, optimization of the absorber parameters to that particular angle of incidence can minimize the reflections at the interface. The thickness of the absorbing layer is chosen to provide enough attenuation of the waves in the absorbing material and to minimize the computational domain. The following steps summarize the design procedure for the isotropic absorber for a given mode (TEmn or TMmn) and a specific guided frequency (fg). The a and b (a > b) are the waveguide dimensions and cr and Pir are the relative permittivity and permeability, respectively, in the waveguide. 1. Choose an arbitrary thickness of the absorber (St). 2. Compute relative intrinsic impedance of the waveguide material (r = Ar/r). Usually, Ag/10 is enough for St. 3. Let Re(CAA) = R(cr) and Im(eAA) such that Im(iAA) > ln(ca), where a is the maximum magnitude of the reflected wave which can be tolerated at

46 the interface. 14 4. Calculate IAA using the following equations, derived from the interface matching conditions for TE and TM wave incidence with angle i0 * for TEmn ({m,n 0 < m,n < oo and m = n ~ 0}) modes, 2AAr Cos2 0- 2er- ( * for TMmn ({m,n 1 < m,n < oo}) modes, PAA -= [e2 sin2 i + AA cos2 i] (2.49) Cr {AA where = cos-1( l - (f/f)2) fc = - 2 /()2 + ( )2 and co is the speed of light in vacuum. It is noted that even though the absorbing layer is designed for a specific mode at a given frequency, and that as a result it has narrow bandwidth characteristics, its performance under single mode operation can be maximized by adaptively changing the absorber parameters at each frequency. As a result, the application of the absorber to the frequency domain FEM can be easily optimized. For the evanescent modes below the cutoff frequency, the performance of the absorber in the context of the FEM cannot be easily evaluated. However, in practice, any instability effect has been observed in the evanescent mode region with the absorber. It should also be pointed out that the performance of the isotropic absorber near the cutoff frequency is degraded, as happens with other types of PML. 2.11 Passive Lumped Elements In high-frequency high-speed circuits, the passive and active lumped elements constitute indispensable ingredients along with distributed circuits, while the sizes 14 Even though one can choose a much larger value of Im(CAA) than the criterion, it is not always advisable since it will require much higher discretization to model the highly varying fields inside the absorber.

47 of the lumped elements are negligible compared to that of the distributed network. In spite of their small size, in many cases, the lumped elements strongly interact with the surrounding environment, including the distributed circuits. As a result, it is very important to have a three-dimensional full-wave electromagnetic field characterization of the circuits having both distributed and lumped elements, and several investigations have recently been reported with time domain PDE techniques [132][143]. The incorporation of the lumped elements into the finite element method, however, has not been well understood, nor its effectiveness and possibilities assessed despite its importance. In this study, the passive lumped elements, such as resistors, inductors, and capacitors, are modeled in three-dimensional finite element techniques in two different manners. The effect of the lumped elements can be implemented starting from the inhomogeneous wave equation: V X ur71(V X E) - k2er E = -jwJZ. (2.50) Similar to the method given in section 2.3, Galerkin's procedure is applied to Eq. (2.50) on Hilbert space and the weak form is found as: for j = 1,...,n, ( x E,V x Qj) = -jkoZo(Qj,Ji) +jkoo (H x Qj) nidS (2.51) where the weighting functions, Q = {Q1, Q2,..., Qn} C fl2, are defined on a Hilbert space and eventually the electric field, E, will be also expanded on a Hilbert space, 7-l2, with the same weighting function. Now, the lumped elements could be implemented through the two different methods discussed below, using the right hand side terms of Eq. (2.51). In the following two subsections, the detailed mathematical formulations will be discussed for the so-called volume current approach and surface impedance approach.

48 2.11.1 Volume Current Method In this method, the first term in the right hand side of the Eq. (2.51) is utilized to incorporate lumped impedance elements: = -jk0Z0o(Qj Ji) for j=1,...,n. (2.52) where ZO is the free space intrinsic impedance. In the above equation, the volume current density can be represented in terms of electric field and conductivity as: sI 1I I ZL J, = vE (2.53) where oa = l/(ZLS) and ZL, s, and / represent, as shown in Figure 2.12, load impedance in ohms (Q), cross-sectional area, and length of the load, respectively. After Eq. (2.53) is combined into Eq. Figure 2.12: Impedance load. (2.52), we obtain 1 =1 = -jkoZo Zs (Q3' E) ZLS for j=l,...,n. (2.54) Now, the electric field is expanded in Hilbert space based on tetrahedral elements as 6 E = ExiQ. i=1 and as a result we have 1 6 Z1 = -jkoZo z X (Qj, Qi) for j=l,...,n. ZLS Z=1 i=1 (2.55) (2.56)

49 If the impedance load is located on the ith edge in a tetrahedron. the inner product between the basis and weighting functions becomes constant: (Qj, Q) = AI 'Qj, where AV is the volume of the load given as AV = sl. Therefore, in most practical situations, Eq. (2.56) is simplified to 12 Z1 = -jkoZo-x,. (2.57) ZL This term only contributes to the diagonal terms, A.-, in the FEM matrix [A] in Eq. (2.20). 2.11.2 Surface Impedance Method The incorporation of the passive lumped elements into the finite element method can be also implemented using the surface integral term, the second term of the Eq. (2.51), 12 = -koZo (H x Qj) nidS for j=1,...,n. (2.58) and the resistive sheet condition [144]: n x (n x E) =-R n x (H+- H-) (2.59) where H+ and H- are magnetic fields defined on the surface s+ and s-, respectively, as shown in Fig. 2.13. For a thin dielectric layer of thickness t, the equivalent resistive sheet having surface resistivity R (Q/O) can be derived from the equivalent polarization current [145, 146]: R = zo (2.60) jko(cT - l)t' For a thin dielectric layer, 12 can be split into two terms for the upper and lower surfaces: 2 = -koZo Qj (n x H+) dS + I Qj (n x H-) dS (2.61) +s s

50 n rS+ er Q n R Figure 2.13: Resistive sheet and equivalent thin dielectric layer. and on the assumption that kot < 1, it is further reduced to 2 = -jkoZo j Qj [n x (H+ - H-)] dS. (2.62) Now, Eq. (2.59) is substituted into the above, yielding 2 = -koZo (n x Qj) (n x E) dS (2.63) and its discretized form on Hilbert space with the edge basis function is found as: 2 = -JkoZo xi (n x Qj). (n x Qi) dS for j=l,...,n. (2.64) This term is evaluated numerically for a given surface and its resistivity. To a first order, this equation can be used to simulate the presence of a thin dielectric layer by choosing an equivalent resistive sheet having resistivity R. Alternatively, a resistive sheet may be equivalently replaced by a thin dielectric layer

51 having thickness t and a relative permittivity of, 1 - 3kotR (2.65) Note that with the latter approach, for resistive lumped elements, a small volume of dielectric material is allocated to have the equivalent dielectric constant. This approach can be also expanded to capacitive and reactive elements by substituting the ZL for load impedance for the resistivity R in Eq. (2.65). In this study, both approaches are implemented, and the performance and accuracy are compared.

CHAPTER III VALIDATION AND PARALLELIZATION 3.1 Introduction Presented in this chapter is validation of the edge-based finite element method by means of the generalized eigenvalue problems for lossless cavities. The motivation behind of using the eigenvalue problem to validate the FEM is the fact that the excitation and the eigenvalue problems share the identical matrices, and the exact solutions can be found for certain eigenvalue problems. As a result, by examining the accuracy of the eigenvalues for various problems, the developed edge-based FEM code can be validated. Furthermore, the validated FEM code is parallelized on the distributed memory parallel computer (the IBM SP2) using the message passing paradigm (MPL). Two types of parallelization schemes are examined in this study and the performance of each scheme is investigated in depth. Having done the validation and parallelization of the FEM, a design guideline for artificial absorbing materials is introduced for waveguide and MMIC applications, and its effectiveness is verified with several examples. The use of absorbing materials is crucial for simulation of open boundary problems as proved whole throughout the thesis. 52

53 3.2 Eigenvalue Problems Seeking the solutions of eigenvalue problem using source-free Maxwell's equations allows basic understanding of the electromagnetic characteristics of a given geometry. In particular, the eigenvalues and the corresponding eigenvectors provide complete signature of the fundamental electromagnetic properties of the structure. Any field distributions in a cavity or guiding structure can be represented as a linear combination of eigenvectors. Strictly speaking, any closed non-radiating geometries support discrete eigenvalues and eigenvectors, in contrast to continuous counterparts of the open or radiating problems. For example, the eigenvalues and the eigenvectors of a lossless cavity, whether it is homogeneous or not, represent resonant frequencies and their mode distributions. On the other hand, the eigenvalues of lossy cavity become complex numbers and thus correspond to complex modes having attenuating as well as non-attenuating components. When there are no excitation mechanisms in domain Q, the generalized eigenvalue equation can be derived from Eq. (2.20) as: [U][x] = A[V][x] (3.1) where M M [U] U[Ue] and [V] = U[V e=l e=l and the matrix entries of the submatrices [Ue] and [Ve] are determined as given in Eq. (2.21) and Eq. (2.22): Ui = (r7VxWVxW,) 3f = 2W.W)

54 Here, the eigenvalue A is given as kA which is a square of the propagation constant. Note that for lossless isotropic material the matrix [U] is symmetric, and [V] is symmetric and positive definite since the basis functions Wi are independent. Using the fact that any positive definite matrix can be decomposed into two triangular matrices, [V] can be written as [V] = [C][C]T (3.2) where [C] is the triangular matrix. Now, the original eigenvalue equation is transformed into [U][x] = A[] (3.3) where [U] = [C]-1[U][C]-T and [k] = [C]T[x]. Therefore, for lossless isotropic medium, all the eigenvalues A become real numbers, since the matrix [U] is symmetric. The solution of the generalized eigenvalue equation could be found using the iterative method, such as Lanczos method, to take advantage of the sparsity of the FEM matrix. In this study, however, all of the eigenvalues including zeros are found using the standard factorization method. 1 Herein, the eigenvalues are calculated for several different types of homogeneous as well as inhomogeneous cavities, and the numerical results are compared with the analytic solution, if available, to validate the developed edge-based finite element method. When analytic solution is not obtainable for inhomogeneous problem, the excitation problem, having exactly the same configuration except the source, is solved to compare the resonant frequency. In many practical applications, the solution of the excitation problem can be found more accurately and efficiently than that of the eigenvalue problem with aid of efficient linear equation solver. Furthermore, 1There are standard EISPACK routines, such as RSG and QZHES, for generalized eigenvalue problems.

55 to characterize unknown dielectric material having complex dielectric constant, the excitation problems are solved and compared with the eigenvalues to extract real and imaginary parts of the permittivity. 3.2.1 Lossless Cavities (a) Empty Rectangular Cavity As a first step toward validating the developed edge-based FEM code, the generalized eigenvalue problem is solved for rectangular cavity shown in Fig. 3.1. The resonant frequency, fonp, of each mode in the cavity can be also obtained as fmnp =_ 1 (rn 72 + (_1)2 + (pr)2 27rVV6 a y b c where {m,n,p I m > 0 n > Op > O,m = n 7& 0} for TEZ and {m,n,p I m > 1,n >_ 1,p > 0} for TMZ mode. As is summarized in Table 3.1 and Table 3.2, the FEM-computed eigenvalues are accurate within ~1.0 % error bound for all the cases listed in the tables. Note that the degenerated modes having the same resonant frequency are also predicted very accurately in both cases. Furthermore, it is also observed that the total number of zero eigenvalues is equal to the number of internal FEM nodes as discussed in Section 2.7.3, and more importantly, no spurious modes are detected in the eigenvalue spectrum. With the remarkable accuracy obtained by the edge-based FEM for the empty cavities, eigenvalues for inhomogeneous cavities having dielectric material are computed in the below. (b) Partially-Filled Cavity Having been proved the accuracy of the edge-based FEM for the simplest cavity structures, the accuracy of the FEM is tested for more complicated geometries. In this study, the eigenvalues of the partially-filled cavity, as shown in Fig. 3.2, with

56 I/ I I s — -r c ---- a ----- I / Figure 3.1: Rectangular cavity having PEC walls and dimensions a, z-direction, respectively. b, and c in x-, y-, and mode (m,n,p) Analytic FEM error [%] (1,1,0) 21.2133 21.2090 0.02 (1,0,1) 21.2133 21.2090 0.02 (0,1,1) 21.2133 21.2090 0.02 (1,1,1) 25.9808 26.1173 -0.53 (2,1,0) 33.5409 33.3714 0.51 (0,2,1) 33.5409 33.3714 0.51 (1,0,2) 33.5409 33.3714 0.51 (1,1,2) 36.7124 36.7696 -0.07 (2,1,1) 36.7424 36.7696 -0.07 (1,2,1) 36.7424 36.7696 -0.07 (2,2,0) 42.4266 42.2972 0.30 (2,0,2) 42.4266 42.2972 0.30 (0,2,2) 42.4266 42.2972 0.30 Table 3.1: First few eigenvalues of an empty rectangular cavity having a = b = c = 10 mm. fo [GHz] mode (m,n,p) Analytic FEM error [%] (1,1,0) 6.2500 6.2543 -0.07 (1,0,1) 8.3852 8.3766 0.10 (2,1,0) 9.0141 9.0040 0.11 (0,1,1) 9.0141 9.0336 -0.21 (1,1,1) 9.7627 9.7799 -0.18 (2,0,1) 10.6064 10.5696 0.35 (1,2,0) 10.6799 10.7353 -0.52 (2,1,1) 11.7261 11.7313 -0.04 Table 3.2: First few eigenvalues of an empty rectangular cavity having a = 40, b = 30, and c = 20 mm.

57 J iC':- - --- -- ----- - -- - - - --, / a --- Figure 3.2: Partially-filled rectangular cavity having dielectric material with er and thickness c'. two different dielectric material are computed and the results are compared with the analytic data. As tabulated in Table 3.3 and Table 3.4, the first few of the FEM-computed eigenvalues for Cr = 2.5 and 12.9 are listed and revealed excellent agreement with the analytic data. It is observed that the error for the higher dielectric constant material shows slightly higher than that of the smaller one mainly due to the numerical discretization error. Similar to the previous empty cavity examples, there are no spurious modes in the whole frequency spectrum. (c) Loaded Rectangular Cavity As a last example of the lossless eigenvalue problem, dielectric cube having dimension a', b', and c' with dielectric constant Cr is placed at the center of the cavity bottom plane as shown in Fig. 3.3. In this section, to quantify the effect of the size of the dielectric cube and its dielectric constant on the resonant frequency of the dominant mode, the dielectric constant Cr is varied from 1 to 50, and the filling

58 fo [GHz] mode (m,n,p) Analytic FEM error [%] LSM111 5.6736 5.6469 0.47 LSElol 7.8617 7.8578 0.05 LSM211 7.9493 7.8826 0.84 LSEo11 8.4350 8.4087 0.31 LSM121 8.8660 8.8654 0.01 LSE11 9.1125 9.1031 0.10 LSM121 9.2044 9.1104 1.02 LSE201 9.8677 9.8513 0.17 LSM311 10.3579 10.2504 1.04 LSM311 10.4876 10.3985 0.85 Table 3.3: First few eigenvalues of a partially-filled rectangular cavity having a = 40, b = 30, c = 20, c' = 5 mm, and ~r = 2.5. fo [GHz] mode (m,n,p) Analytic FEM error [%] LSM111 3.9809 3.9302 1.27 LSM211 4.5625 4.5074 1.20 LSElol 4.6704 4.6360 0.74 LSM121 4.8837 4.8416 0.86 LSEoll 4.8977 4.8418 1.14 LSElll 5.1405 5.1169 0.46 LSM311 5.2048 5.1684 0.70 LSM221 5.2462 5.1988 0.90 LSE201 5.4039 5.3750 0.53 Table 3.4: First few eigenvalues of a partially-filled rectangular cavity having a = 40, b = 30, c = 20, c' = 5 mm, and e, = 12.9.

59 = '/a) = (/b) = (c'c). (3.4) // — a / ----^ wheFigure o is the3.3: Rectangular cavity having dielectric cube anant mothe fcenter of the empty cavity bottomwhile having dielectric constant loade and a' b an d c' in x-, y-, and z-dimension, respectively. factor r ha s been varied from 0.1 to 0.25, where r is de fined as: r = (a'/a) = (b'/b) = (c'/c). (3.4) In a ddition, the frequency shift due to the dielectric materia l is define d as: of as a function l [%]cearly revealed for several different size of the material. As can be observed, as the permittivity resona or filling factor r increases, the resonant frequency of the loaded cavity f. becomes lower and in turn sta increases. he responds tore, for a givenempty cavity, whie r = 1.0 t raightforully warfilled cavity.o predict the resonant frequency of the loaded cavity. In fact, for any unknown dielectric

60 20........ -o — r=0.25 15 - -- r= 0.20 --- r= O.15 ----- r= 0.10 0 0 10 20 30 40 50 Er Figure 3.4: Variation of the resonant frequency of the dominant mode for various size and Er of the dielectric material. material, it is possible to accurately determine the dielectric constant of the material by comparing the resonant frequencies of the empty and loaded cavities, and looking up the curves in Fig. 3.4. When the material has complex permittivity, the real part can be extracted by measuring the shift of the resonant frequency, while the imaginary part can be deduced from the quality factor of the loaded cavity as will be discussed in the next section. It should be pointed out that the accuracy of the resonant frequencies of the dielectric loaded cavity has been found within 2 % by comparing the data obtained from the excitation problem. 3.2.2 Lossy Cavities When a lossy dielectric material having complex permittivity Cf is placed inside of a cavity, the eigenvalues become complex numbers and the quality factor of the cavity becomes finite as illustrate in Fig. 3.5, in contrast to infinity for lossless case.

61 Figure 3.5: Typical resonance curve of a lossy cavity. The quality factor Q is defined based on the center frequency, fo, and two adjacent frequency points, fi and f2. In general, the quality factor (Q-factor) of a lossy cavity is defined as: Q fo (3.6) BW where the frequency Jo is the center frequency. As illustrated in Fig. 3.5, the bandwidth BW is defined as BW = f - fi, where fi and f2 are the frequency points correspond to 1/v/2 of the maximum. In this study, various lossy dielectric materials having different level of loss are modeled using the excitation problem and the intensity of the electric fields are measured to determine the resonant frequency and Q-factor. As can be observed in Table 3.5, the resonant frequencies of the lossy dielectric loaded cavities are mainly determined by the real part of the permittivity, while the imaginary part contributes to finite Q-factors. Note that the Q-factor is linearly proportional to the loss tangent, tan 8, of the material.

62 r r tanc fo [Gliz] BW' [MHz] Q 0.1 6.0 + J0.1 0.0167 6.2140 1.09 5700 0.25 2.25+ O0.0 0.0 6.0717 0.0 oo 0.25 2.25 + 0.001 0.00044 6.0717 0.222 27350 0.25 2.25 + 0.01 0.0044 6.0717 2.24 2710 0.25 2.25 + j0.1 0.044 6.0717 22.3 272 0.25 6.0 + jO.0 0.0 5.8044 0.0 o 0.25 6.0 + 0.01 0.00167 5.8044 0.93 6241 0.25 6.0 + jO.1 0.0167 5.8044 9.29 624.8 0.25 12.9 + j0.0 0.0 5.6032 0.0 oo 0.25 12.9 + j0.001 0.000078 5.6032 0.038 147454 0.25 12.9 + jO.01 0.00078 5.6032 0.38 14745 0.25 12.9 + 0.1 0.0078 5.6032 3.9 1455 Table 3.5: Resonant frequencies and Q-factors for lossy dielectric loaded cavities having various amount of losses. 3.2.3 Summary In this section, the validity and accuracy of the edge-based finite element method is proved using the generalized eigenvalue problems. In particular, the eigenvalues of the empty as well as dielectric loaded cavities are accurately preurdicted, and more importantly, no spurious modes are observed. In addition, it is found that the number of zero eigenvalues are equal to the number of internal tetrahedral nodes. For inverse problem of predicting the dielectric properties of a material, variety of dielectric configurations are examined. Finally, it is demonstrated that when there are lossy dielectric material in a cavity, the real and imaginary parts of the permittivity can be extracted from the resonance curves. 3.3 Parallelization 3.3.1 Motivations The three dimensional edge-based vector FEM is well established and applied for the characterization of many 2D and 3D circuits and interconnects which contain complicated dielectric and metallic configurations [33, 34, 28]. It is proved to be very

63 accurate, computationally inexpensive and highly parallelizable [20]. In this section, two different types of parallelization schemes, parallelization of the linear equation solver and task parallelization, based on distributed memory machine (IBM SP2) are presented for the 3D-FEM code anid the performance of each scheme is closely examined. For the implementation of efficient parallel FEM code, message passing paradigm using MPL (Message Passing Library) has been utilized. In the following section, the parallelization schemes for BiCG routine and the task parallelization are presented and the performances of two approaches are examined. The first strategy, the parallelization of the BiCG routine, is based upon the observation that more than 90% of the execution time of the FEM code is taken by the routine. As a result, development of an efficient linear equation solver can reduce overall computation time. On the other hand, the second approach, the task parallelization scheme, is evolved from the fact that the FEM is a frequency domain technique. When the frequency domain technique is used to characterize any high frequency interconnects, many repetitive calculations at each different frequency point are required, which is time consuming. However, the FEM system of equations for different frequency points are independent each other and can be solved simultaneously as presented in the next section. It has been found that the task parallelization strategy can provide an almost linearly scalable FEM code on a distributed memory parallel computer and it renders fast and efficient computation of the circuit parameters. 3.3.2 Parallelization Strategies (a) Parallelization of the Linear Equation Solver: [A][x] = [b] In most of the practical problems, the generation of the matrix [A] requires much less computational effort compare to solving the linear or eigenvalue equation. As

64 Frequency Sweep fl... fn Element Matrix Computation Matrix Assembly: Ax=b r --- —------- ---- ----- Linear Equation Solver: BiCG I Parallelization I --- —---------- -------- S-Parameters Output Figure 3.6: Parallelization of the linear equation solver, BiCG method, in the FEM scheme. a result, efficient and fast method of solving the equation can speedup the whole solution process. In this study, many efforts are devoted to implement parallelized BiCG method as illustrated in Fig. 3.6. As is well known, the BiCG method with diagonal preconditioning requires one matrix-vector multiplication in every iteration. The matrix [A] is computed once in the beginning of the solution process and used repeatedly during the iterations, while the unknown coefficient vector [x] has to be renewed in every iteration. Furthermore, to maximize the usage of the storage space and to minimize the number of floating point operations in the iterative solver, such as the BICG, only non-zero matrix entries and its integer coordinates are stored. For parallelization of the matrix-vector multiplication, the matrix [A] computed

65 by the master node is partitioned into smaller matrices [Ai] such that [A] = U [1 [A?] (3.7) where N is the number of nodes (CPUs), and each partitioned matrix [A,] is distributed to all the other slave nodes in the beginning of the BiCG procedure. Since for a given frequency the matrix [A] is fixed during the BiCG iteration process, it can be broadcasted once and for all. This is the step 1 in Fig. 3.7. On the other hand, the vector [x] which will be multiplied to the matrix [A] has to be distributed in the every iteration (step 2). As indicated in Fig. 3.7 each node performs partial matrix-vector multiplication simultaneously (step 3) and the results are collected to produce the final result (step 4). The following shows an example of broadcasting the matrix [A] using the message passing paradigm (MPL): CALL MP-BCAST(A,MsgLngth,Source,AllGroup) In a typical distributed memory machine, every node is connected each other using fast internal network, called switch in IBM SP2, and the internal network has its capacity (bandwidth) and resulting communication overhead. In our SP2 system, the bandwidth of the switch is 35 Mbyte/sec and the communication overhead due to the latency is 40 Lsec. Communication buffer size is also limited. As a result, when the size of the matrix or message exceeds a certain limit, the performance of the switch degraded remarkably. To overcome the communication overhead and acquire maximum speed, the size of the matrix should be sufficiently big. In the following example, multiple Send and Block Receive method is used to overcome the bandwidth limitation: IF (TaskID.EQ.O) THEN e- master node

66 Step I: Master Broadcasting [A [.A] [AN] 2... slave nodes Step 2.3.4: Master [X]\ Broadcasting [x] [ slave 2 N nodes [A1][x] [A2][x]... [AN][X] X A Gathering Master [A][x] Figure 3.7: Parallelization of the matrix-vector multiplication. Step 1: broadcast the partitioned FEM matrix [Ai]. Step 2: broadcast the vector [x] in every iteration. Step 3: perform partial matrix-vector multiplication [Ai][x] on each machine, where [A] = U [AJ] for i = 1,..., N. N is the number of slave nodes. Step 4: gather the results to form [b].

67 DO I = 1,NumOfTask-1 Destination = I Ntype = I CALL MP-SEND(A,MsgLngth,Destination,Ntype,MessageID) ENDDO ELSEIF(TaskID.NE.O) THEN +- slave nodes Source= 0 CALL MP.BRECV(A, MsgLngth,Source, TaskID,MessageID) ENDIF Note that to ensure the synchronism between the source (master node) and destination (slave nodes), the block receive subroutine (MP-BRECV) has been used. After the [A] has been broadcasted, the vector [x] is also broadcasted in a similar way in every iteration. Now, each small portion of the matrix-vector multiplication is carried out in each node simultaneously and the final results are gathered to form the resulting vector [b]. In the following an example of the parallel pseudo-Fortran code for matrix-vector multiplication is shown: N = INT(MatrixSize/NumberOfTask) Nstart = TaskID*N + 1 Nstop = Nstart + N - 1 DO I = Nstart,Nstop -- matrix-vector multiplication at each node VI = [A,][x] ENDDO CALL MP-GATHER(VI,b,BlockLength,Destination,AllGroup) In the above example, the sparsity of the matrix [A] can be fully utilized by using the

68 Speed up 1600 1400 --,- M= 38109 1200 --- --- M= 76709 1000 - ----------— o ---- ---- 800 600 - 400 200 0 2 4 6 8 10 12 14 16 18 20 Number of Processors Figure 3.8: Speedup curves of the parallel matrix-vector multiplication for two different problem sizes. M indicates the number of unknowns. indirect index scheme with possible adverse effect such as inefficient memory/cache access. Fig. 3.8 illustrates the performance of the FEM code with the parallelized linear equation solver for two different problem sizes M. As can be observed in the figure, the larger problem size has better performance improvement compare to the smaller one. This is mainly due to the communication overhead and task size, that is, the advantage of the "divide-and-conquer" strategy is offset by the communication latency for M=38190 case. (b) Parallelization of the Tasks The previous approach has gained marginal performance improvement due to the limited bandwidth and communication latency. To overcome these difficulties,

69 task/frequency parallelization strategy is explored in this section. This strategy is based upon the property of the frequency domain method such that the system of linear equations of each frequency point is independent each other. As a result, for a given problem the set of linear equations for different frequency points can be constructed and solved independently and simultaneously to obtain the complete spectrum of the scattering parameters over the broad range of frequency. In this scheme the computational tasks which correspond to different frequency points are distributed among the different nodes, and the solutions are found simultaneously. Indeed, as shown in Fig. 3.9, each slave node is assigned to have different set of problem corresponding to different frequency and the element matrix generation, matrix assembly and the solution process of the matrix equation using the BiCG are performed concurrently on each node. This type of parallelization scheme is often called "embarrassingly parallelized" version of the FEM code, since it only requires minimum parallelization skill, but can provide perfect and almost linearly scalable parallelization performance as illustrated in Fig. 3.10. It is worth mentioning that this approach is equivalent to using many separate computers in a fast network without requiring any communication between them. This would not be considered as a real parallelization in a strict sense, but it can be beneficial to the users of the parallel computers and provides a truly scalable parallelization strategy. Following shows an example of the pseudo-Fortran code for the task parallelization scheme: DO I = 1, NFRQ, NUMTSK IFREQ = TASKID + I IF(IFREQ. GT.NFRQ) STOP

70 Geometry Definition Pre-processing: Mesh Generation — " Frequency Sweep> Parallelization 7- - - -- - - - - - - - - - - - - -- ------- - ----- - - ^'- fl... fn.,/ ------------- - ------ f1i f2 Element Matrix Computation fn-i fn ------------- r --- —------- ----------- ------------------- I I - - - - - - - - - - - - - I - - - ----- I ',L i S-Parameters ' i, t t L -O.u-...... Output Figure 3.9: The concept of the task parallelization.

71 Speed Up F0 E 40.0;35.0 30.0 25.0 20.0 15.0 10.0 5.0.0 0 10 20 30 40 Number of Processors 50 Figure 3.10: Speedup curve for the task parallelization. (Serial FEM code for each frequency point) ENDDO 3.3.3 Discussion To find more effective parallelization strategy for the frequency domain finite element method, two different approaches are studied and the performance of those approaches are measured in terms of total amount of time required to solve a given problem. With the first approach, the parallelization of the matrix-vector multiplication, two different problems having size M equals 38109 and 76709 are solved. As shown in Fig. 3.8, simply increasing the number of CPUs does not guarantee

72 reduced execution time, especially for the smaller problem size. On the other hand. increasing the number of nodes for the larger problem size renders slight performance improvement with diminishing return after a certain point. The number of nodes for maximum speedup depends on the problem size as clearly shown in this study. Note that using more slave nodes over the optimum number causes degraded performance due to the heavy communication overhead and limited li buffer size, and thus offsets the advantage of using multiple nodes. Hence, the optimum number of CPUs for a given problem size depends on the capacity of the communication network, CPU speed, the performance of the switch (bandwidth), and the buffer size. Having the marginal perofomance improvement with the first approach, the frequency parallelization strategy is devised and its performance is fully examined. As discussed in the previous section, this approach is proved to give truly linearly scalable parallel FEM code due to the negligible amount of communication between the processors. Its simplicity and minimum required parallelization effort are the main advantages of this approach, and even more importantly, this strategy allows scalable parallel code and can be applied to many other applications. 3.4 Performance of the Artificial Absorbers In this section, the performance of an artificial absorber presented in section 2.10 is examined in the context of the finite element method for characterization of microwave and millimeter-wave circuits having inhomogeneous dielectric configurations. As an artificial absorber, a medium with designed complex permittivity and permeability is used to introduce a controllable amount of loss. Such an absorber can be effectively placed near circuit discontinuities to absorb the excited dominant and higher order modes as well as surface waves propagating in the substrate. In order to

73 minimize the unwanted reflections at the dielectric-absorber interface, the absorber is designed to provide impedance matching at the interface allowing the waves to be effectively transmitted through the interface and subsequently attenuated while they propagate into the absorber. To validate the concept of impedance matching between an inhomogeneous dielectric half space and an artificial absorbing layer in various situations, six different circuit configurations are examined including microstrip, CPW and waveguide structures. It has been found that a judicious choice of the artificial absorber makes it possible to reduce reflections less than -40 dB. 3.4.1 Microstrip Thru Line A 50 Q microstrip line is designed on 635 em thick substrate and the computational domain is terminated with PEC surfaces as shown in Fig. 3.11. The microstrip line is terminated with two absorber layers to simulate infinite transmission line. The position of the mesh termination surfaces is chosen to prevent the excitation of a waveguide mode in the frequency range of interest (low frequency region) and the side walls are placed 4H1 apart from the microstrip edge. As a result, the microstrip line has only dominant quasi-TEM mode below 22 GHz. Above 22 GHz, however, the overall structure starts to support the dominant as well as higher order waveguide modes depend on operating frequency in addition to the microstrip mode. The absorbers are designed to have two parts. The upper portion of the absorber is chosen to have c1 = jul = 1.0 + j3.0 so that it is matched to the air side. The imaginary part of the absorber is chosen to introduce enough losses for the wave propagating in the absorber. For the lower layer of the absorber, the material constants are chosen as 62 = 3.4 +jlO.0 and u2 = 1.0 +j2.9412 so that its intrinsic impedance is matched to the substrate impedance (772 = o/W= br)

74 0 -10 -20 H-30 - -40 -5O -60 I.... 10 15 20 25 30 Frequency [GHz] Figure 3.11: Computed return loss for the microstrip line terminated with artificial absorbers. Geometrical and material factors: L = 6.38 mm, W = 1.3 mm, D = 5.0 mm, H1 = 0.635 mm, H2 = 3.365 mm, ~r = 3.4,,r = 1.0, 61 = 1.0 + j3.0, /1i = 1.0 + j3.0, E2 = 3.4 + j10.0, and J2 = 1.0 + j2.9412 As observed in Fig. 3.11, the reflection from the absorbing layer is smaller than -40 dB up to 22 GHz. Considering the electrical thickness of the absorbing layer (D = 0.27Ag at 10 GHz) and the hybrid nature of the excited fields, the performance of the absorber is excellent in very broad frequency range. The gradual increase in the IS111 is mainly due to insufficient discretization of the FEM mesh at high frequency end of the spectrum. The dependency of the reflections on the discretization is examined further in later sections. 3.4.2 Microstrip Junction Terminated with Absorber To measure the performance of the absorbing layer in the close vicinity of circuit discontinuities, the absorber is placed at three different positions from the microstrip

7 impedance junction as shown in Fig. 3.12. The thickness of the absorber T is fixed to be 0.16Ag (Ag at 10 GHz) and the return loss is computed while varying the distance D from 0.08A9 to 0.32Ag. As presented in Fig. 3.12, the computed IS111 data for three different absorber positions agree very well with the data derived by the finite difference time domain method within 1.0 dB in the whole frequency region. It is quite remarkable to obtain the above accuracy with isotropic absorbing layer having only 0.16Ag thickness and placed at a distance less than O.lAg from the discontinuity. It is worth mentioning that for a given circuit if only IS111 data is needed to be computed, all the other ports can be terminated by placing appropriate absorbers and measure the reflection coefficient at; the designated port, resulting in significant reduction of required computational cost. 3.4.3 CPW Series Stub In this section, the performance of the isotropic absorber is examined under nonTEM type field distribution. In particular, a CPW series stub as shown in Fig. 3.13 is designed and absorbers are placed at four different places to evaluate the performance of the absorber. Since CPW supports non-TEM mode and its discontinuity generates evanescent and surface waves, placing the absorber in the close vicinity of the discontinuity can indicate the performance of the absorber under non-ideal condition. To simulate an infinitely long CPW thru line, absorbing layers are designed and placed at the end of the line as a matched load. The thickness of the absorber is fixed to 0.24Ag (Ag at 20 GHz). The return loss is calculated for the CPW series stub as a function of the position of the absorber (D = 0.093Ag, 0.213Ag, 0.334Ag, and 0.574Ag. As shown in Fig. 3.13, the calculated reflection coefficients converge

76 0 -5 I I I I I I I I.. -.I -10 F. --- —-— ] T 0 FDTD FEM: D = 1.5 mm - FEM: D = 3.0 mm FEM: D = 6.0 mm.... I.... I.... I. 15 -20 -25 10 15 20 Frequency [GHz] 25 30 Figure 3.12: Computed return loss for the microstrip junction terminated with artificial absorbers. Absorbers are placed at the side and top walls also in addition to the termination end. Geometrical factors: W1 = 1.3mm, W2 = 0.6mm, and T = 3.0mm. The thickness and the dielectric constant of the substrate are H = 0.635mm and ET = 3.4. The absorbers at the end are chosen as: E1 = 1.0 + j15.0 and f,u = 1.0 + j15.0 for the air side and E2 = 3.4 + j15.0 and UL2 = 1.0 + j4.412 for the dielectric side.

I i 1.0.8.6.4.2.0 16 17 18 19 20 21 22 Frequency [GHz] 23 24 25 Figure 3.13: Computed return loss for the CPW series stub terminated with artificial absorbers (closed structure) at 4 different positions. Geometrical factors: G = 0.225mm, S = 0.45mm, L = 1.35mm, W = 0.225mm, and T = 2.0mm. The thickness and the dielectric constant of the substrate are H = 0.635mm and Er = 9.9. The absorbers at the end are chosen as: e1 = 1.0 + j15.0 and l1 = 1.0 + j15.0 for the air side and 62 = 9.9 + j20.0 and /2 = 1.0 + j2.02 for the dielectric side. rapidly as the distance between the CPW series stub and the absorber increases. Compare to 0.1lAg of the previous microstrip case, the absorbers have to be placed slightly further (0.2Ag) in CPW case due to the non-TEM type field distribution. It should be noted that a proper position of an absorber to simulate an infinite line or non-reflecting surface is dependent on the field distribution as can be seen in the case of the microstrip junction and CPW series stub. 3.4.4 CPW Series Stub: Open Structure For the simulation of an open structure, the absorbing layers are placed on top and side walls of the circuit in addition to the termination end as shown in Fig. 3.14.

78 1.0.9.8.7.6.5.4.3.2.1.0 16 17 18 19 20 21 22 23 24 25 Frequency [GHz] Figure 3.14: Computed return loss for the CPW series stub terminated with artificial absorbers (open structure). Geometrical factors: G = 0.225mm, S = 0.45mm, L = 1.35mm, W = 0.225mm, D = 1.775mm and T = 2.0mm. The thickness and the dielectric constant of the substrate are H = 0.635mm and ~e = 9.9. The absorbers at the end are chosen as: 1 = 1.0+j15.0 and /jl = 1.0+ j15.0 for the air side and 62 = 9.9 + j20.0 and [2 = 1.0 + j2.02 for the dielectric side.

79 In this example, the absorbing layer having thickness T (= 2.0 mm = 0.32A9 at Jo = 20 GHz) is placed in the distance D (= 1.775 mm = 0.28A) from the dis — continuity. The distance D is chosen to minimize the computational domain while maintaining the required accuracy (refer to Fig. 3.13). The computed results reveal excellent agreement with the data obtained from the space domain integral equation technique [147]. Even though the thickness of the absorber is 2.0 mm in this case, it can be reduced further by increasing the amount of loss in the material. Increasing the mesh density also contributes to the improvement of the accuracy due to the smaller numerical errors. 3.4.5 Waveguide Absorber To verify the concept and the performance of the isotropic waveguide absorber, the design equations of the absorber presented in section 2.10 is implemented. The designed isotropic material is placed in a waveguide as shown in Fig. 3.15 and the dominant TE1o mode is excited. In contrast to the absorber used in the previous sections, the permittivity and permeability of the waveguide absorber are a function of frequency and excited mode in the waveguide.. The amount of reflections from the absorbing layer is measured using the standing wave patterns in the waveguide similar to the microstrip and CPW cases. As shown in the figure, the reflection coefficient in the whole frequency range is below -30 to — 40 dB except in the immediate vicinity of the 3.15 GHz which is the cutoff frequency of the dominant TE10 mode. The gradual increase of |S111 in the higher frequency region is mainly owing to the insufficient discretization. Since the amount of reflections at the absorber interface is zero, one can further reduce the numerical reflections at the interface by increasing the mesh density and by adopting smooth

80 transition of mesh in the computational domain especially at the material interface. In our simulations the mesh sizes in x-, y-, and z-direction (Sx, 8y, and Sz) are varied in the absorber region as indicated in the figure. As can be observed, the amount of reflected waves at the absorber interface is decreasing as the mesh size in each direction decreases, in particular, in x-direction (Sx). As another example, to investigate the performance of designed waveguide absorber in close vicinity of discontinuities in a waveguide, a 3-dimensional obstacle is placed inside of a waveguide as illustrated in Fig 3.16. As shown in the figure, the width (d1), length (d2), and height (d3) of the conducting cube are 9.0 mm, 9.0 mm, and 6.0 mm, respectively, and it is placed at the center of the waveguide. The waveguide has been operated in the dominant mode region (fTE10 = 6.56 GHz and fTE20 = 13.12 GHz) and the absorber is placed at the end of the waveguide at 11.0 mm (= 0.168 Ag at 8 GHz) away from the cube. The computed result using the FEM is compared with the measurement in [148] and reveals very good agreement. Note that the thickness of the absorber is chosen to be 5.0 mm (= 0.076Ag at 8 GHz) and Ax, Sy, and Sz in the absorber layer are 0.5 mm, 1.386 mm, and 1.04 mm, respectively. 3.4.6 Summary In this section, the application of the artificial absorber to microwave and millimeter-wave structures have been investigated. The design concept of the artificial absorber through the matching of the intrinsic impedances has been found simple and effective. Due to the loss mechanism of the absorber, the artificial absorbing layer can suppress higher order modes as well as the dominant mode in the guided wave structures. Furthermore, incorporation of the concept of the absorber

81 x ~rl, trl Er2, lr2 E, H I p | B tE, 6t 0 -10 -20 -30 -40 -50 ~ * ~ ~~I I I I[I I I I I I I E * 5x = 0.5 mm, Sy = 1.7 mm, 5z = 1.38 mm ----- = 0.5 mm, Sy = 1.19 mm, 5z = 1.1 mm - x = 0.2 mm, 6y = 1.19 mm, Sz = 1.1 mm fTE 10 TE20 fTE fC fC Cf o.... [..I.... ]..... 3 4 5 6 7 Frequency [GHz] Figure 3.15: Computed reflection coefficient of the artificial absorber in a waveguide designed for absorption of the dominant TE1o mode as a function of mesh size (Sx, 5y, and 6z) in the absorber region. Geometrical and material parameters of the waveguide: a = 47.6 mm, b = 22.0 mm, rl1 = 1.0, and url- = 1.0. The thickness of the absorber, 6t, is 0.041A9 (5.0 mm) at 4.0 GHz. fTElo, fTE20 and fTE01 are the cutoff frequency of each mode.

82 E, H b 0 -2 -4 -6 -8 -10 8 9 10 11 12 13 Frequency [GHz] Figure 3.16: Computed transmission coefficient for a conducting cube of dimensions dl = d2 = 9.0 mm and d3 = 6.0 mm. The width and height of the waveguide: a = 22.86 mm and b = 10.16 mm. The thickness of the absorber at the end of the waveguide is fixed to 5.0 mm which is 0.06Ag at 8.0 GHz.

83 into the existing FEM code is straightforward and does not require additional effort. The use of the isotropic absorber does not impose any degradation of the matrix condition number in contrast to the usual ABC's. Application of the absorbing layer also renders the possibility of alleviating the computational burden by reducing the computational domain. 3.5 Conclusions This chapter has presented the validation of edge-based finite element method using the generalized eigenvalue problems for lossless cavity. The eigenvalues for homogeneous and inhomogeneous cavities are found accurately within 1% error, and more importantly no spurious solutions are observed. After the validation, the FEM code is parallelized on the distributed memory machine, the IBM SP2, and two different types of parallelization strategies are examined. In particular, the task parallelization scheme is found to give linearly scalable performance improvement. Finally, a design procedure of the isotropic artificial absorbing layer is presented and its performance in various MMIC and waveguide applications are thoroughly investigated. The absorbing material is crucial for open boundary problems frequently encountered in many MIC/MMICs and waveguide applications.

CHAPTER IV PLANAR/NON-PLANAR HIGH FREQUENCY INTERCONNECTS 4.1 Introduction As has been pointed out by Montgomery [149] in 1989, Microwave Integrated Circuit (MIC) and Monolithic Microwave Integrated Circuit (MMIC) packages remain a performance and cost challenge in many applications. In the past, the microwave industry has treated the packaging of circuits as a necessary and rather obvious step to be taken after the completion of the circuit design cycle. Such an approach was heavily concerned with the thermal and mechanical properties of the package, but ignored important electrical effects which could dominate circuit performance. As a result, package interference made design cycles very costly and many times unsuccessful. It has been only recently that current R&D in industry and government laboratories has initiated an effort to address the critical problems of packaging in a way which links electrical performance with thermal and mechanical properties [150]. MIC and MMIC packages that are capable of good performance at frequencies as high as 60 GHz need to have small volume, low weight, microstrip or coplanar line compatibility, and most important, they should exhibit negligible electrical interference including low insertion and return losses. 84

85 For these packages to acquire some of the characteristics described above, special provisions need to be made during circuit design and layout stages. Specifically, in order to achieve lower weight/volume and reduce cost, designs lead to block and layered configurations where each layer performs a separate function. This design approach results in high density circuits where a large number of interconnects are printed on electrically small surface areas and communicate through the substrate in a direct through-via fashion or electromagnetically through appropriately etched apertures. In a circuit environment of this complexity, parasitic effects such as radiation and cross talk are intensified, thus making the interconnection problem very critical [151]. In order to reduce these effects with circuit layout compensations, accurate modeling tools are needed which can provide the exact correlation between circuit parasitics and geometrical characteristics. The establishment of a multi-layered multi-task monolithic circuit technology presumes a computer-aided design capability far more powerful and versatile than any presently commercially available. Furthermore, it requires a combination of tools varying from static to dynamic and from time to frequency domains, and the needs for a dynamically reconfigurable computational environment becomes more and more important as advancing technologies lead to increasingly complex geometries. As a response, this study presents frequency domain characterizations of a few of the most commonly found planar/non-planar interconnect elements in a multilayered circuit environment. Vertical transitions through apertures have been characterized in the past [151]-[156] and have been found to provide an efficient broadband bandpass RF connection. However, these transitions, due to their bandpass character, cut any low frequency or DC component which may need to propagate vertically to other layers, along with the RF signal as happens in detectors and mixers. In

86 such cases, direct interconnections through via holes are introduced along with their parasitic effects which, if not well understood, can degrade electrical performance to unacceptable levels. One example of such a simple but yet vital interconnection is a via used for ground potential equalization in monolithic circuits where more than one ground planes is utilized. In many configurations, grounded CPW lines are unintentionally created and thus lead to the excitation of parallel plate waveguide modes which, strong as they are, couple various parts of the circuit and introduce high levels of parasitic radiation [157]. One technique to suppress these modes is the use of vias which connect the ground planes printed on opposite sides of the substrate layer. In practical implementations, more than one via is needed, perhaps as many as five to ten per guided wavelength, leading to long structures of inductive elements periodically placed around the circuitry. While this configuration nearly equalizes the ground potentials, it can also excite Floquet waveguide modes that tend to propagate along the via structure and couple to the rest of the circuitry strongly and destructively. Such waves can turn a well designed band pass filter into an element with unrecognizable characteristics. Many recent efforts have emphasized the analytical and numerical evaluation of via inductance and capacitance using a variety of quasi-static and full wave techniques [131], [158]-[169]. In many cases, an equivalent circuit for the via is constructed through a scattering parameter evaluation, while in other cases, the inductance is computed from the total magnetic energy stored around the via or from the use of an excess inductance. While all of these techniques provide reasonable solutions, they disagree in the predicted values by as much as 10 to 30%. Unavoidably, this brings up the issue of which definition is more valid and. for that manner, more effective in

87 providing representative values. Herein, vertical transitions of via-type are analyzed in the frequency domain and the performance under sinusoidal excitation is studied to provide a complete understanding of their electrical characteristics at low as well as high frequencies. In this study, the problems of printed circuit board via holes are analyzed through both the scattering parameters and the flux/current methods, and the results are compared in an effort to establish the relationship between the two techniques and evaluate their contributions [21]. In addition, several different types of planar/nonplanar transition geometries including microstrip-to-CPW, CPW-to-slot, and CPWto-waveguide are fully characterized. The modeling of these structures is performed with the edge-based 3-D finite element method and compared with the time domain FDTD data. 4.2 Microstrip Short Circuit Fig. 4.1 shows a microstrip line printed on a 250,m GaAs substrate. This line is short-circuited through a via hole which is etched in the substrate and then metallized to provide a DC connection to the ground. For the sake of simplicity the geometry of the via is considered cylindrical and not pyramidal as it is in high frequency monolithic applications. Fig. 4.2 shows the electrical characteristics of the microstrip short circuit vs. frequency evaluated by the FEM and demonstrate very good agreement with the FDTD data. The slight discrepancy between the values can be attributed to numerical errors associated with both techniques. As shown in Fig. 4.2, the via provides good ground to the line up to 5 GHz and fails gradually at frequencies above that. Using the derived fields in the frequency domain, an equivalent inductance for

88 Top View -------------------------- -- 2 3 Y Side View Cross Sectional View - - - - - - - - - - - - - - - - - - - - - - - - - - -1 r - - - - - - - - - - - 600, 182 H 2 1, z" I "- — 3 -- W3 Figure 4.1: Geometry of a microstrip short circuit via hole at the center of the line. The dashed lines represent PEC walls to terminate FEM meshes. W, = 85um, W2 = 600,m, W3 = 1785,m, H, = 250,m, H2 = 600OOm, c1 = 12.9, and 2 = 1. The dotted lines around the geometry represent the outer boundaries for FEM mesh termination.

89 -10 - E - -5 ' -20 -30 - / o S2: FEM ---- S,: FDTD -40 - S21: FDTD -50.... 0 5 10 15 20 25 30 Frequency [GHz] Figure 4.2: Scattering parameters for the via. the via can be derived using two different approaches. One approach is based on the flux/current definition of inductance and is applied, herein, to the fields derived by the FEM. The other approach is based on the derivation of an appropriate equivalent circuit by fitting the derived scattering parameters in the frequency range of interest [168, 170]. For the geometries under study, a cascade of T-equivalent circuits, having two series capacitors and a shunt inductor, is used to model the discontinuity. It is found that the values of the capacitors are negligible (- 10-14 F), thus reducing the equivalent circuit to a shunt inductor. This equivalent circuit approach is applied to scattering parameters derived by the FDTD method. The inductances derived by these two approaches are shown in Fig. 4.3 and exhibit very good agreement. Let us now view the structure of Fig. 4.1 as two microstrip lines connected to the same ground via. It can be observed from Fig. 4.2 that [S211 of the microstrip short circuit increases with frequency, thus leading to unwanted coupling between the two

90 200 C '. -o-.- FEM: Flux/current Method 50........ FDTD: 0.1 - 20 GHz Fit c 100 50 0 I I 0 100 200 300 400 500 Height of Via Hole [gim] Figure 4.3: Inductance of the microstrip short circuit with via hole. The inductance is calculated using two different methods. In FEM, the flux density around the via is used (f = 10 GHz). In FDTD, however, the via inductance is extracted by fitting the derived S-parameters to a T-equivalent circuit. microstrip line sections. In particular, the magnitude of the coupling coefficient S21 increases from -30 dB at 1.0 GHz to -4 dB at 30 GHz. The failure of the via to provide a good DC connection at higher frequencies is intensified by increasing the via height, as shown by the variation of the equivalent inductance and of the scattering parameters vs. via height (see Figs. 4.3 and 4.4). The frequency dependence of the via characteristics shows up also in the distribution of the electric fields and the associated current flow, as illustrated in Fig. 4.5. At lower operating frequencies, where the via presents an almost perfect ground, the current flows from the microstrip line to ground through the wall of the via facing the feeding line. Under these frequency conditions, an almost zero current flow is observed toward the other connecting line, providing very high isolation between the

91 0.,.. I...... I '.'.'.' -5 0 5 - 10 -<r - 21 -15 S2 H = 0.3 mm: FEM 15 21 1 --- S21 H = 0.2 mm: FDTD S21: H = 0.3 mm: FDTD -20.... 10 15 20 25 30 Frequency [GHz] Figure 4.4: Scattering parameters of the microstrip short circuit for two different substrate heights (H1 = 0.2, 0.3 mm) with via hole grounding at the center of the line.

92 (a) (b) Figure 4.5: Magnitude of the electric field distribution at f = 20 GHz under the microstrip short circuit (z = 0.19 mm plane). The edge singularity and the standing wave patterns of the electric field distribution are clearly revealed.

93 x Top View IW 1 ------------------------------------------------ Side View ----------------------------------------------- ~2 l 2 I l I I 00W I I::I H1 Figure 4.6: Geometry of a 1-port microstrip ground pad. This is a board level via. Wo = 2.30 mm, W1 = 1.15 mm, W2 = 2.88 mm, W3 = 10.0 mm, H1 = 1.0 mm, H2 = 2.0 mm, el = 3.4, 62 = 1. The dotted lines around the geometry represent the outer boundaries for FEM mesh termination. two lines. As the operating frequency increases, current flows to ground from the two side walls of the via, with considerable current flowing towards the other microstrip line, resulting in substantial coupling between the two lines. 4.3 Microstrip Ground Pad In many high frequency applications, ground connections on substrate surfaces opposite to the ground plane are commonly called for and need to be carefully characterized. These IC/RF ground connections can be realized with the use of ground pads, as shown in Fig. 4.6. In general, the characteristics of the vertical via hole

94 600 I. 500 400 I 300 CZ / 200 100 0.....0.5 1.0 1.5 2.0 2.5 Height Hl [mm] Figure 4.7: The via inductance calculated from the flux/current definition at 10 GHz as a function of via height. and the ground pad can be represented by an appropriate lumped element. In this study, the vertical via hole is modeled by an inductive element in view of its role as a vertical current path from the microstrip line to the lower ground plane. The inductance of the via hole has been computed in a wide frequency range and, without loss of generality, the inductance is computed at 10 GHz as a function of via hole height (see Fig. 4.7) using the well known flux/current relation applied to magnetic field distributions. One important issue associated with the high frequency characterization of distributed components, such as the ground pad, is the use of an appropriate reference plane in the evaluation of the scattering parameters. The choice of this reference plane greatly affects the derivation of a meaningful equivalent circuit. As shown in

95 210. 180 - " ' ' -"' -- o — - -o — - -o --- -.... 150 - -- " - - 120 - 90 -60 - -o- Dy=-3.24mm -:: -—. ---. Dy=-2.88 mm 30-. — 3 — Dy=-1.44mm 0 - - -- Dy=-0.44mm - -0 — Dy = 0.0mm -30 ----- Dy = +0.41 mm -60,, 4 6 8 10 12 Frequency [GHz] Figure 4.8: Phase of the S11 of the microstrip ground pad computed using the FEM at several different reference planes (Dy = +0.41, 0.0, -0.44, -1.44, -2.88, and -3.24 mm from the top to the bottom of the graphs). Dy is 0.0 mm at the geometrical center of the via. H1 = 1.0 mm. Fig. 4.8, the phase of the scattering parameter S1l always varies with frequency, irrespective of the position of the reference plane, in a way that makes it very hard to find an equivalent inductance which is both frequency and reference plane independent. One might need to consider a more complex circuit arrangement to successfully fit the derived parameters over a wide frequency range. As a result, the use of the scattering parameters for the derivation of a simple equivalent circuit becomes much less successful, as compared to the previous approach which is based on the flux/current definition. The variation of the inductance, which is computed based on the flux/current

96 700 ' 600 l H =1.5 mm 500.... --- H = 1.0 mm -— A --- H = 0.5 mm; 400 o l l 300 200 - --------— E --- -------------------— E -- 3 5 7 9 11 13 15 Frequency [GHz] Figure 4.9: The via inductance of the microstrip ground pad calculated from the flux/current definition for three different via heights (Hi = 0.5, 1.0, 1.5 mm) as a function of frequency. definition, vs. frequency for various via heights is shown in Fig. 4.9. As seen from these results, for substrate thicknesses which are approximately less than one fifteenth of the guided wavelength, the inductance is independent of frequency, indicating that the quasi-static predictions could be very accurate. The advantage of using partial differential equation techniques, such as the finite element method, is found in their ability to compute fields and power/energy distributions around the circuit and its discontinuities. The visualization provided by the FEM technique can enhance the qualitative understanding of the electrical perfrmane of the r rcomponents under study. For example, Fig. 4.10 reveals the electromagnetic energy distribution under the microstrip line at a plane 0.9 mm above the ground. As seen from the these

97 0,9 0.8 0.7 0.8 =l za ( = 0 u e hle = 0.8 0.4 0.3 0.2 0.9 0.8 0.7 field distribution are clearly revealed0. Figur...1

98 Magnetic Field Density on the Ground Plane 8O 80 20 40 60 80 100 10 GHz. The via inductance is mainly determined by the magnetic energy surrounding the via. the sides and back of the vertical via holes. In the same figure, the standing wave to the pad are clearly revealed. From the magnetic field distribution around the vertical via hole, the frequency dependent via inductance has been calculated. As hole, along the shortest current path from the line to the ground plane. 4.4 CPW-to-Microstrip Through-Via Transition 4.4.1 Without Via Hole substrate incorporate printed circuit components which may contribute to similar

99 Top View w Side View (Cut along the center) ----------------------------------- -We CPW on top surface - CPW/MS ground, Microstrip on bottom surface ---------------------------------------------- H2 H1 3 H3 Figure 4.12: Geometry of the transition from CPW to microstrip on different layers without via hole. Wg = 50um, W, = 75,/m, Wo = 75bm, We = 200/m, W = 1000m, H1 = 1Om, H2 = H3 = 400/um, and c, = 12.9.

100 I I I I 0 * 000000 -100 I-~ ~0- C) -- S2 FEM -20- --------:FEM o S21: FDTD * Sl: FDTD -30. —. —.-.-.-..-,... 20 40 60 80 100 Frequency [GHz] Figure 4.13: Scattering parameters for the transition from CPW-to-microstrip on different layers without via hole connection. electric functions or tasks. When this approach is incorporated into circuit design, the two substrate surfaces host complementary geometries to reduce cross-talk and to increase the degree of circuit integration. In this section, as an effort to use both sides of a substrate, a coplanar waveguideto-microstrip configuration is investigated. As shown in Fig. 4.12, the CPW line is printed on top of a substrate and the microstrip line is fabricated on the other side of the substrate surface. With this arrangement, the ground of the CPW can serve as the ground of the microstrip, while the CPW apertures and microstrip run in parallel and are separated by some distance to reduce parasitic cross-talk. Given the geometrical details, the CPW transmission line has characteristic impedance 43 Q, while the microstrip line reveals 50 Q. These lines are fabricated on different sides of a 100 jum thick GaAs substrate. As shown in Fig. 4.13,

101 the computed cross-talk (|S211) between the lines is near -15 dB at 20 GHz and decreases further as frequency decreases. In other words, more than 95% of the in — put power is reflected back to the input port in the low frequency region, revealing relatively good isolation between the ports. However, as the operating frequency in — creases up to 100 GHz, this cross-talk increases monotonically up to approximately -4 dB. The increased cross-talk could be ascribed to several facts: high frequency radiation effects from the open end of the CPW and microstrip, or the higher order and substrate modes generated between the two transmission lines. 4.4.2 With Via Hole The CPW-to-microstrip transition structure studied in the previous section can be modified to have much less insertion loss by placing a via hole between the two transmission lines, as shown in Fig. 4.14. When those two lines are interconnected through the metallized via hole between the end of the microstrip line and the center conductor of the CPW, a wide band direct DC/RF connection is obtained. As a result, the transition structure can operate in an ideal fashion, at least in the low frequency region. The computed S-parameters of the transition are shown in Fig. 4.15. The return loss remains under -20 dB for frequencies as high as 40 GHz. The FEM calculated values agree very well with that of the ideal transmission line. Note that the theoretical value of the reflection coefficient of an ideal transmission line is given as -22.5 dB for 43 Q-to-50 Q impedance step. As the operating frequency increases up to 100 GHz, the performance of the transition starts to degraded, giving a return loss of about -9 dB. Even though the transition shows degraded performance in the high frequency region, the overall characteristics closely follow that of the ideal transmis

102 Top View W Side View (Cut along the center),- - - - - - - - - - - - - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - We CPW on top surface K-1 CPW/MS ground, e........... / Microstrip on bottom surface Rectangular via I --- —--------- ---------------------------- ------ H2 H1 H3 Figure 4.14: Transition from CPW to microstrip on dlifferent layers through via hole. Wg = 50/um, W, = 75/m, Wo = 75,m, We = 200y/m, W = 1000m H1 = 100m, H2 = H3 = 400/,m, and r =: 12.9. The vertical via is an empty rectangular cavity with PEC walls.

103 -10, * -20 - -" e= -,,,- -. * 2 -30 FEM -30 -------- Sl:FEM -40 S21: FDTD S*:FDTD -50,,... 20 40 60 80 100 Frequency [GHz] Figure 4.15: Scattering parameters of the CPW-to-microstrip transition with vertical via hole connection. sion lines in the low frequency region. This type of transition structure can be used in multi-layer MIC or MMIC interconnects, the hermetic sealing of a package, and antenna feeding networks. 4.5 Channelized CPW-to-Microstrip Transition In this section, another type of CPW-to-microstrip transition geometry, which is planar instead of the non-planar configuration in the previous section, is carefully designed and characterized for broadband applications. The planar CPW-to-microstrip transition has found many important applications in the area of hermetic packaging, high frequency probing, etc. Fig. 4.16 shows the transition from a channelized CPW [171] to a microstrip line and in the transition, both CPW and microstrip lines are printed on the same substrate surface, as opposed to the previous CPW-to-microstrip transition.

104 Top View ------------------- I — r I i I I I I II I I I I -Twf, I w 1W I. '19 I I II I I I I I I I ------------------- J. Side View (Cut along the center) I --- —-------------------------------------— I I I I I I I I I I I I I I I I Er H2 Hi Figure 4.16: Geometry of a channelized CPW to microstrip transition. Wg = 50gm, W, -- 75pm, W = 1000gm, Hi = 400gm, "2 = 500gm and cr = 12.9.

105..... I.... I.I.v I. I I I I. I I v. I I I I I... E v 0 O O O G ( O --- S: FEM -5 -- S,:FEM o S2: FDTD b Sl FDTD -10 ~ * * * * * -- --—. --- — ------ * —. - --- -15......... 0 10 20 30 40 Frequency [GHz] Figure 4.17: Scattering parameters for the channelized CPW-to-microstrip transition. The scattering parameters of such a transition geometry printed on GaAs substrate of thickness 400 ym are shown in Fig. 4.17. The insertion and return losses remain almost constant in the frequency range up to the 40 GHz region. The agreement between the FEM and FDTD results is very good, and further the results of both techniques agree very well with those obtained using the space domain integral equation method (SDIE) [172]. Also, compared to the value -12.7 dB computed from the ideal transmission line step discontinuity, 50-to-80 Q, the computed data give very close results regardless of its non-TEM propagating modes. Similar behavior has been found for the case of 100 /im GaAs substrate, except that in this case the return loss decreases to about -20 dB. Note that the return loss of the ideal transmission line having a 43-to-49.5 Q discontinuity is -22.5 dB. The reason for this good performance irrespective of the substrate thickness is that the field distribution around the center conductor of the chanellized CPW becomes

106 I Ex I distribution at h = 0.39 mm: f = 20 GHz 2-. 1-.. 0^ 0 A Fi p-0.5 8 10 -0.5 X [mm] Y [mm] I I I I I I I I I V.v I 0 I I I I I I I I I I I I I I IIIIIa.n,i i OI 0 1 2 3 4 5 6 7 8 9 10 Figure 4.18: Magnitude of transverse component of the electric field at z = 390am plane.

107 very similar to that of the microstrip (i.e., most of the fields are confined under the center strip). Fig. 4.18 shows the magnitude of the tangential (to the surface) component of the electric field at a plane just underneath the dielectric-air interface at 20 GHz. The standing wave patterns along the CPW and microstrip lines can be clearly observed. 4.6 Asymmetric 2-port Via with Plated-Through-Hole In a multi-layer circuit environment, the vertical interconnections between various signal lines in different layers are frequently encountered and need to be modeled for a thorough understanding of low as well as high frequency effects. As a matter of fact, the modeling procedure is focused on the stage of finding an appropriate lumped equivalent circuit, due to its usefulness and simplicity, despite its limited accuracy. In general, the vertical interconnections are fabricated using a plated-through-hole (PTH) in the ground planes located between the signal layers. Even though it is possible to have more than one ground plane and signal layer, in this section, we considered a simplified geometry which has two signal layers and one ground plane between them, and the two signal lines are connected through a vertical via hole. Fig. 4.19 has the schematic view of an asymmetric 2-port via hole with a PTH in the middle of the substrate. The PTH is filled with air, while the two substrates have the same relative dielectric constant (Cr = 4.5) and thickness (H1 = H2 = 0.2 mm). To examine the effect of PTH size on the transition performance, three different PTH opening sizes have been examined while the other geometrical factors remained unchanged. At first glance, it is expected that as the size of the opening D increases, the capacitance between the vertical via and the ground plane becomes smaller and the

108 Side View Through Via Upper Microstrip Line Top View PTH Lower Microstrip Line \ i0- -li ii:- ii i:~ f 0 -:: \0 —E1;0 -I- -I 0. - 0I.f -: - i i 0-:-0 ii if - - D -.40W-i~ft4Lyf f X 0 fiD 7 0. | E tEdE. y\ t i 44j 40 f 0 0 X:; - - - - - - X4 - —: ---::;i::: t -:: -: -: —:;;.i - r.i. ^l- j. -:.,:-^ i —Ar Up-permMit 2 I Li: Iw, Figure 4.19: Geometry of asymmetric two port via hole passing through a plated-throughhole (PTH) in the middle of the substrate. H1 = H2 = 0.2 mm, W1 = 0.4 mm, W2 = 0.8 mm, Erl = rE2 = 4.5, D = 0.2, 0.4, 0.7 mm (3 cases).

109.0 -5.0 -10.0 -15.0 on n.I I. I. I ' - -6-' - - i.-o ' D ~:- --— " —:" -"". -,-~ 4 I -)-S-I- 511: case2 -> D=0.2 mm ------ S 1: case3 -> D=0.4 mm ---- S11: case4 -> D=0.7 mm.1[~~. I. 1. - 5. 7. 9. 11. 13. 15. Frequency [GHz] Figure 4.20: Comparison of the S-parameters for 3 different PTH sizes. resulting reflection coefficient would be reduced. However, as shown in Fig. 4.20, the reflection coefficient increases as the distance between the via hole and the opening in the ground plane increases, which means there is more reflection with larger opening size. Considering the fact that increasing D implies less and less interaction between the via hole and the ground plane, the above results seem to be contradictory. However, as seen in Fig. 4.19 the Inicrostrip lines are connected to the vertical via through its pads, and there are impedance mismatches between the microstrip lines on the substrate and lines crossing the air gap between the via pads and the microstrip lines. While the characteristic impedance of the microstrip lines is 48 Q, that of the line across the gap has a much higher value. As a result, as the length of the high impedance line section becomes longer, the effect of the discontinuities becomes intensified. Therefore, it is inappropriate to use the above configuration to examine the effect of the PTH opening size. It is recommended, instead, that for the

110 100.0 I ' I~ I -o — H = 0.1 mm 75.0 -.. ---- H=0.15mm = 500 -.........0. - 25.0 -.0. I I I 5.0 7.0 9.0 11.0 13.0 15.0 Frequency [GHz] Figure 4.21: Via inductances of two different via height: H1 = H2 = 0.1 and 0.15mm and D = 0.4 mm. Inductances are evaluated from the magnetic flux density surrounding the via. characterization of the effect of the PTH sizes, the air region between the via and PTH has to be filled with the same dielectric material to minimize the transmission line discontinuities. The equivalent via hole inductances having two different heights (Hi = H2 = 1.0 and 1.5 mm) are computed from the magnetic flux density, as shown in Fig. 4.21. It is observed that the equivalent inductance of the via hole inside of the PTH has a slight variation as operating frequency increases and shows deviation from linear dependence on the via hole height. Fig. 4.22 shows a simple equivalent circuit (so called, type I) of the asymmetric 2-port via hole geometry. While the inductance values (Lvia) of the via hole are extracted from the flux model, the shunt capacitance is found by matching the scattering parameters, both magnitude and phase, as computed from the FEM with the equivalent circuit.

111 Type I Lvia Lvia T0 Figure 4.22: Type I equivalent circuit used for the asymmetric 2-port via hole geometry. Note that if the PTH opening becomes smaller, C increases and finally becomes a short circuit. Fig. 4.23 shows an example of the matching procedure. As shown in the figure, the magnitude and phase data reveal some degree of discrepancies due to the oversimplified equivalent circuit. For the transition with D = 0.2 mm and H1 = H2 = 0.1 mm, the inductance and capacitance values of the type I equivalent circuit are determined as Lvia = 25 pH and C = 0.28 pF for the best agreement between the FEM and the equivalent circuit data. It is interesting to observe that increasing the distance D causes a smaller capacitance value and finally renders zero capacitance, leading to an open circuit as D becomes infinity. This observation is consistent with the physical structure, since infinite D implies no ground plane. Similarly, decreasing D causes a larger capacitance and eventually makes the shunt capacitance a short circuit, which is what is happening in reality. 4.7 Double Via Hole In many microwave and millimeter-wave circuits, more than one via hole is used in close proximity to provide certain functionality, such as ground path, shielding, or vertical transition. However, these via holes can cause undesirable effects due to their

112 S11 0c 'O -1 c: -3 -4 S21 I - ----- - - -- -- -- -- T- -- ------------- II. — --- --- --- T^ --- —----------- --------------— t --- —----------- --------------— * --- —-----— _ ---I ) en 10 Frequency [GHz] S11 15 5 10 Frequency [GHz] 15 S21 I 100 ---------------------- (I) -0 0 --- ---— 1 --- —------ I >oso~io~opoaoasI 'a (I) 0n 100 q I I I - - - u o-o c (-o o-o - o-o -o oo 4 --- —---- -100: ) -100 5 10 Frequency [GHz] 15 5 10 Frequency [GHz] 15 Figure 4.23: Matching the S parameters using the type I equivalent circuit (Circle: FEM data, Line: Equivalent circuit). H1 = -2 = 0.1 mm and D = 0.2 mm. Lvia = 25 pH, C = 0.28 pF with type I equivalent circuit.

113 Top View W5 +W3 W4 D WI Side View D H-H - -^.. p Via Holes Figure 4.24: Geometry of double via hole. H = 1.6 mm, W1 = 1.4 mm, W2 = 1.0 mm (85 Q), W:3 = 0.2 mm, W4 = 1.0 mm, W5 -- 1.4 mm, En1 = 4.5, D = 0.2, 0.4, 0.6, 0.8 mm (4 cases) mutual interactions, that is, electromagnetic coupling, depending on the operating frequency and geometrical factors. As a result, these closely spaced via holes have to be modeled and used carefully to provide the desired properties and prevent spurious phenomena. For successful design of high performance high frequency circuits, it is imperative to accurately predict the mutual coupling of multiple via holes. In view of this, in this section, the effect of the electromagnetic coupling between two identical ground via holes are studied and an appropriate equivalent circuit is developed to depict the coupling effects. As shown in the schematic view of the structure in Fig. 4.24, two via holes are placed face to face each other with separation distance D. To model the coupling from one via hole to the other, the separation distance is changed from 0.2 to 0.8 mm

114.0 -i -i --- S21:D = 0.8 mm -— e --- S21: D = 0.6 mm ----- S21: D = 0.4mm -10.0 -- S21 D=0.2mm ' — - - — ~:.: -20.0 — v. - " -...." ---'"! -30.0 4.0 6.0 8.0 10.0 Frequency [GHz] Figure 4.25: Comparison of the S-parameters for 4 different separation distances. and the corresponding scattering parameters are computed as presented in Fig. 4.25. As can be observed in the figure, the cross talk, [S21, between two via holes gradually increases as the operating frequency increases and becomes stronger as the distance becomes smaller. The electromagnetic field distributions around the transition are shown in Fig. 4.26 revealing the standing wave patterns on the microstrip lines and its edge singularity. Also, as shown in the figure, more electromagnetic energy is coupled from one line to other as operating frequency increases. Derivation of an appropriate equivalent circuit representing the coupling effects between the via holes could provide a useful tool for design and optimization of high performance MIC and MMIC circuits, in addition to intuitive understanding of the electromagnetic phenomena. As a response, in this study, the type II equivalent circuit shown in Fig. 4.27 is derived for the geometry. Note that as the separation distance D becomes very large, the mutual inductance L12 and capacitance C12

115 0,. 0.8 0(a)7 0.5 0.4 0.2 (b) 0.9 0.8 0.7 0,4 0.3 0.2 0.1 0 (b) 10 GHz under the microstrip line at z - 1.5 mm plane. The edge singularity and the standing wave patterns of the electric field distribution are clearly revealed.

116 Type II Lvia Lvia o~ ------- -------— ~o L12 IC12 T p Figure 4.27: Type II equivalent circuit used for the double via hole geometry. The L12 and C12 represent the mutual coupling between the via holes. Note that if the separation distance D increases, the L becomes smaller and finally becomes ideal short circuit. become negligible, and finally become two separate circuits having Lvia connected to the ground plane directly. It is also interesting to observe that there are resonances in the high frequency region due to the shunt impedance elements and these frequencies are tabulated in Table 4.1. Fig. 4.28 shows an example of the S-parameter matching procedure, and in Table 4.1 the equivalent circuit parameters are listed for the four different separation distances. For all the cases, the via inductance Lvia is forced to have the same value, while the mutual inductance L12 and capacitance C12 are searched to find the best fit of the FEM data. As shown in the table, the capacitance values remain constant in all the cases, but the mutual inductance L12 changes as a function of separation distances. As the distance D increases, L12 decreases from 120 to 80 pH for D = 0.2 and 0.8 mm, respectively. The constant coupling capacitances and changing inductance reflect the fact that the coupling is mainly coming from the vertical current flow on the via holes and their resulting magnetic fields. The derived equivalent circuit can be used to model and predict the electromagnetic coupling among the multiple via

117 Si1 S21 A u -0.2 0) -0.4 'O. — c -0.6 0) -0.8 -1 (I... I UUOrc-,I 4 I ) v m'-10 0) 3 -20.,-. c 0) -30 An I r I I r I I "'1"""'T"""r'""'t"""~ I I I I.-,J-,,-,,-l,,,,,,, L,,,, —lIII I IV I I "-'T"""T"-"'T'""'~ 1 I I 1 r ) 2... 4 6 8 Frequency [GHz] -'tU 10 I,2 4 6 8 1C Frequency [GHz] S21 0) a) 0) s3 100 0 S11 ----- --. --- —--------. -- ___J... —._ —_..-_L-...... -_-_... i I I I 100 ) a) (D a) 0C O3 Q. 0 I I I - I I I I I I ' --- —------------------------—. I I I I I I 4 I II -100 -100 2 4 6 8 Frequency [GHz] 10 2 4 6 8 Frequency [GHz] 10 Figure 4.28: Matching the S parameters with the type II equivalent circuit. Lvia = 380 pH, L12 = 80 pH, and C12 = 0.1 pF. H = 1.6 mm, W1 = 1.4 mm, W2 = 1.0 mm, W3 = 0.2 mm, W4 = 1.0 mm, W5 = 1.4 mm, erl = 4.5, D = 0.8 mm. Circle: FEM data, Line: Equivalent circuit.

118 D [mm] I Lia [pH] L12 [pH] C12 [pF] fres [GHz] 0.2 380 120 0.1 45.94 0.4 380:105 0.1 49.12 0.6 380 90 0.1 53.05 0.8 380 80 0.1 56.27 Table 4.1: The equivalent circuit parameters of the double via hole structure and corresponding resonance frequency. The via inductance Lvia is extracted from the flux model, while the shunt terms, L12 and C12, are determined to match the FEM data with the equivalent circuit. holes often encountered in the MIC and MMICs. 4.8 CPW-to-Slotline Transitions With the advance of the slot antenna and its array technology, leading to a high efficiency and narrow beam antenna system [173], the effective feeding network of the antenna plays a critical role for its performance [174]-[179]. In this study, two different types of slot antenna feeding networks, planar and non-planar, are modeled using the finite element method. The first transition is designed for direct coupling of RF power from the CPW to the slotline, and printed on the same side of the substrate, which is easy to fabricate. In contrast, the second structure is devised for electromagnetic coupling from the grounded CPW to the slotline on the opposite side of the wafer without using bonding wires, and is thus more reliable for harsh environments. 4.8.1 Planar Transition The schematic diagram of a planar CPW-to-slotline transition with an air bridge is illustrated in Fig. 4.29. The characteristic impedances of the CPW and slotline are 60 and 70 Q, respectively. The input CPW line is terminated in an open circuit and an air bridge is used to connect the open end of the CPW center conductor

119 Top View A B (Symmetry Plane) A' B' A-A' cut B-B' cut W2 1 wL W2 -~*-<~-<-=~ 77 w3 __H___ Sr IITH Figure 4.29: Geometry of the CPW-to-slotline planar transition with air bridge: L, = 4.751 mm, S1 = 0.114 mm, S2 = 0.129 mm, W1 = 0.305 mm, W2 = 0.254 mm, W3 = 0.178 mm, H = 0.635 mm, and er = 10.5.

120 0.... I.... I.... 5 - E | ' 0.. I ~ -2~0. -30. SlI FEM --—.. S21:FEM -40..,. 5. 6. 7. 8. 9. 10. Frequency [GHz] Figure 4.30: Scattering parameters for the coplanar waveguide-to-slotline transition with air bridge. to the opposite edge of the slotline. The air bridge plays a role of equalizing the potential between the two conducting surfaces and suppressing any spurious modes. The thickness H of the substrate is 0.635 mm with Cr = 10.5, and the other side of the substrate is not metallized. The circular bend at the end of the slotline is intended to provide a smooth transition and the length Ls is designed to be - Ag(slotline)/4 at the center frequency, 7 GHz. As a result, the electric field has a minimum at the end of the slotline (short circuit) and a maximum at the transition region. In the modeling procedure, the circular bend of the slotline section is approximated to a rectangular geometry having the same length, and the curved air bridges are also replaced with infinitely thin piecewise linear PEC wires. Note that the FEM modeling is performed under a shielding environment. The performance of the transition is shown in Fig. 4.30. It exhibits very narrow band characteristics due to the

121 rectangular approximation and also shows periodicity of frequency response. 4.8.2 Non-Planar Transition A non-planar CPW-to-slotline transition, illustrated in Fig. 4.31, is studied in this section. The characteristic impedance of each transmission line section is 70, 80, and 70 Q for the CPW, microstrip (elongated CPW center conductor) and slotline, respectively. To provide a good short circuit effect, the open end of the CPW ground planes has been extended by length W4 which is Ag(microstrip)/4. In addition, to suppress the parallel plate waveguide mode between the two conducting surfaces, the CPW ground planes are trimmed to finite. The overlapping lengths of the elongated CPW center conductor and the underlying slotline are designed to be Ag(microstrip)/4 and \g(slot)/4, respectively, at the center frequency, 6 GHz. Since the normal electric field has a maximum value at the open end of the microstrip line, the magnetic field has its maximum at Ag(microstrip)/4 from the open end, thus providing maximum coupling from the microstrip to the slotline. For numerical modeling, the curved section of the slotline is approximated to a rectangular geometry for simplicity and the whole structure is placed inside of closed PEC surfaces (i.e., a shielding environment). On the other hand, the measurement was performed under an open environment. As shown in Fig. 4.32, the modeled and measured data reveal a reasonable correspondence, including the prediction of the two resonances in the frequency region of interest, even though the modeled geometry has narrow bandwidth and inferior performance due to the rectangular approximation. From the above observation, it becomes clear that to accurately predict the performance of the transition structures having circular geometries, the curved structures need to be modeled exactly. Compare to rectangular shaped geometries,

122 Top View A B w4 t1 wi 6 - - - - - - --- - - - I- - I I I I I L1, I I I I W2 - W3.Ws,l_ --- —__ --- —----- ~ ~ ~ ~ ~ ~ ~ ~ ~ _. ~ ~ __ ~ ( ~ A' B' A-A' cut B-B' cut W2 3 CPW Plane Wi 4 I I H W5 -Slot Plane Figure 4.31: Geometry of the CPW-to-slotline non-planar transition: L1 = 2.54 mm, L2 = 2.54 mm, L3 = 5.0 mm, L4 = 6.229 mm, W1 = 2.54 mm, W2 = 0.152 mm, W3 = 0.318 mm, W4 = 2.54 mm, W5 = 0.178 mm, H = 0.635 mm, and (, = 10.5.

123 0. * -.!, *4 -- ML -10. 7:o 0 0, -20. - -30 ---- S1 1:FEM -30. - o --- S21: FEM 0 S 1: Measurement * S21: Measurement -40.... 5.0 6.0 7.0 8.0 9.0 Frequency [GHz] Figure 4.32: Scattering parameters for the coplanar waveguide-to-slotline transition under shielding environment. The design frequency is 6 GHz. the smooth curved transitions exhibit wide band characteristics as revealed in this study. 4.9 CPW-to-Waveguide Transition In many millimeter and submillimeter-wave waveguide mounting structures and waveguide probe stations (refer to Fig.4.33), passive and active circuits are embedded into the waveguide [180]-[185]. In these structures, accurate prediction of the driving point impedance and input reflection coefficient is of critical importance to designers. However, the interactions between the lumped or distributed printed circuits and the waveguide are not well understood, and furthermore, the design of such structures heavily depends on experimental data and experience. Thus, the design of a waveguide probe station with maximum coupling from the probe into the waveguide is a tedious and difficult task. In this study, a CPW-to-waveguide transition typically

124 found in the waveguide mounting or probe structures is characterized using the open boundary FEM and its results are verified using the FDTD technique. As shown in Fig. 4.34, the waveguide probe is fed by a shielded coplanar waveguide and has the shape of a rectangular patch antenna. The probe is inserted into the waveguide through a slot. Note that the front end of the probe inserted into the waveguide is not connected to any waveguide walls. For the simulation of the infinite waveguide structure having more than one mode as frequency increases, an isotropic absorbing layer optimized for a dominant waveguide mode is applied1. The dimensions of the probe as well as the thickness and the dielectric constant of the substrate are very important to achieve maximum energy transfer from input port to the rectangular waveguide. In this study, the thickness of the dielectric substrate carrying the probe is chosen to be 2.0 mm with cr=12, and the width and length of the probe are 3.8 and 10.8 mm, respectively. This probe is designed to feed a WR-187 rectangular waveguide and the scattering parameters are computed in the 3.1 to 7.4 GHz range. Even though the probe is printed on top of dielectric substrate, the dielectric material under the probe can be removed, resulting in a membranesupported structure, to maximize the coupling and to minimize the dielectric loss, especially in the high frequency region. In the simulated frequency range, three different modes, TE1o, TEz and TEz1, are excited inside the waveguide with cutoff frequencies at 3.15 GHz, 6.30 GHz and 6.82 GHz, respectively. To characterize the probe performance, the magnitude of the reflection coefficient ISn11, as well as of the coupling coefficient /1 - |S|2 for the dominant TE'1 mode have been calculated, as presented in Fig. 4.35 (a) and (b). For validation purposes, the results for the reflection and coupling coefficients lsee Section 2.10

125 D.U.T. o -t / i - I I: - UT D.U.T.: Device Under Test ) Figure 4.33: An example of a CPW-to-waveguide probe station.

126 Top View A -- A a w waveguide side wall 1 W, = }W2 W3 It b Side View waveguide back short r L m 1! h L I -.7-777M 7, Er -k i h3 h2. I I. I Figure 4.34: CPW-to-waveguide transition: a = 47.6 mm, b = 22.0 mm, L1 = 10.8 mm, L2 = 0.8 mm, W1 = 0.5 mm, W2 = 0.8 mm, W3 = 0.5 mm, W4 = 11.8 mm, W5 = 3.8 mm, hl = h2 = 2.0 mm, h3 = 1.2 mm, h4 = 13.2 mm, and c, -= 12.0.

127 0 -10 m PC v lE 'a Cu.20 -30 3 4 5 6 7 Frequency [GHz] (a) 0 4 es a la / o: ^S1S1 FEM I. I..... I I I I I I -5s -10 3 4 5 6 7 Frequency [GHz] (b) Figure 4.35: Scattering parameters of the CPW-to-waveguide transition.

128 have been compared with the FDTD result, revealing excellent agreement between them. For the geometry analyzed in this study, the 10 dB bandwidth is very small (< 1 GHz) which is unacceptable as a wide band probe station. However, further study [28] reveals that a probe having material with lower dielectric constant, such as a membrane-supported structure, can provide much better performances. Also, for a given dielectric material, a wider patch (near square) exhibits improved performances. 4.10 Numerical Convergence The convergence characteristics of the finite element method largely depend on the electrical parameters under study. Specifically, it has been found that more than 30 cells per guided wavelength are required for the accurate evaluation of a via inductance as indicated in Fig. 4.36, while fewer cells are needed to compute electric fields and resulting scattering parameters. This is due to the form of the utilized basis function, which is linear and has compact support, and the approach followed herein for the computation of the inductance based on the flux/current definition. The electric field edge basis functions used in the study impose a linear variation for the electric field within each tetrahedron, but at the same time, provide a lower order variation for the magnetic field, which remains constant within the same volume Qe. The use of higher order basis functions, while it could complicate computations, may provide faster convergence for both the electric and magnetic fields. In parallel with checking the convergence of a solution with respect to mesh density, we need to make sure that the iterative solver converges. Even though there are many ways to define the convergence criterion for the linear equation solver, in

129 300,..., 250 1,500 5,000 200 2,000 10,000 15,000 ', 2,000 ct 150 10 Number of Unknowns 100 50 0.............. I........ I....... 0 5 10 15 20 25 30 35 40 Number of cells / Guided wave-length Figure 4.36: An example of the convergence test of the FEM. The inductance of the microstrip ground pad shown in Fig.4.6 as a function of the number of cells per guided wavelength. this study, it is specified as the relative residual norm, L2 norm,2 to be less than 10-4 to 10-6. The reason for choosing this error criterion is the sensitivity of the Sparameters with respect to the maxima and minima positions of the voltage standing wave pattern. The scattering parameters are found from the standing wave pattern of the voltage or current waves, and as a result, accurate values for the positions of the extrema are important for stable results. 2Let V be a finite-dimensional vector space. For every real number p > 1, the function 11[ IIp defined by is norm. In particular, it is called L norm, when p =. is norm. In particular, it is called L2 norm, when p — 2.

130 4.11 Conclusions The finite element method in the frequency domain has been applied to characterize several commonly found discontinuity structures in microwave and millimeterwave circuits. Specifically, microstrip short-circuit, microstrip ground pad, CPW-tomicrostrip through-via transition, channelized CPW-to-microstrip transition, asymmetric 2-port via with plated-through-hole, double via hole, planar/non-planar CPW-to-slotline, and finally CPW-to-waveguide transition have been analyzed and their electrical performance has been studied. The derived results from the FEM agree very well with that of the FDTD. Even further, the field visualizations in the discontinuity regions provide deep insight into the electromagnetic phenomena and their coupling mechanism. It has been found that a via hole at the center of a microstrip line fails to provide a good DC connection at higher frequencies, resulting in substantial coupling between the two sections of the line. It has also been found that the equivalent inductance for the microstrip via hole and its ground pad is sensitive to the scattering parameters due to ambiguity in the choice of the reference plane. Furthermore, the electromagnetic field distributions around via holes confirm the fact that the current tends to flow through shortest path from the line to ground plane. Three different types of the CPW-to-microstrip transitions have been found to operate efficiently in the low frequency region for the vertical transition and over a wide frequency range for the channelized transition. In particular, the latter case reveals its frequency independence (up to 40 GHz). In the asymmetric 2-port via hole and double via hole structures, lumped equivalent circuits have been developed for various geometrical parameters to provide guidelines for design of high density

131 MIC and MMICs. As a feeding network of slot antennas, two types of CPW-toslotline transitions have also been characterized in a back-to-back fashion. Finally, a preliminary study of the scaled W-band waveguide probe structure with open boundary FEM has been presented.

CHAPTER V HERMETIC TRANSITIONS The finite element method in frequency domain technique is applied to characterize the microstrip hermetic transition geometries in an effort to investigate highfrequency effects, and the derived results are compared with that of the the finite difference in time domain technique. Measurements are also performed on these transitions and compared favorably with theories. In this chapter, two different types of transitions structures have been analyzed from 10 to 25 GHz and have been found to be limited in performance by higher return loss as frequency increases. It is shown that microstrip-through-CPW hermetic transitions in the shielded environment may suffer from parasitic waveguide modes which, however, can be eliminated with the use of vias at appropriate locations. The hermetic wall on top of the CPW section shows a relatively small (< 2 dB) effect on the original circuit performance. Similarly, the hermetic bead transition shows good performance at a lower frequency region while it degrades as frequency increases. This indicates the need for very careful characterization of transitions intended for use in microwave and millimeter-wave applications. 132

133 5.1 Introduction Hermetic packages are frequently utilized in microwave integrated circuits to provide protection from a hostile environment in addition to reduced electromagnetic interference and radiation leakage. Although hermetic packages are designed to isolate oe one or more parts of the circuitry and protect them from parasitic electromagnetic noise generated by other analog or digital components, they should not affect the functionality of the packaged circuit. Despite these widely acknowledged and desired attributes of packages, the state of the art in packaging is hindered by design inefficiencies and shortcomings. As mentioned in many studies [149],[150],[186]-[198], the electrical performance of hermetic packages is not considered during the first circuit design stage and, in many cases, it turns out that the presence of the package drastically disturbs circuit characteristics. Recently, more rigorous analysis of the effect of hermetic walls and packages has been reported [170],[199]-[209] and reveals a complex and vague design problem which is largely intensified by uncontrollable package interference effects on circuit performance [193]. Inevitably, the circuit in a package needs to be connected to the outer world by means of effective transitions, and the performance of those transition structures needs to be carefully studied to obtain desirable overall circuit response. In most applications, a hermetic wall, which allows access to the interior of a package, is made of a perfect electric conductor placed on top of a low dielectric constant material (ceramic). The package is placed on top of the circuit in such a way that the input and output lines are sandwiched between the ceramic and the substrate in a stripline type of configuration. The transition of the feeding line, from an outside microstrip to an inside stripline through a hermetic wall, locally alters the characteristics of the

134 line, thus leading to reflections and many times localized resonance effects. There have been several studies on the effect of hermetic packages on circuit performance based on the use of an appropriate equivalent circuit which replaces the package walls. The element values of this equivalent circuit are computed through a combination of measured data and ideal transmission line theory. A hermetic package with a coplanar waveguide type transition was used in [188] to achieve input and output port insertion loss of less than 0.25 dB. In [189], a microwave surface mount package was modeled as a combination of coplanar waveguide and microstrip line sections with via holes to simulate the walls and equalize ground potential. In this situation, the circuit components and transmission lines were replaced by electric equivalent circuits whose element values were extracted from measurements. For relatively high-frequency applications (up to 20 GHz) with low transmission line loss and high port-to-port isolation, a surface-mounted MIC package with plated through holes to a microstrip transition structure has been reported [190]. Finally, in [191], a 70 mil microstrip package for a high-frequency high-power circuit application was modeled by equivalent RLC circuit and the values of those equivalent circuit elements were determined by measurements augmented with computer calculations. In this study, the effects of hermetic walls on the characteristics of two different types of transitions are studied. The first type is an intra-chip transition whereas the second type is an inter-chip transition in which the packaged circuits and connecting lines are printed on different substrates. As a result the transition geometries of the second type include finite gaps between the chips that may cause increased interference or even resonances and may substantially affect performance. The main characteristic of the first type transition is the use of a CPW configuration at the

135 transition region. This choice is based on the observation that the characteristic impedance of the microstrip line is determined mostly by the substrate thickness and conductor width, whereas that of the CPW is mainly affected by the ratio of the center conductor width to gap width. This implies that a microstrip transition could be more sensitive to the surrounding environment than a transition of CPW type. In the following sections, a 50 Q to 50 Q microstrip-to-CPW transition is studied and the effect of the hermetic wall which is laid on top of the CPW section is investigated. As another type of hermetic transition, a hermetic bead transition structure is also analyzed. Even though the equivalent circuit approaches based on measurement data can provide an understanding of the effect of the hermetic packages on the original circuit, rigorous theoretical investigation is required for the thorough understanding of the electrical performance of these transitions. In this work, the finite element method with edge-based vector basis function is used to characterize various types of hermetic transitions and the finite difference time domain method as a validation tool. The use of both techniques indicates that there is very good comparison with the exception of moderate differences at frequencies close to resonances. 5.2 Intra-chip Transition: Microstrip-through-CPW Hermetic Transition For an effective transition through a cavity wall, the microstrip line can be gradually changed into a coplanar waveguide which extends beyond each face of the hermetic wall as shown in Figure 5.1. The employed section of the coplanar waveguide is back metallized by the ground of the microstrip line and gives rise to parallel plate waveguide modes, parasitic radiation, and unwanted resonances that occur within

136 Microstrip CPW section section t II 4W Y _____ — ______ W5 Microstrip Section CPW Section W1 D i. W2 H~HW2 t r X l. i8b. w > Figure 5.1: Geometry of the Microstrip-to-CPW transition without hermetic wall. W = 0.48 mm, W2 = 0.14 mm, W = 02 mm, W4 = 2.0 mm, W5 = 3.0 mm, H = 0.635 mm, r = 12.5. Case 1: D = 2.26 mm. Case 2: D = 1.76 mm. the range of the operating frequencies. These resonances are pronounced when the structures are shielded, but are considerably damped when these transitions operate in an open environment. To understand these effects, the hermetic transition is studied in several steps. At first, a simple microstrip-through-CPW transition, operating in open environment, is analyzed using FDTD and compared with experiments to understand the parasitic effects such as radiation and resonances introduced by the CPW section. The transition structure is then placed inside a closed cavity and characterized by using the FEM and the FDTD. Because the side walls can affect the performance of

137 the circuit, one needs to be very careful in the choice of cavity width. In this study. to examine the effect of the cavity width, two different cavity sizes are considered without changing the other geometrical factors. In addition, vertical wires, rect — angular via holes, extended via holes, and PEC walls are employed to connect the CPW ground planes with the CPW backside metalization to eliminate the excited waveguide mode type resonances, and the effects of these connections on the original circuit are investigated. By doing this, the effects of the resonances due to the cavity and the waveguide-like section under the CPW section are closely examined. Finally, a hermetic wall is introduced to complete the transition geometry. 5.2.1 Open Transition Through a CPW Section To quantify the performance of the CPW section utilized in the transition, the microstrip-through-CPW planar transition in an open environment without via holes (refer to Fig. 5.1) is studied first. For the computation of scattering parameters using the finite element method, artificial absorbing layers are designed to simulate the open boundary problem as described in Section 2.10. The study in this section has revealed that the transition itself introduces a resonance around 17 GHz with a very low Q due to high parasitic radiation. Furthermore, the same structure is measured 1 experimentally using the 8510 ANA and the Wiltron test fixture. The measurements are performed using a TRL de-embedding procedure and time domain gating is employed to suppress multiple reflections due to the finite substrate. As shown in Fig. 5.2, the overall agreement between the FEM, FDTD and experimental results is is satisfactory. The discrepancy observed in the higher frequency end of the spectrum is due to the structural difference between the modeled and measured geometries and the performance of the absorbing boundary conditions utilized in 'The measurement was performed by Mr Eray Yasan in the University of Michigan

138 0 — Z- * -- - — ___-:_ -- -- 6 -— 7 -0 -4 -- / 200 3 o ---- ~ S1l: Meas. f o 0 E S21: Eeas. -30 - S1l:FDTD. --- —- - S21: FDTD -40 ---- S 1:FEM S21: FEM -50 10 15 20 25 Frequency [GHz] Figure 5.2: S-parameters for the open transition shown in Fig. 5.1. In this case, duroid substrate, which has a relative dielectric constant er of 10.8, is used. The CPW ground planes and the PEC underneath the substrate are connected in the measurements. both FEM and FDTD techniques. In particular, the FDTD approximated the exact geometrical factors of the transition to accommodate a uniform mesh and assumed an infinite substrate with the upper CPW ground and the lower conductor at exactly the same potential. However, in the experiment we had to use finite-sized substrate and, even further, the CPW ground planes and lower conductor were connected using conducting strips to provide the same potential. 5.2.2 Closed Transition Through a CPW Section Following the treatment of the open transition, the same geometry is placed inside a metallic cavity and is characterized for two different cavity sizes. Fig. 5.3 shows the computed S-parameters for D = 2.26 mm and 1.76 mm. As seen from the figures, the presence of the cavity, even though it is designed to be below cutoff at the frequencies of interest, drastically affects the performance by introducing high Q resonances. As

139 0. -10. -20. -30. -40. -50. 10. 15. 20. 25. Frequency [GHz] (a) 0 *o 3 ~a., c& rt 0. -10. -20. -30. -40. -50. -60. 10. 15. 20. 25. Frequency [GHz] (b) Figure 5.3: Comparison between FEM and FDTD data for (a) case 1 (D = 2.26 mm) and (b) case 2 (D = 1.76 mm).

140 expected [199]-[201], these resonances are shifted to higher frequencies as the cavity size decreases. For the case D = 2.26 mm, there are two resonances in IS211i which are clearly observed at 21 GHz and around 9 GlIz (not shown in the figure). These strong resonances are due to TE-like modes [210]-[211] excited in the region between the ground planes of the CPW section, its back metalization, and the side walls. In particular, the resonance around 10 GHz corresponds to the TE10 mode of that waveguide-like section as clearly shown in Fig. 5.7(a). The dips in S111 around 15 GHz in Fig. 5.3 are due to the length of the narrow microstrip line (refer to W5 in Fig. 5.1) as confirmed through another numerical experiment in Fig. 5.4. That is, when the effective length of that section is Ag/2, total transmission can occur with the same input and output transmission lines as in our cases. Also, similar argument can be applied to the drops at 20 and 21 GHz in Fig. 5.3 with the length of the CPW section. The return loss of case 1 is less than -10 dB in the frequency range of 11 to 20 GHz showing a good transition performance. Theoretical data derived by the FEM show very good agreement with that of the FDTD method and reveal resonance phenomena indicating the need for a mechanism to suppress them. The next subsections will be devoted to the investigation of variety of mechanisms needed to suppress the excitation of the parasitic waveguide modes responsible for the previously observed resonances. 5.2.3 Vertical Wire Grounding As a first attempt to eliminate the effect of the excited waveguide modes under the CPW section, vertical wires are inserted between the CPW ground planes and the PEC back-metalization under the substrate as shown in Fig. 5.5. On either side of the longitudinal symmetry plane, four wires (400,um apart) are inserted at a

141 O - 0ED- E-0-D e0 e D1-l-a e - — Eea I-, B 0o D-PP 0 0 -10 0 00, ^ 0 Q 0 Q ~ ~ ~ ~ 0 0 0 \ /0 Sll:W4=1.4mm o0 S21 W4 = 1.4 mm 4 e SI:W4 = 2.0 mm ------- S21:W4= 2.0 mm -50 —. 10 15 20 25 Frequency [GHz] Figure 5.4: S-parameters for two different size of CPW ground plane (W4 in Fig. 5.1). Other geometrical factors are remain the same. distance of 466 ytm from the symmetry plane. These vertical wires are implemented as infinitely thin PEC wires for convenient numerical simulations. The scattering parameters of the modified structure are computed using the finite element method and show very good agreement with the FDTD results as shown in Fig. 5.6. The results shown in the figure indicate that, in the low frequency region (up to 13 GHz), the transition circuit performs fairly well, while its performance degrades as the operating frequency increases beyond 13 GHz. This degradation indicates that the wires play a role in increasing the resonance frequency, but do not suppress it completely. This is due to a waveguide-type resonance (TEio-like mode) which develops between the wires and the side walls under the CPW section. They will always exist in some form as long as there is an enclosure. Fig. 5.7(b) shows the vertical electrical field distribution around the transition region at 10 GHz.

142 Top View (1/4 of the structure) I I I I ' ', 1 = "!I i: i. 1.: ) ~i i ii i ii.!. - i4i:i?i.iiii0i0iii-i~ii.iRi W1 W2 - - D1 D2 Figure 5.5: Geometry of the microstrip-through-CPW transition without hermetic wall: vertical wire grounding. W1 = 0.466 mm, W2 = 2.034 mm, D1 = 0.40 mm, and D2 = 0.40 mm. Other geometrical factors are the same as those in the previous example.

143 0 —. 0-o-,-.D-C.-G.. Dm * E1B OsE-8 —,c. D: -10. -30. - S: FEM S21: FEM -40. SI I:FDTD —... --- —- S21::FDTD -50. ~~,,I.. I.. 10. 15. 20. 25. Frequency [GHz] Figure 5.6: S-parameters of the vertical wire grounding case. 5.2.4 Rectangular Via Hole Grounding In this subsection, rectangular PEC via holes are used to connect the CPW ground planes to the bottom PEC surface as shown in Fig. 5.8. There are four vertical via holes, and the size of each via hole is 600 x 452,pm and the height is 635 um, which is the same as the substrate thickness. Also, the ground plane of the CPW is reduced to half of the original size to alleviate the strong resonance under the CPW grounds. The calculated S-parameters are shown in Fig. 5.9. The overall frequency response of the structure is similar to the wire grounding case, but the resonant frequency is shifted to a further higher frequency region. This is due to the reduced distance (W3 + W4 = 1.582 mm in Fig. 5.8) between the PEC via walls and the cavity side walls. Note that the aforementioned waveguide-like structure has incomplete side and top walls as shown in the figure, and this explains the reduced Q at the resonance of approximately 23 GHz.

144 (a) (b) - (c) Figure 5.7: Electric field (E,) distribution around the transition region at f=10 GHz. (a) Microstrip-through-CPW transition, (b) Microstrip-through-CPW transition with vertical grounding wires, and (c) Microstrip-through-CPW transition with extended rectangular vias.

145 A-A' Cut - - A|:::.; 7::::. j f Of: 0 m it. \: \ - Dy\. \ W1 ~ w2 W3 W4 7IFF"T- ]' --- / 0 —. 4-, - PEC Via Figure 5.8: Geometry of the microstrip-through-CPW transition without hermetic wall: Rectangular via hole grounding with finite size CPW ground plane. The size of the via hole is 0.452 x 0.6 mm. W1 = 0.226 mm, W2 = 0.452 mm, W3 = 0.452 mm, W4 = 1.13 mm. Other geometrical factors are identical to the first example.

146 0...,..-... *.. * -10.. - -20. - S21.FEM - Sll:FDTD ------- S21 'FDTD -30........ 10. 15. 20. 25. Frequency [GHz] Figure 5.9: S-parameters for the rectangular via hole grounding case. 5.2.5 Extended Rectangular Via Hole Grounding To avoid excitation of resonant waveguide-like modes under the conductor-backed CPW section, extended via holes are placed under the CPW ground planes as shown in Fig. 5.10 and their effect is investigated. Analysis of this structure confirms that the extended rectangular via holes suppress unwanted resonances (refer to Fig. 5.11). The field distribution around the transition region also verifies the suppression of the resonant field under the CPW ground pads. As shown in Fig. 5.11, the overall return loss is improved by about 2 dB by extending the via holes while the performance of the transition structure degrades as frequency increases. Specifically, over 15 GHz, it has greater than -10 dB return loss. The effect of the extended vias may be simulated by an array of vias.

147 Top View (1/4 of the structure) W7 W8 _ - i iii! li!i |__ W2 WE;B~& '^ W3 I Iiiii 1?i:? W 9 W4 t-1 W6 W5 Figure 5.10: Geometry of the microstrip-through-CPW transition without hermetic wall: extended rectangular via hole grounding. The top view shows one quarter of the whole geometry. W1 = 0.1 mm, W2 = 0.14 mm, W3 = 0.226 mm, W4 = 0.2 mm, W5 = 0.6 mm, W6 = 0.2 mm, W7 = 1.5 mm, W8 = 0.5 mm, W9 = 2.034 mm.

148 '. ' ' — *-4 — ** *- --....., ~4. v _...,, 4 -10..= o/' 00: SI I: FEM -20. X /* ~ S21: FEM S 11: FDTD -.- ----- S21: FDTD -30. ------------------ 10. 15. 20. 25. Frequency [GHz] Figure 5.11: S-parameters of the extended rectangular via hole grounding case. 5.2.6 PEC Wall Grounding In this subsection, the effect of further extension of the via holes to form solid walls under the CPW ground pads is studied. The schematic view is shown in Fig. 5.12 and Fig. 5.13 shows the computed frequency response of the transition with solid walls under the CPW ground pads. As expected, the transition characteristics are improved over the previous structure having extended via holes. After the PEC walls are inserted, the resonance phenomenon is completely vanished, and as a result the S111 is decreased. In particular, in the high frequency region, the circuit performance with the solid PEC wall is improved by 2 to 3 dB in return loss. Up to this point, the effects of the cavity size, vertical wires, and various shapes of via holes which connect the CPW ground to the PEC surface under the substrate have been examined. It is worth mentioning that the vertical via holes having various shapes can be implemented in practical circuits by removing the substrate material

149 A-A' Cut W1 Er PEC Walls ) Figure 5.12: Geometry of the microstrip-through-CPW transition without hermetic wall: PEC wall grounding. W1 = 2.034 mm and other geometrical factors are identical to the first example. followed by back metallization, so called micro-machining technique. In the next subsection, the effect of the hermetic wall on top of the circuit is studied. 5.2.7 Microstrip-through-CPW Transition With Hermetic Wall In this subsection, the effect of the hermetic wall on top of the microstrip-throughCPW transition having extended vias is investigated. The hermetic wall is formed by a metallic wall bonded on top of a ceramic or dielectric material which forms the upper part of the wall. Fig. 5.14 shows the geometrical details of the hermetic transition structure with extended via holes. The height of the ceramic is equal to the

150 0. 0B-B-9 00 0 - 0-D o 0 0 —B - 0 — 0 4 6- 0B-0-D B G -0-0-0 C 10. Sll:FEM -20. r S21: FEM S 11: FDTD ------ S21: FDTD -30..... E... I..i 10. 15. 20. 25. Frequency [GHz] Figure 5.13: S-parameters of the PEC wall grounding case. substrate thickness while its relative dielectric constant is chosen to be Cr2 = 2.3, in contrast to the 12.5 of the substrate. From a circuits point of view, the ceramic on top of the transmission line may alter the characteristic impedance of the line, and hence it may affect the resulting S-parameters. Since there is a conducting plane on top of the ceramic, the CPW with hermetic wall can be considered as an inhomogeneously sandwiched CPW section and its characteristic may be slightly affected. In Fig. 5.15, the magnitude of the scattering parameters obtained using the FEM and FDTD are compared. After placing the hermetic wall on top of the CPW section, the overall return loss has increased by 2 to 3 dB while its frequency dependency has remained unchanged. In the whole frequency region considered, from 10 to 25 GHz, the return loss is less than -5 dB. However, in the low end of this frequency region, up to 13 GHz, the return loss is less than -10 dB. The effect of the hermetic wall is small as expected from the field distribution in the CPW structure.

151 vias B-B' Cut AA' Cut Figure 5.14: Geometry of the microstrip-through-CPW transition with hermetic wall on top of the extended PECndding case. The width of the hermetic wa is 1.6 mm. WI = 0.48 mm, W2 = 0.14 mm, W3 = 0.2 mm, W4 = 1.6 mm, H1 = H2 = 0.635 mm, crl = 12.5, and C.2 = 2.3. H1 = H2 = 0.635 rnm, el = 12.5, and e,,2 = 2.3.

152 20. a- S21: -E}{-FEME o -10. '-- Sil:FDTD -30..... |-... 10. 15. 20. 25. Frequency [GHz] Figure 5.15: Comparison of the S-parameters for the transition with hermetic wall on top and extended PEC via under the circuit plane. 5.2.8 W-Band Application of the Hermetic Transition Based on the previous results, a transition from the microstrip-through-CPW section has been designed to investigate the possibility of using the hermetic wall transition in the W-band application. The hermetic wall is built on top of the CPW section of the microstrip-through-CPW transition, and the width and thickness are chosen to be 250 pm and 100 Fm, respectively, with a dielectric constant of 2.5. The S-parameters computed in the shielded environment using the FEM are shown in Fig. 5.16. It is observed that the transitions with and without the hermetic wall show very similar characteristics. In particular, in the high frequency region, over 70 GHz, the transition reveals very good performance, having less than -20 dB of return loss.

153 -0 -} — E 1- l-&0 W -h l -0l -I3-G - - -El-[D- D -3-0 -E-a -- -( - -20 6- J.0 0 o0 -30 - -40 E S: Without Wall ---- S21:Without Wall -50 0 SIl:With Wall 000 E S21: With Wall -60........ 50 60 70 80 90 100 110 120 Frequency [GHz] Figure 5.16: S-parameters for the transition with and without hermetic wall on top of the CPW section. The geometrical factors: W = 70 prm, W2 = 30 pIm, W3 = 40 prm, W4 = 300 /m, and W5 = 350 gm (refer to Figure 2). The width and height of the dielectric brick under the PEC wall are 250 fm and 100 pm, respectively.

154 5.3 Inter-Chip Transition: Hermetic Bead Transition Another type of commonly encountered feed-through, the hermetic bead transition structure, is analyzed in this section. Unlike the previous intra-chip transitions, the hermetic bead transition, so called, inter-chip transition, can be used as an effective interconnection method between chips and other circuit components through PEC wall. For optimum design of the transition structures, several different configurations, including square and circular dielectric rings, are carefully investigated in this study to provide design guidelines. 5.3.1 Squire Dielectric Ring As shown in Fig. 5.17, a circular coaxial line is approximated with rectangularshaped stripline and the dimensions are chosen to have 50 Q characteristic impedance. The dielectric constant of the ring is also chosen as 10.8, which is the same as the substrate. The thickness of the hermetic wall and the air gap spacings are chosen as 1.5 mm and 0.4 mm, respectively. S-parameters computed using the FEM are presented in Fig. 5.18 showing very good agreement with the FDTD results. The hermetic bead transition performs quite well in the low end of the frequency spectrum as expected, while the return loss degrades as frequency increases. To improve the performance of the transition and to provide useful design information, the air gap spacing (L2) is varied from 0.4 mm to 0.1 mm. The computed S-parameters are shown in Fig. 5.19 for three different cases. As can be observed, the overall transition performances are improved as the air gap spacing is decreased. Although the geometrical factors are not optimized for better transition characteristics, the overall performance of the circuit gives an essential understanding of the performance of this type of transition.

155 A IH1 Top View Cross-Sectional View (A-A') W5 l- -I H2 - I I ' W2 1. W2 W1 W4 Figure 5.17: Geometry of the hermetic bead transition. W1 = 0.55 mm, W2 = 0.21 mm, W3 = 1.27 mm, W4 = 2.225 mm, Ws = 5.0 mm, H1 = 0.635 mm, H2 = 4.0 mm, L1 = 1.50 mm, L2 = 0.40 mm, and Er = 10.8.

156 0 --20 S: FEM -0 S21:FEM -40 SI:FDTD.X... ----. S21: FDTD -50 10 15 20 25 Frequency [GHz] Figure 5.18: S-parameters of the hermetic bead transition with square dielectric ring. 5.3.2 Circular Dielectric Ring For an investigation of the effect of the rectangular approximation of a circular ring, the hermetic bead transition with circular dielectric ring is studied in this subsection. The radius of the circular ring is chosen as H1 so that the circular ring is fit into the square hole. As shown in Fig. 5.20, the performance of the circular ring is very similar to that of the square one except at the lower frequency region. The slightly degraded performance of the circular ring is mainly due to the effect of the conductor placed closer than the square one. 5.4 Conclusions The finite element method has been successfully applied to characterize the effect of the hermetic wall/bead transitions commonly found in MIC and MMIC packaging. The derived FEM results agree very well with the FDTD data in low as well as high frequency regions. It has been shown that the grounding structures (wires, vias, and

157 la C. 4 — Ca m 0 -10 -20 -30 -40 -50 10 15 20 Frequency [GHz] (a) 25 1 0 m "0 3 d, bD c3 -1 -2 -3 -4 -5 1. - I - - -- -1 1 — - -- S21: L2 =0.4 mm S21:L2=0.2mm S21:L2=0.Imm.... I.... I 10 15 20 25 Frequency [GHz] (b) Figure 5.19: Comparison of the S-parameters for the hermetic bead transition of three different air gap width L2. (a) |S11| and (b) IS211.

158 0 -- - —.- -- -.-. _ - -10 -20 / Square Circular __ Ring Ring 30 / S11l:Qularr g o/ S21: iculardWng;:.: -40 ---- S11:SquaeIg I - S21: Squarering -50 10 15 20 25 FIequency [GHz] Figure 5.20: Comparison of the S-parameters of the hermetic bead transition with square and circular dielectric ring of radius equals H1. PEC walls), which are placed under the CPW ground pads to equalize the potential between the ground planes on the upper and lower surfaces, play a key role in the performance of this type of transition. The hermetic wall placed on top of the CPW transmission line section degrades the return loss of the original circuit only by a couple of decibels. In addition, a hermetic bead transition structure was designed and analyzed for the circular- and rectangular-shaped rings, and the effect of air gap spacing is also studied. The performance of this transition is found to gradually degrade as frequency increased.

CHAPTER VI K/Ka-BAND MMIC PACKAGE In this chapter, electromagnetic leakage and spurious resonances in a K/Ka-band (18 - 40 GHz) MMIC hermetic package designed for a phase shifter chip are studied using the finite element method and the numerical simulation results are compared with the measurement data. Both the measured and calculated data indicate several spurious resonances in the 18 to 24 GHz region, and the origin of this phenomenon is identified by virtue of the modeling capability of the FEM. Moreover, the effect of I)C bias lines, bond wires, shielding and the asymmetry of the package on electrical performance are closely examined. In addition, the effect of adding a resistive coating to the inside surface of the package lid and also the use of dielectric packaging materials with very high loss tangent are studied in view of the suppression of the spurious resonances. Finally, design guidelines for the improved package are presented. 6.1 Introduction In general, packages for low frequency integrated circuits (ICs) have been designed without consideration of the high frequency effects and have revealed satisfactory performance. However, high performance packages, especially for microwave and millimeter-wave integrated circuit applications, should satisfy stringent mechan 159

160 ical, electrical and environmental requirements. From a mechanical and environmental point of view, a package should provide protection to the internal circuits from the surroundings. Furthermore, packages are required to exhibit minimum insertion loss, good isolation between the ports as well as electromagnetic shielding for minimum interference [193, 195]. Another important electrical requirement of a package is non-invasiveness with respect to circuit performance [212]: a package fails electrically if parasitic cavity resonances substantially deteriorate circuit performance. As a result of all of these requirements, successful development of a high frequency package requires careful design strategies. Recently, low cost high performance MMIC packages have been developed by using approximate equivalent circuit models or experimental data [186, 188]. However, due to the limited accuracy of the modeling tools, the designed packages exhibit a serious performance degradation at higher frequencies. To overcome the above difficulties in designing high-frequency high-performance package, frequency domain full-wave electromagnetic tool, FEM, is applied [31, 208, 209, 213] for various applications. The goal of this study is to characterize the electrical performance of a K/Kaband hermetic package designed for a MMIC phase shifter (refer to Fig. 6.1) and comprehend the parasitic effects introduced by the package geometry. For thorough understanding of the package performances, the effects of the various features of the package, such as metal filled vias, I)C bias lines, bond-wires, structural symmetry/asymmetries, and even the effect of the test fixture on the circuit performance are extensively investigated. Furthermore, the use of dielectric packaging materials with high loss tangent and the influence of a resistive coating on the inside surface of the package lid are also examined. For the EM characterization of the K/Ka-. band MMIC package, a parallelized three dimensional finite element method, which

161 1 mm Figure 6.1: Top view of the K/Ka-band MMIC package. is optimized on the distributed memory machine with task parallelization strategy, is applied [20, 29]. The parallelized 31) FEM code exhibits near linearly scalable performance improvement due to the frequency independent nature of the frequency domain FEM. To the best of the author's knowledge, this is the first comprehensive high frequency full-wave treatment for an existing millimeter-wave package. The modeling effort described herein is divided into three parts. In the first part, the input and the output microstrip lines are symmetrically located with respect to the two planes of symmetry of the package. In addition, the package is considered to be free standing, i.e., totally isolated. Under these assumptions only one quarter of the package needs to be considered, and as a result the modeling effort is conlputationally simpler. In the second case, the input and the output microstrip lines are asymmetrically located with respect to one of the plane of symmetry. This is of

162 more general interest since in real applications the transmission lines on the MMIC chip need not be symmetrically located. In the last part, the microstrip lines are asymmetrically located and the package is placed inside a test fixture. By placing the package inside the text fixture, the interactions between the fields leaking through the package dielectric walls and the surroundings can be modeled. This model is of great practical interest since it points towards ways and means to improve future package performance. The modeled characteristics are compared with experimental results and show very good agreement. In addition, the effects of the three different types of packaging features, such as additional metal-filled vias, resistive coating, and lossy packaging materials, are carefully examined to improve the package performance. This modeling procedure has the potential to predict the performance of other types of packages such as those used in wireless communications which include multi-chip modules. 6.2 Modeling of the Package A schematic is shown in Fig. 6.2 for the K/Ka-band MMIC package fabricated by Hughes Aircraft Company for the NASA Lewis Research Center. The package has 50 Q microstrip input/output feed lines for the RF signal and two sets of five metal lines on either sides for DC bias and control. In addition, a set of twelve metal filled vias tie the top and bottom PEC ground planes to provide mechanical strength and also serve as an EM shield. The package is fabricated from alumina (92 % pure, ~r = 9.5) using the high temperature co-fired ceramic (HTCC) process. To characterize the package, a 50 Q through line is placed in the recess as shown in Fig. 6.2 and is wire-bonded to the input/output microstrip feed lines. The peripheral dimensions of the package are 7.112 x 7.112 x 0.762 mm. It can be observed from

163 Ground Plane Figure 6.2: Schematic diagram of the K/Ka-band hermetic package designed and manufactured for a MMIC phase shifter chip. The diameter of the vertical metal filled vias are D, = 0.203 mm and the distance between these vias are Ds = 1.016 mm. The upper and lower alumina layers are each 0.381 mm thick (Dt).

164 the figure that the metal filled vias and the input/output microstrip feed lines are displaced towards one side of the package introducing a package asymmetry. This asymmetry is attributed to the specific geometry of the MMIC phase shifter which the package intends to house. While this package provides hermeticity, it does not provide perfect electromagnetic shielding, and as a result the packaged circuits are exposed to a semi-open environment. For the simulation of this environment, artificial absorbing layers have been designed using lossy isotropic dielectrics as shown in Fig. 6.3. The performance of the absorber is controlled by assigning a certain amount of losses in the absorbing material and by specifying its thickness 1. In this study, we designed two different types of isotropic absorbing materials for air and dielectric side terminations. The material parameters are chosen to be 6 =a = iJ = 1.0 + j10.0 for the air side, and diA = 9.5 + j15.0 and dil 1.0 + jl.5789 for the dielectric side. The thickness t of the artificial absorbers is assigned to be 0.70 mm to allow enough damping of the fields inside of the absorbers and to minimize the computational domain. It is well known that this type of absorber performs well for near-normal incident fields. 6.3 Numerical and Experimental Results and Discussions 6.3.1 Modeling of Isolated Symmetric Package First of all, an isolated symmetric package is designed and modeled, in which the structural details are very close to that of Fig. 6.2 except a few minor modifications for symmetric arrangement. The input/output microstrip lines and the through line are adjusted to the center of the package, and the two asymmetric metal filled vias are aligned in line. The computed scattering parameters are shown in Fig. 6.4. The return loss is less than -20 dB over the entire frequency range and the insertion 1refer to Section 2.10

165 Artificial absorber Top of the package air t ldiel l l. Coaxial connector Hf height for RF in/output PEC wall PEC wall Full height (recess depth) PEC wall (test fixture) Figure 6.3: Schematic showing the location and extent of PECs in addition to the artificial absorbers which simulate the test fixture and open environment. Two different types of absorbers (iAA = A = 1.0 + j10.0, diel 95 + j15 5.0, and pl 1.0 + jl.5789) are designed and placed in either side of the package facing the DC bias lines (not fully shown in the figure for simplicity).

166 0 r- -0-D-a -E —D- -E{- -- -El -0 B B- -B-[e-0- Q-E} -(-a-a mI0 SI II: -—.-.-. I S211 0 -20 - H -30 -40 -50...,...,...,I...,......I...a...,..... 18 20 22 24 26 28 30 32 34 36 38 40 Frequency [GHz] Figure 6.4: Computed S-parameters for the isolated symmetric package. loss remains within -1.0 dB. As will be shown in the later sections, the symmetric arrangement of the package impacts on the performance of the package greatly. 6.3.2 Modeling of Isolated Asymmetric Package In this section,, the performance of the isolated asymmetric package is modeled and the computed S-parameters are presented in Fig. 6.5. As can be observed in the figure, the package reveals relatively good performance with better than -10.0 dB return loss over the entire frequency range. It is also noted that no cavity or spurious resonances are observed, even though the size of the overall package becomes larger than half guided wavelength at frequencies above 15 GHz. The absence of any resonance is also evident from Fig. 6.6, which shows the computed vertical electric field distribution in the package, and is attributed to the imperfect EM shielding. Even though the package provides excellent mechanical/environmental protection and hermeticity, the virtual side walls formed by the 12 vertical metal filled vias

167 0 -— E) — -& ---El ----[3 —o-0-l-t- — ~-[-Q-3-o-I-~ — -- -— Q -— ~-E- - -10., -20 -30 -40 — -... I S21I -50.... 18 20 22 24 26 28 30 32 34 36 38 40 Frequency [GHz] Figure 6.5: Computed S-parameters for the isolated asymmetric package. between the DC bias lines do not provide a solid EM shield. Leakage from the cavity semi-open walls drains the energy from the cavity and in turn contributes for the suppression of the occurrence of cavity resonances. The effects of the bond wires and DC bias lines on the characteristics of the package are also studied. This is accomplished by computing the scattering parameters of the modified package: the four bond wires are eliminated between input/output RF microstrip lines and the 50 Q through line remaining 40,um gap between them, and 5 DC bias lines are also removed from both sides of the package. The computed scattering parameters for the simplified structure shown in Fig. 6.7 reveal no discernible differences from the original package. This result implies that the presence of the bond wires and the DC bias lines is not critical to the performance of the package. Furthermore, it can be argued that the shape of those structures, namely bonding wires and DC bias lines, should not disturb the overall characteristics of the

168 -lo 1 110Xlx-ZO 13~~:::1130 (a) 11 EII E As-50 -80 (lb) (b) f-= 29.5 GHz in the isolated asymmetric hermetic package.

169 IS 111 - original -------- IS211 - original o IS Ill- simplified I IS211 - simplified, I I I I I I I I - o. ----— EEr — ----- r — ---- -—. --- -- --- - -10. -20 w Z-30 -40 -50... 18 20 22 24 26 28 30 32 34 36 38 40 Frequency [GHz] Figure 6.7: Comparison between the computed S-parameters for the isolated asymmetric package with (original) and without (simplified) DC bias lines and bond wires. package. Compared to the performance of the symmetric package shown in Fig. 6.4, the return loss of the asymmetric package surprisingly increases about 10 dB. 6.3.3 Modeling of Asymmetric Package Placed in a Test Fixture The asymmetric package is now placed inside of a test fixture and the whole structure is modeled to quantify the susceptivity of the package to the environment and to identify the electromagnetic field leakages and spurious resonances. The test fixture is modeled by using PEC walls along the ends of the package with an opening at the center for the input/output coaxial connectors as illustrated in Fig. 6.3. The

170 half height PEC walls are designed for the simulation of the recess depth of the test fixture, and the absorbers are placed around the package on top of the half PEC walls. The measured and computed scattering parameters for the package placed in the test fixture are presented in Fig. 6.8 revealing very good agreement. As can be observed, the overall structure including the package and test fixture suffers from spurious resonances in the low frequency region (18 to 24 GHz), even though it exhibits very good performance in the rest of the frequency range. The resonance phenomenon can be also understood by investigating the EM field distributions at various frequencies as shown in Fig. 6.9. For example, at f = 20.5 and 22.0 GHz, ISn11 has peaks (high insertion loss) indicating strong mismatch and energy leakages through the input/output RF microstrip lines as well as between the metal filled vias as clearly shown in Fig. 6.9 (a). It is also observed that at these frequencies the excited EM fields are strongly concentrated in the package frame, having the shape of a dielectric ring, placed on top of the input/output microstrip lines and the DC bias lines. Even though the strong resonances in the package should be suppressed by allowing energy leakage, it should be noted, however, that the poor EM hermeticity could cause a serious EM compatibility (EMC) problem. Fig. 6.9 (b) shows the field distribution at 29.5 GHz, where IS111 has a dip indicating small insertion loss. The calculated results also show the voltage standing wave pattern of the EM field on the microstrip line. 6.4 Features Which Can Enhance Package Performance 6.4.1 Asymmetric Package with Additional Metal Filled Vias To investigate the role of the vertical metal filled vias on the EM hermeticity of the package, and to study the energy leakage from the input/output ends of the

171 r- " +-j W) C7 -10 -20 -' \: -30 - -40 -50 -—.-. 18 20 22 I \ -1 24 26 28 30 32 34 36 38 40 Frequency [GHz] (a) - II I -2 - I -4- I -6 -8...-.... I 18 20 22 24 26 28 30 32 34 36 38 40 Frequency [GHz] (b) Figure 6.8: Measured and computed S-parameters ( (a) IS11 and (b) IS211 ) for the asymmetric package residing in the test fixture.

172 30 -50 -0 -... - O............ - 3 0 -20 -70 -40 -SO -.1.............;.'''' -60 -70 (b) metic package placed in the test fixture a~t (a) f = 20.5 GHz and (b) f =29.5 GHz in the isolated asymmetric hermetic package. 1|-20 | #0'IIti'.. -..........i 1 -40 GAnth slte sm etrica hemti ac0e

173 50 n / T Thru Line Additional fillvMetal-filled vias RF In/Out "4* T^ Ground Plane Figure 6.10: Schematic diagram of the asymmetric hermetic package showing four additional vias at the input and output ends of the package. package, four additional vias are placed as shown in Fig. 6.10. The dimensions of the additional vias are the same with those of the previous set of vias. As shown in Fig. 6.11, the frequency response of the package is quite degraded in the most of the frequency region, indicating the fact that providing additional vias could aggravate the package internal resonances by suppressing the leakage of the EM fields. Interestingly enough, the enforcement of strong EM hermeticity by placing additional vias causes strong unwanted internal resonances, while poor EM hermeticity makes the package more susceptible to the surrounding environment. To overcome the above contradictory difficulties, the effect of using lossy packaging materials are examined in the following sections for the suppression of the unwanted cavity resonances as well as energy leakages.

174 0 -5 r- " a) 73 ~es. T-.-10 -15 -20 -25 — 30 — 35 18 20 22 24 26 28 30 32 34 36 38 40 Frequency [GHz] Figure 6.11: Computed S-parameters for the asymmetric package with four additional vias at the input and output ends of the package.

175 6.4.2 Asymmetric Package with Resistive Coated Lid As an effort to suppress the unwanted cavity resonances in the 18 to 24 GHz region (refer to Fig. 6.8(b)), a resistive sheet backed by PEC plane is placed on top of the package. The value of the resistivity R of the sheet is chosen to be 100 Q/o and 400 Q//, and the thickness t to be 0.1 or 0.3 mm. As shown in Fig. 6.12, the spurious resonances in the low frequency region (18 to 24 GHz) are completely suppressed for all 3 cases, and there are no distinguishable differences between them. This observation indicates that placing resistive material on the inside surface of the package lid can suppress the unwanted resonances by imposing lossy boundary condition on the fields resonating inside the package. 6.4.3 Asymmetric Package with Lossy Dielectric Frame It is very important to take measures for the suppression of the internal resonances during the first design stage in order to avoid performance degradation. Suppression of the internal resonances may be accomplished through a variety of approaches. One approach may suggest the design of a leaky package so that strong resonances are not supported. This approach has been examined in the previous sections and reveals that the package interferes strongly with its surroundings. As an alternate approach, one can fabricate the package frame from a dielectric material having nonzero loss tangent. Fig. 6.13 shows S-parameters of the package with two different dielectric materials: case 1: r = 9.5(1.0 + j0.1), /r = 1.0 + jO.0, and Case 2: = 9.5(1.0 + jO.5), Lr = 1.0 + jO.O. As can be realized from the figure, the unwanted cavity resonances in the low frequency region are completely suppressed by the assigned small amount of loss in the packaging material. It is also important to notice that the insertion loss increases to -2 dB due to the loss in the packaging

176 " -0 6 —.. Vr — r —.4 0 -10 -20 -30 -40 -50 18 20 22 24 26 28 30 32 34 36 38 40 0 -1[ ": C'l -n cV Frequency [GHz] (a) R= 100, t = 0.3 --- R=400, t = 0.3 R=400 t=0.1... I... I... I... I... I... I... I... I... I I...I... -2' -3 -4 - N 18 20 22 24 26 28 30 32 34 36 38 40 Frequency [GHz] (b) Figure 6.12: Computed S-parameters ( (a) IS111 and (b) IS211 ) for the asymmetric package residing in the test fixture with resistive coating on the inside surface of the lid. The resistivity and the thickness of the coating is indicated as R in [9/o] and t in [mm], respectively.

177 0 U~ e —e-6 *~ e@0-4 b' b - '- *- b -* * * * * * l -10 / o\ () -20 - IS211:casel '-0 -30 o -40 ISllI:case2 I S21 I: case2 -50,,,...,I I I I 18 20 22 24 26 28 30 32 34 36 38 40 Frequency [GHz] Figure 6.13: Computed S-parameters for the asymmetric package with frame constructed from dielectric material with high loss tangent. Case 1: Er = 9.5(1.0 + j0.1) and /r = 1.0 + j0.0. Case 2: Ec = 9.5(1.0 + j0.5) and Ir = 1.0 + jO.0. material. The differences in the scattering parameters for two cases corresponding to different loss tangents remain negligible in the whole frequency spectrum. 6.5 Conclusions In this study, a 18 to 40 GHz hermetic package is experimentally characterized and also modeled using the finite element method. The FEM accurately predicts the S-parameters, energy leakage and spurious resonances which degrade the package performance. To the best of author's knowledge this is the first comprehensive study of RF leakage and resonances in a millimeter-wave package. The modeled results show that the unwanted spurious resonances in the package can be suppressed by incorporating a resistive coating on the lid and by the use of dielectric materials with high loss tangent. The modeled S- parameters for a symmetric package predict low insertion loss on the order of -1.0 dB over the 18 to 40 GHz band. From this

178 study, it becomes clear that symmetry in the package construction is crucial for good performance.

CHAPTER VII SYSTEM LEVEL EM MODELING OF DIGITAL ICs AND PACKAGES In this chapter, as an application of the finite element method to system level electromagnetic modeling for digital ICs and packages, a new hybrid technique is developed combining a full-wave EM field analysis tool (FEM) and a circuit simulator (HSPICE). The three-dimensional FEM has been used as a tool to extract equivalent circuits from the detailed structures and the circuit simulator is the used for efficient time domain analysis of the lumped element network. This hybrid methodology provides accurate noise maps' at any place in the circuit topology and guides effective placement of decoupling capacitors. 7.1 Introduction Inductively induced voltage fluctuations in the power and ground lines of digital systems are a source of performance degradation and may pose serious reliability problems. Variously termed simultaneous switching noise (SSN), inductive noise, delta-I noise, and ground bounce, this phenomenon is most severe when a large nunber of high frequency chip drivers switch simultaneously, and cause a large amount 1The noise maps are generated by measuring maximum, minimum, or peak-to-peak voltage fluctuations on any given plane, such as the power and ground planes. In other words, at each electrical node, the desired values are recorded in the whole simulation period. As a result, the values in the noise maps, in general, correspond to different time instants. 179

180 of current to be injected into the power distribution network. The adverse effects of SSN are manifested in a variety of transient and permanent system malfunctions, including the appearance of undesirable glitches on what should otherwise be quiet signal lines and the flipping of state bits in registers and memories. Because of their unpredictability, SSN-related bugs have been reported to be among the most difficult to track down and correct [214]. While the fundamental mechanism underlying SSN, namely current switching, has been widely recognized for many years [215]-[219], the computer-aided design tool kits available to a digital system designer today still lack adequate SSN modeling, analysis, and design capabilities. This lack of modeling tools can be partially attributed to the complexity of the EM modeling required to obtain accurate estimates of the inductive parasitics responsible for SSN, coupled with the need for system-wide global analysis of this effect. Current approaches for dealing with SSN are largely ad hoc, relying primarily on the ability of expert designers to postulate worst-case scenarios for the occurrence of SSN-related errors and to analyze these scenarios using pessimistic estimates of packaging parasitics. The goal of this study is to take a first step towards evolving a systematic methodology for modeling and analysis of SSN in printed circuit boards (PCBs). In this methodology, the seemingly contradictory goals of modeling accuracy and global analysis efficiency are reconciled through a divide-and-conquer process, namely micro- and macro-modeling. Accuracy is insured by performing detailed EM field analysis on appropriately chosen small sections of the PCB, referred to as "tiles", to create lumped electrical equivalent circuit models. PCB tile models are subsequently combined with models for chip current drivers and package leads to produce an electrical simulation model for the PCB power distribution subsystem. A circuit

181 simulator such as SPICE [220] is then used to exercise this model under a variety of current excitation conditions to yield noise maps that indicate the variation in power and ground potential as a function of location on the PCB. It is instructive to note that this approach is similar, at least in spirit, to the approach routinely employed in modeling and simulation of very large scale integrated (VLSI) chips. Accurate circuit simulation of VLSI chips containing tens of thousands of transistors is made possible by separating the detailed physical models of 3-D transistor structures from the global chip-wide electrical analysis. Device modeling involves creating lumped electrical equivalent circuits that are obtained by solving the partial differential semiconductor equations governing the movement of charge carriers within the semiconducting material. These equivalent circuits, appropriately parameterized by material, geometrical and environmental factors, are then combined to produce a complete electrical model for a chip that is amenable to analysis by fast circuit simulators. The rest of this chapter is organized in four sections. Section 7.2 reviews current approaches to the modeling and analysis of SSN and motivates our approach to the problem. The tiling procedure central to our method is described in Section 7.3 along with the equivalent circuits for three different tile types. In Section 7.4, several simulation results are presented for a number of different PCB configurations that are aimed at enhancing the understanding of the SSN phenomena. In particular, the hybrid approach developed in this study is validated by comparing the measurement data provided by Intel for their P6 "gadget". The contributions of this work and directions for future work are summarized in Section 7.5.

182 7.2 Problem Statement and Previous Work A simplified model of the current path through a PCB-mounted IC chip is shown in Fig. 7.1 [215]-[217],[221]. Current enters the PCB through a metal connector from an external power supply (denoted as VDD) and flows through the conducting material of the board's power plane, the chip's package leads, and bonding wires to a VDD pad (point A) connected to the chip's internal power distribution network. Current then flows through one or more output pad drivers C and returns through a Vss pad (point B), bonding wires and package leads, to the board's ground plane and ground connector to the external power supply. In this simplified model, lumped inductances are used to represente the PCB power and ground planes (labeled LPCB) and the package leads (LpIN). It is also assumed that the internal chip power distribution network is mostly resistive and does not significantly contribute to the overall inductance responsible for SSN. Techniques for efficient design of chip power distribution are discussed in [222]. When a pad driver C changes state to charge or discharge its load capacitance, the resulting current surge through the power distribution network produces inductive voltage drops and in turn lowers the potential at point A below VDD and raises the potential at point B above Vss. Such glitching, or noise, in what should otherwise be quiescent power and ground rails leads to harmful effects on the chip circuitry, such as increased driver delay and loss of signal integrity [215, 216]. The increase in driver delay is directly attributable to the temporary lowering (raising) of potential at point A (point B) causing a reduction in the current drive available for charging (discharging) the load capacitance to VDD (Vss). The second effect is more insidious, as it may lead to logic errors. Since the power and ground lines establish the reference

183 VDD LPCB ------ --------— B PIN D C Chip l l LPIN.- --------- - -- - - -- - -- - - - - - - -I LPCB Vss Figure 7.1: Traditional way of modeling the inductances of the PCBs and package bonding wires. levels for all other signals, fluctuations in these levels may be misinterpreted as logic transitions in data and clock lines. If these transitions are inadvertently latched, the circuit may end up in an incorrect logic state. For example, a dip in potential at point A may be seen as a spurious glitch at the output of the quiet driver D that flips the state of latch E. These problems are compounded when a large number of output drivers switch simultaneously on a PCB with many chips. Much of the previous work on SSN can be characterized as seeking (1) the determination of an aggregate effective inductance Leff that accounts for the power and ground leads, the board itself, and the mutual inductances between the chip packaging and the PCB, and (2) the P anhedevelopment of design guidelines to reduce SSN effects. Some of the early work in this area was done in [216], where a detailed electrical model of the package is simulated using ASTAP [223] and a gross value of

184 Leff is determined using: where Vn is the noise magnitude, N is the number of simultaneously switching drivers di and is the rate of change of driver current. II1 this approach, only a single chip is considered and the computed Leff is used to represent the noise effects. In a similar manner, Rainal [218] gives detailed formulas for calculating the self and mutual inductances, assuming that the pin inductances are the principal sources of switching noise on a chip. Using detailed electrical models of the package and the noise tolerance level of a logic family, Katopis [217] applies statistical design rules to determine the maximum number of drivers that may switch simultaneously. Furthermore, EDN's Advanced CMOS Logic ground bounce tests [224] conclude that centrally placed ground pins or surface mounted packages reduce the ground bounce effects. In all these treatments, the parasitic inductance of the PCB is ignored and the lead inductances of the package pins are considered to be the primary contributors to SSN. Recently, however, several attempts have been made to account for the effect of the PCB power and ground plane inductances. Specifically, the current distribution on the power and ground planes is computed by solving the field equations using source and sink points. In [225], the switching noise due to inductances in the board as well as the vias is computed using an appropriate SPICE model. The board is assumed to have four layers: signal, power, ground, and reference layer for the measurements. The model is created by appropriately discretizing the power and ground planes and finding the equivalent circuit for each mesh. More recent work [219] has shown that the induced noise voltage is, in fact, sub-linear in the number of switching drivers. This phenomenon is due to the fact

185 that the ground bounce acts as a form of negative feedback, reducing the drive of the CMOS transistors. In this work Leff is calculated by placing centralized point sources and distributed sinks on the power planes. The static field equations are then solved for the power planes, and the Leff is computed in a modified form of Eq. (7.1) to take into account the negative feedback effect. The level 1 SPICE model of the CMOS transistors used in deriving the noise equations limits the accuracy of the model when it is applied to sub-micron technologies. Moreover, the computed value of Leff varies with location when many chips are placed on a PCB. A detailed EM analysis of the interconnect structure in a computer package is performed in [170, 226] using the FDTD and transmission line matrix (TLM) methods. While accurate, these approaches become computationally infeasible for realistic systems where a PCB may contain several hundred ICs with multiple layers of power, signal, and ground planes. An accurate SSN simulator must consider the inductances due to both the PCB and the pins themselves. The technique developed during the study presented in this chapter uses a combination of both macroscopic and microscopic modeling approaches, as illustrated in Fig. 7.2. The micromodeling provides an accurate electrical representation of the board as an interconnection of N-port networks. These networks are three-dimensional ladders of RLC elements and represent the electrical equivalents of subsections of the board, referred to as tiles. The tile equivalent circuits are defined according to the expected electric field signatures inside the board, and high accuracy is achieved by taking into account the dynamic nature of the electromagnetic fields. Depending on the substructure, these models may be based on quasi-static or full-wave solutions. The resulting board equivalent circuit is augmented with models for chip current drivers and simulated using HSPICE [227].

186 PCB Config EM-Based -' Tile Equiv Ckt [ ~lj ] Library Ii UI I f z Tiling Processor ) U PCB Lumped Electrical Ckt Model." JC 0 Co e'a 0 C) 0.-" rA C) 'e C.. "5L C. Current Source MA A Model 4 I Circuit Simulator ) Voltage Noise Map r 40 r - I -T'" I - | 5 tn7 40 10 T — Figure 7.2: Modeling of simultaneous switching noise (SSN) and simulation methodology.

187 The macroscopic model considers the SSN at the system level and accounts for its variations across the PCB due to the flow of current along the conducting paths on all planes. A noise map of the power and ground planes is obtained as a result of the macromodeling and the regions where the effect of noise is maximum are identified for the possible insertion of bypass capacitors. While the hybrid approach developed in this study is similar to [225], it has the advantage of combining the accuracy of a full-wave field simulator with the large-scale capacity of SPICE electrical simulation, and can be easily generalized to handle any number of power, ground and signal layers. Moreover, since the equivalent circuits resulting from the full-wave analysis consist of lumped R, L and C elements, we can potentially use the AWE-based simulators, which are 3 to 4 times faster than SPICE [228], provided that the drivers are linearized appropriately. In the following sections, detailed descriptions of the micromodeling approach will be presented, while relatively short discussion is devoted to the macromodeling in view of its well-founded feature. 7.3 PCB Tiles and their Lumped Equivalent Circuits In this section., the hybrid modeling methodology is presented starting with the tiling procedure and concluding with a summary of the primitive equivalent circuits structures which are the main building blocks of the printed circuit board geometry. 7.3.1 Tiling Procedure The board is partitioned into a set of tiles, primitive structures, using a systematic discretization, or so-called tiling, procedure. In this procedure, grid points (nodes) are introduced for all the physical details of the structure. These physical nodes correspond to crossings and bends of traces, via connections to planes, corners of

188 the printed circuit boards, connections to connectors and drivers, apertures on the planes, connections to loads, decoupling capacitors, etc. These nodes exist on all planes of the PCB and may vary in number from plane to plane. The intent of the tiling approach is to construct a three-dimensional grid that preserves all the geometrical details of the printed circuit board. The resulting tiles are rectangular and non-uniform, and their grid points provide the electrical nodes for the threedimensional SPICE model which is generated at the end of the tiling process. It is important to notice that the interactions between the neighboring tiles are fully accounted for by electrically connecting them together. Any spurious coupling between them, which is a second order effect, however, is not considered here, unless the tile equivalent circuit reflects it. As a result, to minimize the error due to the tiling approach, all the closely coupled geometries have to be modeled together to yield an appropriate tile and this can be accomplished by applying the finite element method whenever necessary. A simple example of the tiling procedure is shown in Fig. 7.3. During the tiling procedure three different types of tiles can be defined according to their shape: interior, edge, and corner tiles. The interior tile has a cross-shaped lumped equivalent circuit having four branches (or ports), and this tile can be also used to represent the edge and corner tiles by appropriately terminating one or two branches. On the other hand, tiles can be classified into three different primitive structures according to their functionality, as illustrated in Fig. 7.3: power/ground plane tiles, power/signal/ground tiles and power/ground pin tiles. The equivalent circuits for these primitives are derived by using appropriate electromagnetic simulation tools, as will be discussed below. In addition, the equivalent circuit representations for the input current and voltage sources as well as connectors are considered. Throughout

189 PWR / GND Tile D DC' nm Dx A n.m B B L=L C= C'm ~, = Rnm" A' C' I'-! 7 - - - - - - - - - - - - - - - - - I I ile i -; - I [^{ --- —----- Signal ne I I I Dn'm ' I - - r- - -t - - -- i~~i lii!iiiii;i!~ l i l iilii;!ii!=ili l~ii;*:ii =~i i;:~~~iii;~~i i~ii~iii i~iiiii; z { I D~' { I J ~D I I I PWR/GND Pin Tile II I I I I I I +1 - 7 c m B Q PWR n'm' L= L SGL N C= Cnm' R= Rn" m Q' GND PWR / SGL/GND Tile Figure 7.3: Equivalent circuits for three different tiles, such as power/ground tile, power/ground pin tile, and power/signal/ground tile.

190 the present study it has been assumed that the connectors are attached to one side of the PCB to maintain constant potentials on power and ground planes. For PCBs having multiple power/ground planes, as commonly encounted in highspeed digital circuits, the equivalent circuits for each tile can also be stacked vertically to model the three-dimensional nature of the geometry. As a result, the tiling approach has no inherent limitations in terms of number of layers, except an increased computational burden. The information from all of these equivalent circuits are recorded in an HSPICE input file, and then circuit simulation is performed in time or frequency domain. From the time domain circuit simulation, the potential distributions on the power and ground planes are obtained, and from these potential distributions the switching noises can be estimated. The input impedance at the driving port can be also easily computed from the voltage and current information at a given electrical node. 7.3.2 Equivalent Circuits As mentioned above, a variety of modeling approaches are followed in order to identify the equivalent circuits for the primitives defined during the tiling procedure. Descriptions of these approaches are given below. (a) Power/Ground Plane Tiles A power/ground plane primitive is considered as a section of a parallel plate waveguide. From microwave circuit theory, the equivalent circuit for such a structure can be represented as a 4-port network, as shown in Fig. 7.3. Under the assumption that the parallel plate waveguide supports only a TEM wave, while evanescent higher order modes are weakly excited and can be practically considered non-existent, the equivalent circuits can be easily computed using the ideal transmission line theory.

191 The above underlying assumptions are, in general, valid for clock frequencies up to the GHz range. Furthermore, to take into account the floating potentials on the power and ground planes, a special form of an equivalent circuit is devised, as shown in Fig. 7.3. Since an arbitrary surface current can be decomposed into two orthogonal components, the parallel plate waveguide is decomposed into two orthogonal 2-port transmission lines corresponding to each direction of the current, and as a result a 4-port cross shaped equivalent network is obtained. For the 4-port network, the appropriate equivalent inductances L m^ and L' m for two orthogonal directions and the capacitances Cnm are derived from the 2-port ideal transmission line as: 307rh wL = - sin(/3qDn'm) (7.2) WLm s in( nDym) (7.3) Cn,m [Dnm tan(t D nm) n,m tan( Dn) (74) 120rh Y 2 2 Y where DJ and D-,n 1 < i < N and 1 < j < M are the x and y dimensions of the (i j)th tile as shown in Fig. 7.3. N and M are the number of tiles in x and y direction, respectively. Furthermore, h and E. are the thickness and relative dielectric constant of the dielectric material separating the conducting layers of the printed circuit board and,3g is the propagating constant at the frequency of equivalent circuit extraction. In addition to the inductive and capacitive elements, resistive elements are also added to take into account the effect of the ohmic losses and to prevent pure inductive loops in the equivalent circuits. The value of the resistive element can be calculated as: R= W/0 (7.5) bo2a where a is the conductivity of the material.

192 Under the assumption that the electrical dimensions of the tile are much smaller than the wavelength at the operating frequency, that is, O/gD 'm and,SgD 'm < 1, the above equations can be reduced to a more compact form using the small argument approximation of the trigonometric functions: 307'h (On,'n L D, = (Dnm (7.6) c y 307rh Dn ( __ = 7 ~ - n,m (Dnm x Dn,) (7.8) 2407rhc where c is the speed of light in free space. As can be recognized from the above expressions, the inductance and the capacitance become frequency independent as both the size of each tile and the distance between the power and ground planes become relatively small compared to the smallest wavelength of interest in the propagation signals. As frequency increases, however, the inductance and capacitance values may vary accordingly, giving rise to high-frequency effects. For the examples presented in this study, we assume that the inductances and capacitances in all equivalent circuits are independent of frequency. (b) Power/Signal/Ground Tiles The signal line located between two conducting planes is modeled as a stripline (as illustrated in:Fig. 7.3). Similar to the power/ground plane tile, the equivalent circuit for a signal line is developed to take into account the effect of the signal line on the potential fluctuations at the power and ground planes by using the series inductances and shunt capacitances as shown in the figure. Derivation of the equivalent circuit is completed in two steps. First, the characteristic impedance of the stripline is evaluated from the empirical formula [229] which includes the effect of the fringing fields. Second, an equivalent T-network is constructed for a stripline

193 of length Dn'"r and width Dy', as shown in the bottom of the figure, and the values of the capacitance and inductance in the equivalent circuit are determined after taking a procedure similar to that in the previous subsection. The values of the inductances and capacitances for this primitive can be calculated as follows: wLn',m = znm' tan( 3g' ) (7.9) 2 WCn'm' = n3D'm') (7.10) o1,,m, where =T ' 3 _(7.11) Vcn (We/H) + 0.441 We D"'m' 0 if Dn"m' /H > 0.35 H H (0.35- )2 if Dy'm/IH < 0.35 After the equivalent circuit components are determined, those elements are augmented into the equivalent circuit network of the power/ground planes: that is, P - P' is connected to A - A' and Q - Q' to B - B'. By doing this, the interaction between the signal line and the power/ground plane can be fully accommodated. In general, the equivalent inductance L"n m' and capacitance Cn",m' are functions of frequency. However, if the length of the transmission line tile is smaller than that of the wavelength corresponding to the clock frequency, then /3gD"m' << 1, thus leading into frequency independent circuit components. As a result, the above equations are simplified to the following: Ln' - Zm-D'onm'-V (7.13) 2c C"m' O"x r (7.14) Z ',m c where Zo m' is the characteristic impedance of the signal line.

194 (c) Power/Ground Pin Tiles To model the power and ground pins connected to a current source or driver, as shown in Fig. 7.3, a simple equivalent circuit is employed. The power pin, directly connected to the power plane, is modeled as a vertical inductance element, while the ground pin, which is connected to the ground plane through a via hole on the power plane, is replaced by two inductances, one for the line above the power plane and one for the line below. Furthermore, to increase the accuracy of the model, one capacitor is added between the ground pin and the power plane, as shown in Fig. 7.3. This capacitor represents the parasitic capacitive effect between the pin and the opening in the plane. The appropriate inductance and capacitance values are computed for a given geometry from an electromagnetic simulator, the finite element method (FEM) [29, 30], developed in the previous chapters. The FEM is based on tetrahedral edge basis functions and provides a highly accurate model of the fields including their dynamic effects. The values of the inductance and capacitance in the equivalent circuit are extracted by matching the S-parameters derived from the FEM to the S-parameters of the equivalent circuit representing the power/ground pin tiles. A similar type of equivalent circuit could be derived from a simplified analytic expressions, as presented in [169]. It has been found in the previous chapter that the inductance of the pin is a function of the length and diameter of each pin, and in practice, the inductance can be considered frequency independent. The amount of the capacitive effect between the pin and the via hole in the power plane is mainly determined by the diameters of the pin and surrounding hole. Note that as the diameter of the vertical via or pin increases, the inductance of the structure decreases due to the increased width of the current path.

195 Top View Top View Side View P \ P G:r25pF |3.4 Via Hole 50 mm Figure 7.4: Example PCB used in the validation experiments. 7.4 Simulation Results The modeling and simulation methodology described above was used to analyze SSN for a number of configurations of an example PCB whose physical characteristics are shown in Fig. 7.4. These configurations differed in the locations of the power/ground connectors, the number and placement of the current sources and sinks, the amount of current injected by each driver, and the size and location of decoupling capacitors. The effect of a signal trace has been also modeled as a stripline in one of the experiments. Using the flow outlined in Fig. 7.2, a PCB configuration is tiled and simulated in the time domain using HSPICE. Transient simulation results, in the form of voltage waveforms at the power and ground grid points, have been subsequently converted to produce 3-D noise maps. In these maps, the maximum drop in voltage at each grid point on the power plane and the maximum increase in voltage at each grid point on the ground plane have been plotted, unless otherwise specified. In general, depending on physical importance, several different types of noise maps can be generated, such

196 as average, RMS, peak-to-peak, and maximum or minimum values on the power and ground planes. 7.4.1 Convergence Test: Tiling Optimization The first task is to determine how changing the number of tiles used to model the power and ground planes affects the noise maps, and the experimental setup is shown in Fig. 7.5(a). The connector is placed on the left hand side and a current source is connected between the power and the ground planes as shown. The point "p"' denotes the connection to the power plane and the point "g" denotes the connection to the ground plane through a via hole. The current waveform, shown in Fig.7.5(a), is a triangular pulse with rise (tr) and fall time (tf) of 500 ps and a peak value (Imax) of 10 mA. This wave shape is chosen to model the actual current flowing through the n and p-type transistors of a switching CMOS inverter. The noise maps for the various tiling cases are shown in Fig.7.5 (b) and (c) for the power and ground plane, respectively. These maps indicate a maximum voltage drop of 2.73 mV on the power plane (corresponding to point "p") and a maximum voltage increase of 3.69 mV on the ground plane (corresponding to "g") regardless of the number of tiles. As can be seen from the maps, the essential features of the noise maps remain the same, with maps corresponding to finer grids showing better resolution of the noise magnitude. Fig. 7.6 illustrates the CPU and memory usage as a function of the number of tiles for this configuration. The CPU and memory requirement for the circuit simulation increases rapidly as the number of tiles increases. As a result, 625 tiles for the partitioning of the PCB is chosen as a reasonable compromise between the accuracy and simulation time for all other experiments described in this section, except for the Intel P6 board.

197 ii I i(t) A [mA] Imax cG Pt t f t[nsec] (a) 225 tiles 625 tiles 1024 tiles 50 - " 50 - " — " ' I 0 - -'" 0 0 10 mm] 0 10 ] Y [mm]0o 110 2mm]3 Y[m ]100 10 2m3 40 50 (b) 5, 10.. — — '"'" 10 I 5, a, I 4 0 10 Qmm Ymm, ~g 030 4[m3 0m 1 0 210 30m30 40 50 (c) Figure 7.5: (a) Board layout and current waveform for validation experiments. Noise maps on the (b) power plane and (c) ground plane for three different number of tiles.

198 [(1 I --- 4000 (0 a) ' 3000 E D 2000 n1000 0 ~ ~ ~, 50 00 0 nl L,,r 0 500 1000 0 500 Number of Tiles Number of Tiles Figure 7.6: CPU and memory usage versus the number of tiles. 1000 0o c Co connector G P iO) G Pi) G connector Figure 7.7: Variation in the position of the connector.

199 Power Plane Max(Vp~p) = 5.516(mV] Ground Plane Max(Vp) =7.814[MV] 30 -~"20 ~ 10 Y [mm] 40 50 0 0 X [MM] MaxOV-pp) = 23.1 1 [mV.] Max(Vpp) = 21.22[mV] (a) (C) 50 (b) Max(VLpp) = 5.51 6[mV] 50 (d) Figure 7.8: Noise map (peak-to-peak voltage variations) on the power and ground planes with connector in various places: (a) left, (b) top, (c) right, (d) bottom.

200 The position of the connectors is then varied keeping the coordinates of the current injection points on the PCB identical, as shown in Fig. 7.5. The connectors are alternately placed on the left, top, right and bottom edges of the PCB (see Fig. 7.7). As can be observed in Fig. 7.8, the potential fluctuations on the power and ground planes have their maximum values when the connectors are placed on the top and right edges. This is readily explained by noting that the current injection points for these two configurations are farther away from the DC power supply than in the other two configurations, making them more susceptible to noise. 7.4.2 Seven Drivers with Decoupling Capacitors In order to obtain more realistic values of the noise voltages, the PCB with actual transistor drivers replacing the ideal piecewise linear current source has been simulated. Seven large CMOS inverters were connected to different points on the power and ground planes. Transistor sizes and capacitive loads have been chosen to represent realistic VLSI chip pad drivers: 1400 um/1.2 gm for p-channel MOSFETs, 800,um/1.2,m for n-channel MOSFETs, and 25 pF load capacitance. Each driver is excited by a 1.7 ns wide voltage pulse with 0.8 ns rise/fall times. The resulting noise maps are shown in Fig. 7.9 (a) and indicate that the maximum noise levels on the power and ground planes are -439.26 mV and 560.93 mV, respectively. In order to quantify the effect of the decoupling capacitors, 100 nF ideal capacitances (i.e. with no lead inductance) are connected across the power and ground leads of the seven drivers as illustrated in Fig. 7.10. The corresponding noise maps are shown in Fig. 7.9 (b) and indicate a reduction in the switching noise to -149.4 mV and 98.43 mV, respectively, on the power and ground planes. It is observed that the decoupling capacitance connected between the power and ground leads of the

201 Power Plane Min(Vmin) = -439.3[mV] Ground Plane Max(V_max) = 560.9m '4 50 40 (a) Min(V_min) = -149.4[mV] Max(V_max) = 98.43[mV] 0 50 40 X [mm] X [mm] (b) Figure 7.9: Noise maps (minimum on power and maximum on ground plane) (a) without and (b) with decoupling capacitors. Note that the scale of the vertical axis is not identical for both cases. driver allows the most effective noise reductions on the power and ground planes. Variation in noise levels as a function of the magnitude of decoupling capacitance is shown in Fig. 7.11. In general, larger values of decoupling capacitance contribute to reduce the noise level. However, switching noise cannot be reduced to zero regardless of the amount of decoupling capacitance used. Instead, the noise tends to saturate to a specific level beyond a certain critical value of the decoupling capacitance (the "knee" of the curves). While the critical value has been turned out roughly 1.0 nF with the given configuration, it is possible to have different values for different

202 i(t) Ground Plane Figure 7.10: Configuration of the decoupling capacitance. -100 '> E -200 I-I 'i-300 -400 10 I. I I...... I I I,,,,. I I,, rv.. I I.,,..., I. I...,.. I. I... I. I )-2... I.... I........ I........ I. I......... I I....... 10-1 100 101 C_decap [nF] 102 103 104 R o o. I I I I..... I I. I.. I.. 1.1.. I I 111. I UVj > E400 x co E I -2 9cc200 I I)I. I...-........ I...................................... 10-2 10-1 100 101 C_decap [nF] 102 103 104 Figure 7.11: Noise variation tances. on power and ground planes with various decoupling capaci

203 configurations. 7.4.3 Modeling of Signal Traces In this section, the effect of a signal trace on the switching noise is studied. For the experiment, one signal line having width 0.5 mm and length 42 mm is placed in the middle of the dielectric substrate ([r = 3.4) and a current source is connected at the center of the line. The current driver and its excitation pulse waveform are identical to those used in the previous experiment. In addition, one end of the signal line is terminated with 50 Q lumped resistor, while the other end of the line is connected to a trapezoidal voltage source having 0.8 nsec rise/fall time and 1.7 nsec pulse width, and its peak value is 5 V. The PCB configuration used in this experiment and the equivalent circuit of the stripline with current source are shown in Fig. 7.12. The resulting noise maps, minimum potential on the power and maximum on the ground planes, are shown in Fig. 7.13 for the both cases: with and without stripline. As can be observed in the figure, the presence of the stripline or signal trace causes a minor perturbations (i.e., reduced potential difference between the power and ground planes) to the potential distributions in the power and ground planes, especially in the ground plane. When the stripline is present, the maximum noise levels in the power and ground planes become -93.5 mV and 141.0 mV, respectively. However, when the stripline is removed, the noise levels in the power and ground planes are slightly changed to -87.0 mV and 168.6 mV, respectively, due to the reduced overall inductance and current injection into the planes. Note that those peaks on the power and ground planes correspond to the locations of the driver pins connected to each plane. When the width of the signal line becomes narrower, the interactions between the signal line and the power/ground planes are expected to be reduced due to the

204 Top View I.. - E 4mm 42 mm t v ipu 50 ~ to driver input A Side View 0 i u rw cd E E E E i G I, 50mm Source/Stripline Configuration Figure 7.12: Top and cross sectional views (not to scale), and source/stripline configuration. The parenthesized numbers are coordinates in [mm] relative to lower left corner of the PCB. V denotes connection to ideal ground, i.e., node 0 of HSPICE.

205 Power Plane 0. - Ground Plane.*- s c E.-1 > -200 50 4( Y[n 0 E -100 -.> -200 50 50 (a) 50 (b) Figure 7.13: Noise maps on the power and ground planes (a) with signal trace and (b) without signal trace. decreasing mutual capacitances. As a result, for very narrow signal traces, the effect of the traces can be considered as a secondary phenomenon. 7.4.4 Split PCBs As shown in Fig. 7.14, a given PCB layer is split into two or more planes to provide better noise immunization and resulting reduced SSN. This type of structure can be modeled by introducing gap capacitances between two conducting planes in a close proximity, and the appropriate capacitance values, Cg and Cp, are extracted from the equivalent circuit of the microstrip gap geometry [230, 231] as illustrated in Fig. 7.15:

206 Vss3 Figure 7.14: Split PCB and short circuit pins. Cg i c/ (a) (b) Figure 7.15: Gap in a microstrip line (a) and its lumped equivalent circuit (b).

207 for r = 9.8 and 0.5 < w/h < 2, Cp = Ceven (7.15) C9 = (Codd - Ceven)/2 (7.16) where Codd/w = (s ) e [pF/m] w Ceven / = ( S)Me [pF/m] w with m = -w[0.619 log(w/h)- 0.3853] for 0.1 < s/w < 1.0 ko = 4.26- 1.4531og(w/h) for 0.1 < s/w < 1.0 me = 0.8675, ke = 2.043( -)012 for 0.1 < s/w < 0.3 h 1.565 0.03 m = - 1, k-= 1.97- 0 for 0.3 < s/w< 1.0 (w/h)0.16 (w/h) For any other dielectric material in the range 2.5 < cr < 15, the the even and odd mode capacitances can be calculated by using the following factors: Codd(Cr) = Codd(9.6)[9 ]r.8 9.6 Ceven(Cr) = Codd(9.6)[ ] 0.9 9.6 After the equivalent capacitances are found, those elements are augmented into the original network, and the inductance elements in the gap region are then removed. The large void area in the central region is also modeled by removing the electrical nodes corresponding to that area. On the other hand, for the modeling of the short circuit pins as shown in Fig. 7.14, an appropriate inductor is allocated to each pin. To quantify the effect of the gaps and short circuit pins in the split PCB, the 64 mm x 64 mm PCB is split into two sections as illustrated in Fig. 7.16 and the gap

208 i C 11.,C. C f. iC " C e e e gg 9 9 gg9 gg9 gg9 Figure 7.16: Configuration of the split PCB having 25 short circuit pins ('') and 11 current; sources (Imax = 100 mA and 0.5 nsec rise/fall time). The dots ('.') represent electrical nodes on the planes. Note that the connector is placed at the left hand side of the board supplying 0.0 V only for the lower plane. 0.6. 0.3 15 10 ". '. 00 14 1515 0 20 0 20 (a) (b) Figure 7.17: Noise maps on the (a) upper (Vmax = 420 mV) and (b) lower planes (Vmax = 260 mV).

209 spacing s is chosen to be 2.54 mm. In addition, 25 short circuit pins are distributed evenly in the left and right hand sides, and 11 current sources are connected in the central region of the left hand side (denoted as '@' also). The magnitude of the current sources is 100 mA and the rise and fall times are fixed to 0.5 nsec. Fig. 7.17 depicts the potential fluctuations on the upper and lower planes, and reveal the effect of the short circuit pins. The short circuit pins have a tendency to elevate the potential on the lower plane and lower on the upper plane. The potential fluctuation on the left hand side is much smaller than that of the right hand side both on the upper and lower planes because of the connector providing constant potential. Furthermore, the potential difference between the two split planes on the upper plane suggests the effectiveness of the gaps for reducing the overall noise level. 7.4.5 Intel P6 Gadget For the validation of the proposed approach, the noise level is measured under the control environment for the Intel P6 gadget as illustrated in Fig. 7.18. To minimize the effect of the measurement errors, multiple drivers are provided to switch simultaneously, and then the voltage difference between the power and ground planes at various locations are measured and compared with the modeled result. For the modeling, the lower half of the P6 board is discretized into number of tiles (34 x 37 non-uniform tiles) and 64 triangular current sources are injected simultaneously between the power and ground pins. Those 64 sets of power and ground pins form two groups (32 in each group) and each aggregated source current group is located in the upper left and right Conner of the chip area. The input current waveform and its frequency component are shown in Fig. 7.19. Note that the Fourier components in the signal span up to near Gigahertz range.

-dam a m am onmumwmmm -mum mu mu mumV 00 CJo -ummum~mA somm -7C7IW I I I II111iF

211 Current Pulse in Time Domain 0 5 10 15 20 time [nsec] Frequency Component in the Pulse 25 30 2 2.5 3 Frequency [GHz] Figure 7.19: Input current waveform (0.8 nsec rise and fall time) and its frequency component. Fig. 7.20 shows the comparision between the measured and simulated voltage waveforms at 4 different locations denoted as AS41, AP28, AC35, and AL27. These measurement points are located in the upeer right corner of the chip area. As can be observed in the voltage waveforms, the simulated results follow the measurement very closely, even for the small details. With this comparison, the validity and effect of the proposed approach have been proved. Note that in the measurement, voltage sources are used to excite the board, while in the simulation, the equivalent current sources are utilized. 7.4.6 Summary From the experiments performed in this chapter, the following observations can be made. First, for a given amount of current injection, the SSN increases when the current drivers are placed farther away from the equipotential connectors that

212 2.9 - -2.8 < 2.7.... > 2.6...... Measurement 2.5 - Simulation 0 5 10 15 20 25 30 2.9 2.78................. 2.6. - Measurement 2.5...................... 2.5 - Simulation 2.4 -0 5 10 15 20 25 30 2.9 < 2.7 > 2.6 Measurement 2.5 - * * * * * *.... * *........................................... sim ulti o....... 2.5 Simulation........... 2.4 - I I I i 2.9 0 5 10 15 20 25 30 > 2.6 -..........*............................ - Measurement 2.5 - -- - 1 Simulation | 2.4 0 5 10 15 20 25 30 Figure 7.20: Output current waveforms at 4 different points with 64 input current sources switching simultaneously: magnitude lrna = 15 mA and rise/fall time tr = 2f = 0.8 nsec

213 supply DC power to the board. Second, maximum reduction in the SSN occurs when decoupling capacitors are placed in the vicinity of the current injection points (i.e. between the power and ground driver pins). Third, increasing the value of the decoupling capacitance tends to reduce the magnitude of the SSN. Beyond a critical value, however, further increases in capacitance value does not reduce the overall SSN. This critical value depends on the PCB configuration and the amount of current injected by the drivers. Fourth, the presence of the signal traces marginally contribute to noise on the power and ground planes. Limited evidence, however, suggests that this contribution is dwarfed by the current injection from the drivers. Additional experiments also suggested that significant lead inductances introduced by the decoupling capacitors may in some cases nullify the noise reduction benefits of these capacitors. Fifth, the effective of splitting the power or ground planes are examined giving the evidence of suppressed noise. Finally, the actual 8 layer Intel P6 board is modeled and the excellent comparison between the measurement and simulation proves the accuracy of the hybrid approach. 7.5 Conclusions This chapter has presented a systematic methodology, combining the finite element full-wave electromagnetic tool and circuit simulator HSPICE, for the modeling and analysis of SSN in digital ICs and packages. In this hybrid methodology, the seemingly contradictory goals of modeling accuracy and global analysis efficiency are reconciled through a divide-and-conquer process. The accuracy is insured by performing detailed EM field analysis on appropriately chosen small sections of the PCB, referred to as tiles, to create lumped electrical equivalent circuit models. The PCB tile models are subsequently combined with the models for chip current and/or

214 voltage drivers and package leads to produce an electrical simulation model for the PCB power distribution subsystem. A circuit simulator is then used to exercise this model under a variety of current or voltage excitation conditions to yield noise maps that indicate the variation in the power and ground potential as a function of location on the PCB. Furthermore, the effect of the decoupling capacitances and splitting the PCB are simulated, and discovered that the noise levels on the power and ground planes are saturated to a certain value as the decoupling capacitance increases. Finally, a realistic circuit, Intel P6 gadget, is modeled and found that the agreement between the measurement and simulations is very satisfactory, thus proving the effect and accuracy of the proposed hybrid approach.

CHAPTER VIII CONCLUSIONS 8.1 Summary of Achievements The main theme of this dissertation is the development of a rigorous but efficient and flexible three-dimensional full-wave electromagnetic modeling tool for the characterization of high-speed high-frequency interconnects ranging from a simple via hole to an entire millirmeter-wave package. For this purpose, the edge-based finite element method with tetrahedral sub-domain elements has been successfully developed and parallelized on a distributed memory parallel computer, the IBM SP2. It is worth mentioning that the use of the edge-based basis functions renders the FEM solutions free from spurious modes which until now have been a bottleneck of the node-based FEM. To overcome the limitation imposed on general partial differential equation techniques regarding the truncation of the computational domain at a finite size, an isotropic absorbing material is utilized for the simulation of open boundary problems, and its performance is proved to be satisfactory in MMIC applications, which have inhomogeneous material configurations. Furthermore, the parallelization of the developed FEM technique makes the time consuming full-wave analysis more affordable and attractive as a design and optimization tool for high-speed high-frequency interconnects. 215

216 The development of a stable and rigorous FEM technique raises the possibility of real time characterization of complicated microwave and millimeter-wave devices and circuits which was not feasible before. In this research effort, the parallelized FEM technique has been applied to the characterization of various MMIC components: single and double via holes, plated-through-hole, CPW-to-microstrip transitions, CPW-to-slot transitions, CPW-to-waveguide transition, and hermetic wall and bead transitions. Furthermore, a K/Ka-band MMIC package with dozens of vertical via holes, DC bias lines, bonding wires and gaps is thoroughly investigated, and the internal package resonances are identified along with a few possible treatments of the phenomenon. For the via hole geometries, appropriate equivalent circuits are derived to represent the structures as a lumped circuit. These equivalent circuits are then used in system level electromagnetic modeling of digital ICs and packages in Chapter 7 to provide equivalent circuits for tiles. One of the greatest advantage of using the finite element method is its capability of field reconstructions at any computational domain. Field visualization by means of these field reconstructions in a domain offers a thorough understanding of the electromagnetic phenomena occurring in reality and helps to diagnose any potential circuit malfunctions. Hybrid approach combining electromagnetic full-wave analysis and circuit simnulation has also been proposed for the system level modeling of digital ICs and packages. As presented in Chapter 7, remarkable agreement between the simulation results and measurement data has been achieved for an 8-layer INTEL P6 gadget. Considering the fact that there are hundreds of signal traces and vertical via holes, the correlation between the predicted and measured data proves the effectiveness

217 of the proposed technique. Furthermore, the simultaneous switching noise and its suppression by placing decoupling capacitors are also closely studied and visualized. Since the technique does not assume the actual ground planes to be perfect, the potential fluctuations on the power as well as ground planes are easily simulated. 8.2 Suggested Future Work Several possible directions could be taken as future research: * refinement and improvement of the finite element method itself * application of the FEM for the characterization of structures which are not feasible to characterize using other methods * efficient parallelization on shared as well as distributed memory machines * further improvement of the hybrid method presented in Chapter 7 including automatic interfacing between circuit layout files and tiling processor. First of all, the refinement and improvement of the finite element method includes a few interesting topics, such as the development of a more sophisticated and adaptive three-dimensional mesh generator, improvement of absorbing boundary conditions, application of a more efficient linear equation solver, more effective and interactive electromagnetic field visualization, etc. Furthermore, passive and active lumped electrical elements can be also incorporated into the FEM framework. Second, the finite element method can be applied to the characterization of other interesting high-speed high-frequency circuits and interconnects as well as packages. The recent trend of using the electromagnetic simulator for system level applications suggests a possible direction. Since most circuit and system designers heavily

218 depend on their experience and intuitions which might be misleading, rigorous theoretical data can contribute to a more thorough understanding of the electromagnetic phenomena occurring in and around circuits. Third, effective parallelization of the finite element method on different types of parallel computers based on shared or distributed memory can speed up the time consuming full-wave analysis and could lead to the full-wave analysis tool becoming a design and optimization vehicle. Another important task is the effective parallelization of conjugate gradient type linear equation solvers. Since these types of solvers utilize repeated matrix-vector multiplications, efficient parallel matrix-vector multiplication could enhance the overall performance of the FEM technique. Finally, the hybrid approach combining the FEM and circuit simulator for system level modeling of digital ICs and packages can be evolved to its next generation by developing automatic interfacing between circuit layout and tiling processors, constructing a more accurate and extensive equivalent circuit library, and automatic placement of decoupling capacitors to minimize the simultaneous switching noise. Furthermore, for the accurate simulation of the non-linear characteristics of the drivers, more sophisticated driver modeling will be necessary.

APPENDICES 219

220 APPENDIX A Calculation of the Element Matrix Entries The final matrix equation of the FEM could be read as follows: [A]T[x] = [b] where M [A] = Ui[Ue]-k2[Ve]} e=l M [x] = Uxe e=l M [b] = U[b] e=l and M denotes total number of tetrahedrons. The number of linear equations, which is equivalent to the dimensions of matrix [A], is in general different from the number of tetrahedrons in the domain. The entries of the elementary FEM matrices for each tetrahedron could be evaluated by analytical and numerical integrations of the inner product with the aid of basic integration formulas for shape functions'. The final expressions for each matrix entry of the equation are derived as follows: ie e2V/e Aij = Uij - k Vij Jf f fSV( )a(A2)b(A3)'(A4)ddzdydZ = (a+b+c+d+3)! 6Ve

221 where L[j = (/u71VxWVxW3) l^ = (-WVx,WV)) t. ( Tl Wi 7 Wj) and 1 < ij < 6. Considering the fact that the gi and gj are constant vectors, the evaluation of the Ui can be performed directly by using one of the properties of the edge basis function, equation (2.17), and yields: U.j = 4Ve7-l(gi- gj). The V/- may be broken into three terms as ^ =,r (Pi + P2 + P3) where Pi = (f, fj) = Ve(fi fj) P2 = (fi,gj x r) + (g x r, fj) P3 = (gi X r,gj x r). Note that the computation of P1 can be perfomed easily since the fi and fj are constant vectors by definition. Evaluation of the P2 and P3 can also be done analytically: V= P2 - (TXSX + TySy + TzSz) 4 and ve P3= 2 20 + + (giygjy + gizgjz)[Sxx + (Sx)2] (gixgjx + gizgjz)[Syy + (Sy)2] (gixgjZ + gigz)[Syy + (Sy)2]

222 - (gtjxgjy + 9IYjX[X + SrSy] - (g, ygj.z + 9,gyly + Sys'] -(gixgjz + gzg.r[X + SXSZ]} where = 9i 4 4 4 4 4 4 4 4 4 =x ZXiyi SYZ=Zyi, ISx=Z x'zz i=1 i=1 i=1 In the above, the constant vectors T and gj are computed numerically for a given tetrahedron. The -volume of a tetrahedron with four nodal coordinates of the vetexes, (xiyz) 1,...,4, is given by determinant of 4 x 4 matrix as follows: 1X1 Yi Z1 ye = -det X 2Z 6 1X X3 Y/3 Z3 1X4 Y4Z

223 APPENDIX B Pre-conditioned Conjugate Orthogonal Conjugate Gradient Method For the solution of linear equation Ax = b with sparse complex symmetrix matrix A, the conjugate orthogonal conjugate gradient (COCG) method [58] is implemented in this study. The COCG method is different from other projection type methods, such as the congugate gradient (CG) method, conjugate gradient squared method and GMRES, since these methods are based on forming an orthogonal basis for Krylov subspace1, but not in the COCG case. For the non-Hermitian symmetric matrix, the orthogonality cannot be achived in a simple way. The key step in the COCG method is the new definition of the orthogonality2 for the residual vectors: (VmnVn) = 6mn where (x,y)= = xjyj. The subscripts m and n in the vectors, xm and xn, denote the vectors in the mth and nth iterations. Now, with the new definition of the conjugate orthogonality in mind, lThe j-dimensional Krylov subspace KI(A: ro,) for a linear system Ax = b with non-singular A is defined as the span of the vectors {vo, Avo,...,AJ-1Vo} for a given initial residual vo = ro = b - Axo. 2The definition of the Hermitian orthogonality in the inner product space is given as: (Xm, Xn) = 6mn

224 the COCG algorithm can be derived from the standard conjugate gradient procedure [48]-[57] by replacing the complex Hermitian inner product with the bilinear form [x, y] = j xjy as summarized in the following: Initialization Xo: given vO = b-Axo; p-i = 0; /-1 = 0; Wo = K- vo; po= [Vo,,w]; For j = 0, 1,2,... iterations PI uj a3 = Wj- +j-lPj-1 = Apj = [uj,pj]; if pj = 0, then quit. = Failure! - pj/j. Xj+i = xj -+ ojpj; It xj+i Is accurate enough, then quit. =L Convergence Vj+1 = - wj+i = K-1Wj+i; Pj+1 = [vj+i,w j+l]; if pj+l is small, then quit. =i Failure pj = Pj+l /pj Next j Note that K is used as a pre-conditioning matrix and in this study the diagonal entries of the matrix A are used, i.e. [IKij] = [Aj Sij]. As a result, the inverse of K can be computed without additional computations: [Iij]-~ = [A-1jlS]. Considering the fact that the most computional efforts in the CG type algorithm is devoted to the matrix-vector multiplication, the COCG method provides clear advantage over the

225 CG methods, since COCG method requires only one matrix-vector multiplication in each iteration instead two in the CG method. Furthermore, the matrix-vector multiplication of a sparse symmetric matrix can be highly optimized with minimum storage requirements. Parallelization of the matrix-vector multiplication routine on shared and distributed memory machines would be promising.

226 APPENDIX C Network Parameters To characterize a given circuit, several different network parameters are needed to be extracted. For example, for a vertical via hole, connecting two different signal lines on different layers or providing a ground path to a transmission line, the inductance of the via hole is a vital information for its characterization, and this inductance can be computed by two different methods. In the followings, two different methods of extracting equivalent via hole inductance are sketched, and the mathematical background of the extraction of the scattering parameters for one and two-port networks are briefly discussed. In general, the scattering parameters are used to characterize any given VLSI and MMIC interconnects. C.1 Inductance of Vertical Via Hole The inductance of a vertical via hole can be computed from several different methods, however, in this study, two different methods are utilized: current/flux method and scattering parameter method. In the flux/current method, the inductance of a via hole is calculated by taking the ratio of the current flowing on the via hole to the flux surrounding the current path as illustrated in Fig. C.1. The current flowing through the vertical via hole is evaluated from the closed path line integral as I = j H dC [Ampere]

227 I -- I 4 C _.~~' -C Hb. orthogonal plane to the input signal line Figure C.1: Vertical via hole as a current path. and the flux wounding the vertical current path is computed as = Jlj Lr H. dS [Wb/m2]. From the value of the current and flux computed above, the inductance of the vertical current path is computed as: L= [Henry]. This approach has a certain advantage to understand and characterize the effect of the vertical via hole, since the effect of the vertical current can be isolated from the other current components due to the orthogonality property of the electromagnetic fields. Thus inductance computed from the flux/current definition contains the effect of the via only. In the second approach, the scattering parameters and lumped equivalent circuit network are used to extract the equivalent via hole inductance. The scattering parameters are computed from the voltage or current standing wave patterns on the input transmission line as will be discussed in the following sections. After the

228 + V 1-port Network F(X), Zin(x) + I — Vl 0 - 2-port Network.I2 + ---.47y Jl(x), Zin,l (X) Figure C.2: Equivalent circuits for one- and two-port network. scattering parameters are obtained, those values are matched with the frequency response of the chosen equivalent circuit. Note that the obtained inductance and capacitance values comprise the overall effect beyond the reference plane and also depend on the position of the reference plane. C.2 One-Port Network To characterize any one-port network (refer to Fig. C.2), the input reflection coefficient is determined from the voltage standing wave pattern, and to obtain unique values of the phase and magnitude of the reflection coefficient, the position of the reference plane has to be predefined and fixed. In this syudy, all the scattering pa

229 rameters are computed with the reference plane chosen at the position where the discontiniuty starts. First, from the voltage maximum and minimum values the voltage standing wave ratio (VSWR) is determined: W lmaxI I VminI Note that for a lossless circuit, the VSWR has infinite value. After the VSWR is computed, the reflection coefficient looking towards the discontinuity is determined anvwhere on the transmission line as: VSWR-1 -I2 _g(X-ax) (X) - VSWR + 1e where xmax is the position of a voltage maximum. If a circuit is lossless, then the magnitude of r(x) becomes unity, and as a result, the phase of the reflection coefficient r(x) becomes an important factor. Moreover, the normalized input impedance z(x) at any point on the transmission line can be determined from the reflection coefficient F(x) as: (X) =1 + F(x) z - (x) C.3 Symmetric Two-Port Network The network parameters for symmetric two-port networks can be determined by analyzing the voltage or current distributions on the input and output ports with the even and odd mode excitations. In general, any symmetric two-port network can be represented as shown in Fig. C.2 and characterized by the following relationship between the voltage and current at the ports ([232], pp.216-220) V1 Zll z12 11 V2 Z21 Z22 I2

230 PEC or PMC wall gl1 or 4 Ig9 1or Ig2 12 e,o e,o Zin(X) Zn(X) Figure C.3: Symmetric two port network with the even and odd mode excitations corresponding to placing magnetic and electric walls at the center of the circuit. First of all, to extract the scattering parameters the even and odd mode input impedances (Z^, z9) are computed from the two different excitations, even and odd mode. From these values, the impedance parameters are determined according to ze + z~ Ze - Z in n in -in Zll Z 12= 2 2 2 Note that for the passive symmetric two-port reciprical networks, the z11 is equal to Z22, and z21 to z12. When the excitation sources are current as shown in Figure C.3, the even mode excitation (Ig = Ig2 = Io) corresponds to placing a perfect magnetic conducting (PMC) wall at the center of the circuit, since the voltage standing wave pattern has maximum at the center and the voltage maximum is equivalent to open circuit, thus leading to PMC wall. Similarly, the odd mode excitation (Ig, = -Ig2 Io) corresponds to placing the PEC wall at the center of the circuit. The scattering parameters are then obtained from the normalized impedance parameters by using the following equations: 121 - Z12 22 - Z121 + 2Z11 - z22 + 1

231 312 - S21.2 - _ —22 + 1 - 11 — 12 1 where S11 is equal to S22. and S12 to S21, because the 2-port network is reciprocal and symmetric. C.4 Asymmetric Two-Port Network In order to characterize any asymmetric two-port network, three different excitation mechanisms are needed to obtain the impedance parameters z11, z22 and z-1 with the method employed in this study. As three different excitation pairs, the even (Igi = Ig2 = Io), odd (Ig = -Ig2 = Io) and sum (Ig- = 1, Ig2 = 2) excitation mechanisms are utilized in this study. From these excitations, the normalized impedance parameters are evaluated by using the following relations [233] 11 = ls + (Zls -lo)(Zls - Zle) (2 - Z2o) A 22 = Z2s - (Z2s - Z2o)(Z2s - Z2e)(- 1) / Z2e - 22 Z12 = ~(Z )le -Zll) VZle- Z11 where A = (Zls - Zlo)(Z2s - Z2e) - (Zls - Zle)(Z2s - Z2o). In the above expressions, zle, Z2e, zlo1, 20, Zls and Z2s are the normalized input impedances at the port 1 and 2 for even, odd and sum excitation pairs, respectively. The sign of the z12 should be carefully decided as discussed in [233]. The scattering parameters for asymmetric 2-port network can be determined with the following expressions: 2 -1- 22 C + 1 ) + 12 11 1)(22 1) - 2 (Z11 + 1)(Z22 + 1)- 2

232 -2 _1 _ 2 S22 2 (11 2 + 1)(222 +1)- - 2Z12 S12 - 2 (ZI + 1)(22 + 1)Note that for reciprocal 2-port network the S12 is equal to the S21, and for symmetric network, S11 equals S22. For a lossless 2-port network, IS112 + S1212 = 1.

BIBLIOGRAPHY 233

234 BIBLIOGRAPHY [1] R. F. Harrington, Field computation by moment methods. New York: Macmillan, 1968. [2] P. P. Silvester and R. L. Ferrari, Finite elements for electrical engineers. Cambridge: Cambridge university press, 1991. [3] J. Jin, The finite element method in electromagnetics. New York: John Wiley & Sons, inc, 2 ed., 1993. [4] K. S. Yee, "Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media," IEEE Ant. Propagat., vol. 14, pp. 302-307, 1966. [5] M. Krumpholtz and P. Russer, "Two-dimensional FDTD and TLM," Int. J. of Num. Model.: Electron. Net. Dev. Fields, vol. 7, pp. 141-153, 1994. [6] M. Krumpholtz, C. Huber, and P. Russer, "A field theoretical comparison of FDTD and TLM," IEEE Trans. Microwave Theory Tech., vol. MTT-43, pp. 1935-1950, 1995. [7] 0. C. Zienkiewicz, The finite element method: vol.1 Basic formulation and linear problems. London: McGraw-Hill Book Company, 4 ed., 1994. [8] F. Kikuchi, "Mixed and penalty formulation for finite element analysis of an eigenvalue problem in electromagnetism," Computer Methods in Applied Mechanics and Engineering, vol. 64, pp. 509-521, Oct. 1987. [9] B. M. A. Rahman and J. B. Davies, "Penalty function improvement of waveguide solution by finite elements," IEEE Trans. Microwave Theory Tech., vol. MTT-32, pp. 922-928, Aug. 1984. [10] J. P. Webb, "The finite element method for finding modes of dielectric-loaded cavities," IEEE Trans. Microwave Theory Tech., vol. 33, pp. 635-639, July 1985. [11] J. P. Webb, "Efficient generation of divergence-free fields for the finite element analysis of 3D cavity resonances," IEEE Trans. Magnetics, vol. 24, pp. 162-165, Jan. 1988. [12] B. B. Davies, F. A. Fernandez, and G. Y. Philippou, "Finite element analysis of all modes in cavities with circular symmetry," IEEE Trans. Microwave Theory Tech., vol. MTT-30, pp. 1975-1980, Nov. 1982.

235 [13] A. Bossavit. '"A rationale for edge-elements in 3-D fields computations," IEEE Trans. Magnetics. vol. 24, pp. 74-79, Jan. 1988. [14] A. Bossavit, "Whitney forms: A class of finite elements for three-dimensional computations in electromagnetism." Proc. IEE, vol. 135, Pt. A, pp. 493-500. Nov. 1988. [15] A. Bossavit and I. Mayergoyz, "Edge-elements for scattering problems," IEEE Trans. Magnetics, vol. 25, pp. 2816-2821, July 1989. [16] A. Bossavit, "Solving Maxwell equations in a closed cavity, and the question of spurious modes," IEEE Trans. Magnetics, vol. 26, pp. 702-705, Mar. 1990. [17] A. Bossavit, "Simplicial finite elements for scattering problems in electromagnetism," Computer Methods in Applied Mechanics and Engineering, vol. 76, pp. 299-316, Nov. 1989. [18] J. P. Webb, "Edge elements and what they can do for you," IEEE Trans. Magnetics, vol. 29, pp. 1460-1465, Mar. 1993. [19] M. L. Barton and Z. J. Cendes, "New vector finite elements for three-dimensional magnetic field computation," J. Appl. Phys., vol. 61, pp. 3919-3921, April 1987. [20] J.-G. Yook and L. katehi, "Characterization of MIMICs packages using a parallelized 3D FEM code," in Proc. of Advance in Computational Electromagnetics Symposium (ACES), pp. 954-961, 1996. [21] J.-G. Yook, M. Kurk, N. Dib, L. Katehi, and T. Arabi, "Evaluation of ground inductance in printed circuit boards," in Proc. of Topical Meeting on Electrical Performance and Electronic Packaging (EPEP), pp. 119-121, Oct., 1993. [22] J.-G. Yook, V. Chandramouli, L. Katehi, and K. A. Sakallah, "Computation of switching noise in printed circuit boards," in Proc. of Progress in Electromagnetic Research Symposium (PIERS), p. 864, July 1995. [23] J.-G. Yook, V. Chandramouli, L. P. B. Katehi, and K. Sakallah, "Simultaneous switching noise in printed circuit boards: Electromagnetic modeling/simulation," Radiation Lab. Report RL-912, The University of Michigan, Ann Arbor, MI, Feb. 1995. [24] J.-G. Yook, V. Chandramouli, L. Katehi, and K. A. Sakallah, "Computation of switching noise in PCBs for digital packages," in Proc. of Topical Meeting on Electrical Performance and Electronic Packaging (EPEP), Oct.2-4, 1995. [25] J.-G. Yook, V. Chandramouli, L. Katehi, and K. A. Sakallah, "Computation of switching noise in printed circuit boards," IEEE Trans. Components, Packaging, and Manufacturing Technology: Part A (accepted), 1997. [26] J.-G. Yook and L. Katehi, "Characterization of CPW to microstrip transitions using the finite element method," in IEEE Ant. Propagat. Soc. Int. Symp. Dig., vol. 17, pp. 1206-1209, May, 1994.

236 [27] J.-G. Yook, N. Dib. L. Katehi, R. Simons. and S. Taub, "Theoretical and experimental study of microstrip-to-slot line uniplanar transition," in IEEE Ant. Propagat. Soc. Int. Syrnp. Dig., pp. 1206-1209, June, 1994. [28] E. Tentzeris, M. Krumpholz, N. Dib, J.-G. Yook, and L. Katehi, "FDTD characterization of waveguide probe structure," IEEE Trans. Microwave Theory Tech. (submitted), 1997. [29] J.-G. Yook, N. Dib, and L. Katehi, "Characterization of high frequency interconnects using finite difference time domain and finite element methods," IEEE Trans. Microwave Theory Tech., vol. MTT-42, pp. 1727-1736, Sept. 1994. [30] J.-G. Yook and L. Katehi, "Capacitance and inductance evaluation of vertical interconnects using the finite element method," in Int. URSI Symp. Dig., p. 33, June, 1993. [31] J.-G. Yook, N. Dib, E. Yasan, and L. Katehi, "A study of hermetic transitions for microwave packages," in IEEE MTT-S Int. Microwave Symp. Dig., pp. 1579-1582, May, 1995. [32] J.-G. Yook, N. Dib, E. Yasan, and L. Katehi, "Modeling of hermetic transitions for microwave packages," Int. J. of Microwave and Computer-Aided Engineering, vol. 6, pp. 351-368, 1996. [33] J.-G. Yook, L. Katehi, R. Simons, and K. Shalkhauser, "Experimental and theoretical study of parasitic leakage/resonance in a K/Ka-band MMIC package," IEEE Trans. Microwave Theory Tech., vol. 44, Dec. 1996. [34] J.-G. Yook, L. Katehi, R. Simons, and K. Shalkhauser, "Experimental and theoretical study of parasitic leakage/resonance in a K/Ka-Band MMIC package," in IEEE MTT-S Int. Microwave Symp. Dig., pp. 223-226, 1996. [35] X. Yuan, "Three-dimensional electromagnetic scattering from inhomogeneous objects by the hybrid moment and finite element method," IEEE Trans. Microwave Theory Tech., vol. MTT-38, pp. 1053-1058, Aug. 1990. [36] J. Jin, J. L. Volakis, and J. D. Collins, "A finite element-boundary integral method for scattering and radiation by two- and three-dimensional structures," IEEE Ant. and Propagat. Magazine, vol. 33, pp. 22-32, June 1991. [37] A. Chatterjee, J. Jin, and J. L. Volakis, "Computation of cavity resonances using edge-based finite elements," IEEE Trans. Microwave Theory Tech., vol. MTT-40, pp. 2106-2108, Nov. 1992. [38] A. Chatterjee, J. Jin, and J. L. Volakis, "Edge-based finite elements and vector ABC's applied to 3-D scattering," IEEE Ant. Propagat., vol. AP-41, pp. 221-226, Feb. 1993. [39] J. S. Wang and R. Mittra, "Finite element analysis of MMIC structures using absorbing boundary conditions," in IEEE MTT-S Int. Microwave Symp. Dig., pp. 745-748, 1993.

237 [40] K. Hayata. Mi. Koshiba, MI. Eguchi, and MI. Suzuki,! Vectorial finite-element method without any spurious solutions for dielectric waveguiding problems using transverse magnetic-field component," IEEE Trans. Microuwave Theory Tech.. vol. MTT-34. pp. 1120-1124, Nov. 1986. [41] E. Yasan, J.-G. Yook, and L. katehi, "Modeling the interface between fields and devices using the finite element method," in Proc. of Progress in Electromagnctic Research Symposium (PIERS), 1996. [42] J.-G. Yook, N. Dib, and L. Katehi, "Modeling of microwave and millimeter-wave packages," in Proc. of Progress in Electromagnetic Research Symposium (PIERS). p. 460, July, 1995. [43] W. Schroeder and I. Wolff, "The origin of spurious modes in numerical solutions of electromagnetic field eigenvalue problems," IEEE Trans. Microwave Theory Tech.. vol. 42, pp. 644-653, April 1994. [44] A. F. Peterson, D. R. Wilton, and R. E. Jorgenson, "Variational nature of Galerkin and non-Galerkin moment method solutions," IEEE Ant. Propagat., vol. 44, pp. 500 — 503, April 1996. [45] Z. Ren, F. Bouillault, A. Razek, A. Bossavit, and J. C. Verite, "A new hybrid model using electric field formulation for 3-D eddy current problems," IEEE Trans. Magnetics, vol. 26, pp. 470-473, Mar. 1990. [46] J.-F. Lee, "Analysis of passive microwave devices by using three-dimensional tangential vector finite elements," Int. J. of Num. Model.: Electron. Net. Dev. Fields, vol. 3, pp. 235-246, 1990. [47] J.-F. Lee, G. M. Wilkins, and R. Mittra, "Finite-element analysis of axisymmetric cavity resonator using a hybrid edge element technique," IEEE Trans. Microwave Theory Tech., vol. MTT-41, pp. 1981-1987, Nov. 1993. [48] M. R. Hestenes and E. Stiefel, "Methods of conjugate gradients for solving linear systems," J. of Research of the National Bureau of Standards, vol. 49, pp. 409-436, Dec. 1952. [49] C. Lanczos, "Solution of systems of linear equations by minimized iterations," J. of Research of the National Bureau of Standards, vol. 49, pp. 33-53, July 1952. [50] G. H. Golub and C. F. V. Loan, Matrix Computations. New York: The Johns Hopkins University Press, 1985. [51] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in FORTRAN. Cambridge: Cambridge University Press, 1992. [52] R. W. Freund, "Conjugate gradient-type methods for linear systems with complex symmetric coefficient matrices," SIAM J. on Scientific and Statistical Computing, vol. 13, pp. 425-448, Jan. 1992. [53] H. A. van der Vorst, "BI-CGSTAB: a fast and smoothly converging variant of BICG for the solution of nonsymmetric linear systems," SIAM J. on Scientific and Statistical Computing, vol. 13, pp. 631-644, March 1992.

238 [54] T. K. Sakar and E. Arvas, "On a class of finite step iterative methods (conjugate directions) for the solution of an operator equation arising in electromagnetics," IEEE Ant. Propagat., vol. AP-33, pp. 1058-1066, Oct. 1985. [55] R. E. Kleinman and P. M. van den Berg, "Iterative methods for solving integral equation," Radio Science, vol. 26, pp. 175-181, 1991. [56] I. S. Duff, "A survey of sparse matrix research," Proc. IEEE, vol. 65, pp. 500-535. 1977. [57] C. F. Smith, A. F. Peterson, and R. Mittra, "The bi-conjugate gradient method for electromagnetic scattering," IEEE Ant. Propagat., vol. AP-38, pp. 938-940, June 1990. [58] H. A. van der Vorst and J. B. M. Melissen, "A Petrov-Galerkin type method for solving Ax = b, where A is symmetric complex," IEEE Trans. Magnetics, vol. 26, pp. 706-708, March 1990. [59] W. Schroeder and I. Wolff, "A hybrid-mode boundary integral equation method for normal- and superconducting transmission lines of arbitrary cross-section," Int. J. of Microwave and Computer-Aided Engineering, vol. 2, pp. 314-330, 1992. [60] J. A. M. Svedin, "A numerically efficient finite-element formulation for the general waveguide problem without spurious modes," IEEE Trans. Microwave Theory Tech., vol. MTT-37, pp. 1708-1715, Nov. 1989. [61] D. R. Lynch and K. D. Paulsen, "Origin of vector parasites in numerical Maxwell solutions," IEEE Trans. Microwave Theory Tech., vol. MTT-39, pp. 383-394, Mar. 1991. [62] K. D. Paulsen and D. R. Lynch, "Elimination of vector parasites in finite element maxwell solutions," IEEE Trans. Microwave Theory Tech., vol. MTT-39, pp. 395 -404, Mar. 1991. [63] D. R. Lynch, K. D. Paulsen, and W. E. Boyse, "Synthesis of vector parasites in finite element Maxwell solutions," IEEE Trans. Microwave Theory Tech., vol. MTT-41, pp. 1439-1448, Aug. 1993. [64] D. R. Lynch. K. D. Paulsen, and J. W. Strohbehn, "Finite element solution of Maxwell's equations hyperthermia treatment planning," J. Computational Phys., vol. 58, p. 246269, 1985. [65] X. Yuan, D. R. Lynch, and K. D. Paulsen, "Importance of normal field continuity in inhomogeneous scattering calculations," IEEE Trans. Microwave Theory Tech., vol. MTT-39, pp. 638-642, April 1991. [66] I. Bardi and 0. Biro, "An efficient finite-element formulation without spurious modes for anisotropic waveguide," IEEE Trans. Microwave Theory Tech., vol. MTT-39, pp. 1133-1139, July 1991. [67] M. Wong, O. Picon, and F. Hanna, "Spurious-free finite element analysis of multiaxial coplanar waveguide discontinuities," in Proc. of European Microwave Conference, pp. 334-337, 1993.

239 [68] MI. Wong. 0. Picon, and V. F. Hanna, "Three dimensional finite element analysis of N-port waveguide junctions using edge-elements," in IEEE MTT-S Int. Microuave Symp. Dig., pp. 417-420, 1992. [69] M. Hano, "Finite-element analysis of dielectric-loaded waveguides," IEEE Trans. Microwave Theory Tech., vol. MTT-32, pp. 1275-1279, Oct. 1984. [70] A. R. Pinchuk, C. W. Crowley, and P. P. Silvester, "Spectrally correct finite element operators for electromagnetic vector fields," J. Appl. Phys., vol. 63, no. 8, pp. 3025 — 3027, April '!88. [71] T. Angkaew, M. Matsuhara, and N. Kumagai, "Finite-element analysis of waveguide modes: A novel approach that eliminates spurious modes," IEEE Trans. Microwave Theory Tech., vol. MTT-35, pp. 117-123, Feb. 1987. [72] A. Konrad, "On the reduction of the number of spurious modes in the vectorial finite element solution of three-dimensional cavities and waveguides," IEEE Trans. Microwave Theory Tech., vol. MTT-34, pp. 224-227, Feb. 1986. [73] S. H. Wong and Z. J. Cendes, "Combined finite element-modal solution of threedimensional eddy current problems," IEEE Trans. Magnetics, vol. 24, pp. 2685-2687, Nov. 1988. [74] H. Whitney, Geometric integration theory. Princeton: Princeton University Press, 1957. [75] J. P. J. R. Fraser, "Unified approach to problems in electromagnetism," Proc. IEE, Pt. A, vol. 131, pp. 55-61, Jan. 1984. [76] J. D. Baldomir, "Differential forms and electromagnetism in 3-dimensional Euclidean space r3," Proc. IEE, Pt. A, vol. 133, pp. 139-143, May. 1986. [77] G. A. Deschamps, "Electromagnetics and differential forms," Proc. IEEE, vol. 69, pp. 139-143, June 1981. [78] Z. L. Zhu, J. 13. Davies, and F. A. Fernandez, "3-D edge element analysis of dielectric load resonant cavities," Int. J. of Num. Model.: Electron. Net. Dev. Fields, vol. 7, pp. 35-41, 1994. [79] R. S. H. Hoole, S. Jayakumaran, and S. Yoganathan, "Tetrahedrons, edges and nodes in 3-D finite-element meshes," Electron. Lett., vol. 22, pp. 735-737, July 1986. [80] B. H. McDonald and A. Wexler, "Finite-element solution of unbounded field problems," IEEE Trans. Microwave Theory Tech., vol. MTT-20, pp. 841-847, Dec. 1972. [81] X. Yuan, D. R. Lynch, and J. W. Strohbehn, "Coupling of finite element and moment methods for electromagnetic scattering from inhomogeneous objects," IEEE Ant. Propagat., vol. AP-38, pp. 386-393, March 1990. [82] S. D. Gedney and R. Mittra, "Analysis of the electromagnetic scattering by thick gratings using a combined FEM/MoM solution," IEEE Trans. Microwave Theory Tech., vol. MTT-39, pp. 1605-1614, Nov. 1991.

240 [83] S. D. Gedney, J. F. Lee. and R. Mittra. "A combined FEM/MoMNI approach to analyze the plane wave diffraction by arbitrary gratings," IEEE Trans. Alicrouave Theory Tech., vol. MTT-40, pp. 363-370, Feb. 1992. [84] T. Nakata, N. Takahashi, K. Fujiwara, and M. Sakaguchi, "3-D open boundary magnetic field analysis using infinite element based on hybrid finite element method," IEEE Trans. Magnetics, vol. MAG-26, pp. 368-370, Mar. 1990. [85] S. J. Salon and J. D. Angelo, "Applications of the hybrid finite element-boundary element method in electromagnetics," IEEE Trans. Magnetics, vol. MAG-24, pp. 80 -85, Mar. 1988. [86] Z. Ren, F. Bouillault, A. Bossavit, and J. C. Verite, "A new hybrid model using electric field formulation for 3-D eddy current problems," IEEE Trans. Magnetics. vol. MAG-26, pp. 470-473, Mar. 1990. [87] J.-F. Lee and Z. j. Cendes, "Transfinite elements: A highly efficient procedure for modeling open field problems," J. Appl. Phys., vol. 61, pp. 3913-3915,15 April 1987. [88] B. Engquist and A. Majda, "Absorbing boundary conditions for the numerical simulation of waves," Mathematics of computation, vol. 31, pp. 629-651, July 1977. [89] Y.-C. Ma, "A note on the radiation boundary conditions for the Helmholtz equation," IEEE Ant. Propagat., vol. AP-39, pp. 1526-1530, Oct. 1991. [90] A. Bayliss, M. Gunzburger, and E. Turkel, "Boundary conditions for the numerical solution of elliptic equations in exterior regions," SIAM J. on Applied Mathematics, vol. 42, pp. 430-450, April 1982. [91] C. H. Wilcox, "An expansion theorem for electromagnetic fields," Communications on Pure and Applied Mathematics, vol. IX, pp. 115-134, 1956. [92] J. D'Angelo and I. D. Mayergoyz, "Finite element methods for the solution of RF radiation and scattering problems," Electromagnetics, vol. 10, pp. 177-199, 1990. [93] J. D'Angelo and I. D. Mayergoyz, "On the use of local absorbing boundary conditions for RF scattering problems," IEEE Trans. Magnetics, vol. 25, pp. 3040-3042, IJuly 989. [94] A. F. Peterson, "Absorbing boundary conditions for the vector wave equation," Microwave and Opt. Tech. Lett., vol. 1, pp. 62-64, 1988. [95] A. F. Peterson, "Accuracy of 3-D radiation conditions for use with the vector Helmholtz equation," IEEE Ant. Propagat., vol. 40, pp. 351-355, 1992. [96] J. P. Webb and V. N. Kanellopoulos, "Absorbing boundary conditions for the finite element solution of the vector wave equation," Microwave and Opt. Tech. Lett., vol. 2, pp. 370-372, 1989. [97] J. P. Webb, "Absorbing boundary conditions for the finite-element analysis of planar devices," IEEE Trans. Microwave Theory Tech., vol. 38, pp. 1328-1332, 1990.

241 [98] R. L. Higdon, "Numerical absorbing boundary conditions for the wave equation." Mathematics of computation, vol. 49, pp. 65-90. July 1987. [99] R. L. Higdon, "Absorbing boundary conditions for acoustic and elastic waves in stratified media," J. Computational Phys., vol. 101. pp. 386-418, 1992. [100] R. Janaswamy, "2-D radiation boundary conditions on an arbitrary outer boundary," Microwave and Opt. Tech. Lett., vol. 5, pp. 393-395, 1992. [101] G. Y. Kriegsmann, A. Taflove, and K. R. Umashankar, "A new formulation of electromagnetic wave scattering using an on-surface radiation boundary condition approach," IEEE Ant. Propagat., vol. 35, pp. 153-161, 1987. [102] T. G. Moore, J. G. Blaschak, A. Taflove, and G. Y. Kriegsmann, "Theory and application of radiation boundary operators," IEEE Ant. Propagat., vol. 36, pp. 1797-1812, 1988. [103] J. G. Blaschak and G. Y. Kriegsmann, "A comparative study of absorbing boundary conditions," J. Computational Phys., vol. 77, pp. 109-139, 1988. [104] J. B. Keller and D. Givoli, "Exact non-reflecting boundary conditions," J. Computational Phys., vol. 82, pp. 172-192, 1989. [105] C. K. W. Tam and J. C. Webb, "Radiation boundary condition and anisotropy correction for finite difference solutions of the Helmholtz equation," J. Computational Phys., vol. 113, pp. 122-133, 1994. [106] 0. M. Ramahi and R. Mittra, "Finite-element analysis of dielectric scatterers using the absorbing boundary condition," IEEE Trans. Magnetics, vol. 25, pp. 3043-3045, July 1989. [107] L. W. Pearson, R. A. Whiaker, and L. J. Bahrmasel, "An exact radiation boundary condition for the finite-element solution of electromagnetic scattering on an open domain," IEEE Trans. Magnetics, vol. 25, pp. 3046-3048, July 1989. [108] J. F. Imhoff, G. Meunier, and J. C. Sabonnadiere, "Finite element modeling of open boundary problems," IEEE Trans. Magnetics, vol. MAG-26, pp. 1659-1661, Sept. 1990. [109] J. F. Imhoff, G. Meunier, X. Brunotte, and J. C. Sabonnadiere, "An original solution for unbounded electromagnetic 2D- and 3D-problems throughout the finite element method," IEEE Trans. Magnetics, vol. MAG-26, pp. 1659-1661, Sept. 1990. [110] Y. R. Brauer, S. M. Schaefer, J.-F. Lee, and R. Mittra, "Asymptotic boundary condition for three dimensional magneto static finite elements," IEEE Trans. Magnetics, vol. MAG-27, pp. 5013-5015, Nov. 1991. [111] B. Stupfel, "Absorbing boundary conditions on arbitrary boundaries for the vector wave equations," IEEE Ant. Propagat., vol. AP-42, pp. 773-780, 1994. [112] B. Stupfel and R. Mittra, "A theoretical study of numerical absorbing boundary conditions," IEEE Ant. Propagat., vol. AP-43, pp. 478-487, May 1995.

242 [113] A. Khebir, A. B. Kouki. and R. Mittra. "Asymptotic boundary condition for finite element analvsis of three-dimensional transmission line discontinuities," IEEE Trans. Microwave Theory Tech., vol. MTT-38, pp. 1127-1432, Oct. 1990. [114] A. Khebir, A. B. Kouki, and R. Mittra, "Higher order asymptotic boundary condition for the finite element modeling of two-dimensional transmission line structures." IEEE Trans. Microwave Theory Tech., vol. MTT-38, pp. 1433-1438, Oct. 1990. [115] A. Khebir, R. Mittra, and 0. Ramahi, "Numerically derived absorbing boundary condition for the solution of open region scattering problems," IEEE Ant. Propagat.. vol. AP-39, pp. 350-353, Mar. 1991. [116] T. B. A. Senior, J. L. Volakis, and S. R. Legault, "Higher order impedance and absorbing boundary conditions," IEEE Ant. Propagat., to be appear. [117] R. Mittra, O. Ramahi, A. Khebir,, and R. Gordon, "A review of absorbing boundary conditions for two and three-dimmensional electromagnetic scattering problems,' IEEE Trans. Magnetics, vol. 25, pp. 3034-3039, July 1989. [118] J.-P. Berenger, "A perfectly matched layer for the absorption of electromagnetic waves," J. Computational Phys., vol. 114, pp. 185-200, 1994. [119] K. Li, J. T. Johnson, J. J. Akerson, R. T. Shin, and J. A. Kong, "Theoretical and numerical analyses of Berenger's PML," in Proc. of Progress in Electromagnetic Research Symposium (PIERS), July 1995. [120] C. M. Rappaport, "Perfectly matched absorbing boundary conditions based on anisotropic lossy mapping of space," IEEE Microwave and guided wave letters, vol. 5, pp. 90-92, Mar. 1995. [121] J.-M. Jin, J. L. Volakis, and V. V. Liepa, "Fictitious absorber for truncating finite element meshes in scattering," Proc. IEE, Pt. H, vol. 139, pp. 472-476, Oct. 1992. [122] T. Ozdemir and J. L. Volakis, "A comparative study of an absorbing boundary condition and an artificial absorber for truncating finite element meshes," Radio Science, vol. 29, pp. 1255-1263, 1994. [123] S. R. Legault., T. B. A. Senior, and J. L. Volakis, "Design of planar absorbing layers for domain truncation in fer applications," Electromagnetics, vol. 16, pp. 451-464, 1996. [124] R. Mittra and U. Pekel, "A new look at the matched layer (PML) concept for the reflectionless absorption of electromagnetic waves," IEEE Microwave and guided wave letters, vol. 5, pp. 84-86, Mar. 1995. [125] Z. Wu and J. Fang, "Numerical implementation and performance of perfectly matched layer boundary condition for waveguide structures," IEEE Trans. Microwave Theory Tech., vol. 43, pp. 2676-2683, Dec. 1995. [126] Z. S. Sacks, D. M. Kingsland, R. Lee, and J.-F. Lee, "A perfectly matched anisotropic absorber for use as an absorbing boundary condition," IEEE Ant. Propagat., vol. 43, pp. 1460-1463, Dec. 1995.

243 [127] W. C. Chew and XW. H. Weedon, "'A 3D perfectly matched medium from modified Maxwell's equations with stretched coordinates," MAicrowave and Opt. Tech. Lctt.. vol. 7, pp. 599-604, Sept. 1994. [128] W. V. Andrew, C. A. Balanis. and P. A. Tirkas, "A comparison of the Berenger perfectly matched layer and the Linrnan higher-order ABC's for the FDTD method," IEEE Microwave and guided wave letters, vol. 5, pp. 192-194, June 1995. [129] J. R. Wait, "On the PML concept: a view from the outside," IEEE Ant. and Propagat. Magazine, vol. 38, pp. 48-51, April 1996. [130] J. D. Moerloose and M. A. Stuchly, "Behavior of Berenger's ABC for evanescent waves," IEEE Microwave and guided wave letters, vol. 5, pp. 344-346, Oct. 1995. [131] U. Pekel and R. Mittra, "An application of the perfectly matched layer (PML) concept to the finite element method frequency domain analysis of scattering problems," IEEE Microwave and guided wave letters, vol. 5, pp. 258-260, Aug. 1995. [132] R. H. Voelker and R. J. Lomax, "A finite difference transmission matrix method incorporating a nonlinear device model," IEEE Trans. Microwave Theory Tech., vol. 38, pp. 302-312, Mar. 1990. [133] I. Wolff, "Finite difference time-domain simulation of electromagnetic fields and microwave circuits," Int. J. of Num. Model.: Electron. Net. Dev. Fields, vol. 5, pp. 163 — 182, 1992. [134] T. Shibata and H. Kimura, "Computer-aided engineering for microwave and millimeter-wave circuits using the FD-TD technique of field simulations (invited article)," Int. J. of Microwave and Computer-Aided Engineering, vol. 3, pp. 238-250, 1993. [135] W. Sui, D. A. Christensen, and C. H. Durney, "Extending the two-dimensional FDTD method to hybrid electromagnetic systems with active and passive elements," IEEE Trans. Microwave Theory Tech., vol. 40, pp. 724-730, April 1992. [136] R. Gillard, J.-H. Corre, M. Driss, and J. Citerne, "A general treatment of matched terminations using integral equations - modeling and applications," IEEE Trans. Microwave Theory Tech., vol. 42, pp. 2545-2553, Dec. 1994. [137] J. S. Nielsen and W. J. R. Hoefer, "Modeling of nonlinear elements in a threedimensional condensed node TLM mesh," Int. J. of Microwave and Computer-Aided Engineering, vol. 3, pp. 61-66, 1993. [138] V. A. Thomas, M. E. Jones, and R. J. Mason, "Coupling of the PISCES device modeler to a 3-D Maxwell FDTD solver," IEEE Trans. Microwave Theory Tech., vol. 43, pp. 2170-2172, Sept. 1995. [139] B. Toland, B. Houshmand, and T. Itoh, "Modeling of nonlinear active regions with the FDTD method," IEEE Microwave and guided wave letters, vol. 3, pp. 33-35, Sept. 1993. [140] B. Tola a nd. in, B. Houshmand, and T. Itoh, FDTD analysis of an active antenna," IEEE Microwave and guided wave letters, vol. 3, pp. 423-425, Nov. 1993.

244 [141] V. A. Thomas, K.-NM. Ling, MI. E. Jones. B. Toland, J. lin. and T. Itoh, '"FDTD analysis of an active antenna." IEEE Microwave and guided wave letters, vol. 4. pp. 296-298, Sept. 1994. [142] V. A. Thomas, M. E. Jones, MI. Piket-MayV A. Taflove, and E. Harrigan. "The use of SPICE lumped circuits as sub-grid models for FDTD analysis," IEEE Microwave and guided wave letters, vol. 4, pp. 141-143, May 1994. [143] R. H. Voelker and C. G. Sentelle, "FD-TLM modeling of Josephson junction circuits," IEEE Microwave and guided wave letters, vol. 5, pp. 137-138, May 1995. [144] T. B. A. Senoir and J. L. Volakis, Approximate boundary conditions in electromagnetics. London: Institution of Electrical Engineers, 1995. [145] R. F. Harrington, Time-harmonic electromagnetic fields. New York: McGraw-Hill Book co., 1961. [146] R. F. Harrington and J. R. Mautz, "An impedance sheet approximation for thin dielectric shells," IEEE Ant. Propagat., vol. AP-23, pp. 531-534, July 1975. [147] N. Dib, J. W. Harokopus, G. Ponchak, and L. Katehi, "A comparative study between shielded and open coplanar waveguide discontinuities," Int. J. of Microwave and Computer-Aided Engineering, vol. 2, pp. 331-341, 1992. [148] J. C. Oliver and D. A. McNaramara, "Analysis of multiport discontinuities in waveguide using a pulsed FDTD approach," IEEE Trans. Microwave Theory Tech., vol. 42, pp. 2229-2238, 1994. [149] J. D. Montgomery, "Hybrid MIC north american markets," Microwave J., vol. 32, pp. 32-39, April 1989. [150] in Proc. of Topical Meeting on Electrical Performance and Electronic Packaging (EPEP), 1993. [151] N. VandenBerg and L. Katehi, "Broadband vertical interconnects using slot-coupled shielded microstrip lines," IEEE Trans. Microwave Theory Tech., vol. MTT-40, pp. 81-88, Jan. 1992. [152] T. Tanaka, K. Tsunoda, and M. Aikawa, "Slot-coupled directional couplers between double-sided substrate microstrip lines and their applications," IEEE Trans. Microwave Theory Tech., vol. MTT-36, pp. 1752 —1757, Dec. 1988. [153] A. Tran, B. Houshmand, and T. Itoh, "Analysis of microstrip lines coupled through an arbitrary circular aperture in a thick common ground plane," in Proc. of Topical Meeting on Electrical Performance and Electronic Packaging (EPEP), pp. 19-21, 1993. [154] A. Ittipiboon, D. Roscoe, and M. Cuhaci, "Analysis of slot-coupled double layer microstrip lines," in IEEE 1990 Conference Proceedings: Symposium on Antenna Technology and Applied Electromagnetics, pp. 454-459, Aug. 1990. [155] N. Herscovici and D. Pozar, "Full-wave analysis of aperture-coupled microstrip lines," IEEE Trans. Microwave Theory Tech., vol. MTT-39, pp. 1108-1114, July 1991.

245 [156] G. Luong and K. Chang. '"Analysis and development of a broadband interconnect for microstrip lines," Microwave and Opt. Tech. Lett.. vol. 5. pp. 388-393. July 1992. [157] R. Jackson, "Mode conversion at discontinuities in finite-width conductor backed coplanar waveguide," IEEE Trans. Microwave Theory Tech., vol. MITT-37. pp. 1582 -1589. Oct. 1989. [158] T. Wang, R. Harrington. and J. Mautz, "The equivalent circuit of a via" Tranls. Soc. Computer Simulations, vol. 4, pp. 97-123, 1987. [159] A. Djordjevic and T. Sakar, "Computation of inductance of simple vias between two striplines above a ground plane," IEEE Trans. Microwave Theory Tech., vol. MTT33, pp. 265-269, 1985. [160] T. Wang, R. Harrington, and J. Mautz, "Quasi-static analysis of a microstrip via through a hole in a ground plane," IEEE Trans. Microwave Theory Tech., vol. MTT36, pp. 1008 —1013, July 1988. [161] T. Wang, J. Mautz, and R. Harrington, "The excess capacitance of a microstrip via in a dielectric substrate," IEEE Trans. Computer-Aided Design, vol. 9, pp. 48-56, Jan. 1990. [162] M. E. Goldfarb and R. A. Pucel, "Modeling via hole grounds in microstrip," IEEE Microwave and guided wave letters, vol. 1, pp. 135-137, June 1991. [163] P. Kok and D. D. Zutter, "Capacitance of a circular symmetric model of a via. hole including finite ground plane thickness," IEEE Trans. Microwave Theory Tech., vol. MTT-39, pp. 1229-1234, July 1992. [164] J. Quine, H. Webster, and H. Glascock, "Characterization of via connections in silicon circuit boards," IEEE Trans. Microwave Theory Tech., vol. 36, pp. 21-27, Jan 1988. [165] T. Sakar, Z. Maricevic, and J. Zhang, "Evaluation of excess inductance and capacitance of microstrip discontinuities," in Proc. of Topical Meeting on Electrical Performance and Electronic Packaging (EPEP), pp. 104-106, 1992. [166] J. D. G. Swanson, "Grounding microstrip lines with via holes," IEEE Trans. Microwave Theory Tech., vol. MTT-40, pp. 1719-1721, Aug. 1992. [167] R. Sorrentino, F. Alessandri, M. Mongiardo. G. Avitabile, and L. Roselli, "Fullwave modeling of via hole grounds in microstrip by three-dimensional mode matching technique," IEEE Trans. Microwave Theory Tech., vol. 40, pp. 2228-2234, Dec. 1992. [168] P. Harms, J. Lee, and R. Mittra, "Characterizing the cylindrical via discontinuity," IEEE Trans. Microwave Theory Tech., vol. 41, pp. 153-156, Jan. 1993. [169] Q. Gu, Y. Yang, and M. A. Tassoudji, "Modeling and analysis of vias in multilayered integrated circuits," IEEE Trans. Microwave Theory Tech., vol. MTT-41, pp. 206 -214, Feb. 1993. [170] W. D. Becker, P. H. Harms, and R. Mittra, "Time-domain electromagnetic analysis of interconnects in a computer chip package," IEEE Trans. Microwave Theory Tech., vol. MTT-40, pp. 2155-2163, Dec. 1992.

246 [171] R. Simons. G. Ponchak, K. Martzaklis, and R. Romanofskv. -Channelized cpw: D)iscontinuities, junctions, and propagation characteristics," in IEEE MITT-S Int. 11icrowave Symp. Dig., pp. 915-918, 1989. [172] N. Dib, Characterization of Coplanar W1laveguide Transmission Lines and Discontinuities. PhD thesis, University of Michigan, 1992. [173] K. S. Yngvesson, T. L. Korzeniowski, Y.-S. Kimi, E. L. Kollberg, and J. F. Johansson. "The tapered slot antenna - a new integrated element for millimeter-wave applications," IEEE Trans. Microwave Theory Tech., vol. 37, pp. 365-374, Feb. 1989. [174] R. N. Simons, T. D. Perl, and R. O. Lee, "New coplanar waveguide feed network for 2 x 2 linearly tapered slot antenna subarray," Microwave and Opt. Tech. Lett., vol. 5, pp. 420-423, Aug. 1992. [175] R. N. Simons, R. Q. Lee, and T. D. Perl, "New techniques for exciting linearly tapered slot antennas with coplanar waveguide," Electron. Lett., vol. 28, pp. 620-621, Mar. 1992. [176] R. N. Simons, S. R. Taub, and P. G. Young, "Novel coplanar waveguide to slotline transition on high resistive silicon," Electron. Lett., vol. 28, pp. 2209-2210, Nov. 1992. [177] R. N. Simons, R. Q. Lee, and T. D. Perl, "Effect of a dielectric overlay on a linearly tapered slot antenna excited by a coplanar waveguide," Microwave and Opt. Tech. Lett., vol. 6, pp. 225-228, Mar. 1993. [178] R. N. Simons, S. R. Taub, R. Q. Lee, and P. G. Young, "Microwave characterization of slot line and coplanar strip line on high-resistivity silicon for a slot antenna feed network," Microwave and Opt. Tech. Lett., vol. 7, pp. 489-494, Aug. 1994. [179] R. N. Simons, N. I. Dib, and L. P. B. Katehi, "Coplanar stripline to microstrip transition," Electron. Lett., vol. 31, pp. 1725-1726, Sep. 1995. [180] S. C. Shi, J. Inatani, T. Noguchi, and K. Sunada, "Experimental investigation of the 'post mount' structure in waveguide-type SIS mixers," Proc. IEE, Pt. H, vol. 142, pp. 339-344, Aug. 1995. [181] R. B. Keam and A. G. Williamson, "Broadband design of coaxial line/rectangular waveguide probe transition," Proc. IEE, Pt. H, vol. 141, pp. 53-58, Feb. 1994. [182] T. F. McMaster, M. V. Schneider, and J. W. W. Snell, "Millimeter-wave receivers with subharmonic pump," IEEE Trans. Microwave Theory Tech., vol. 24, pp. 948 -952, Dec. 1976. [183] R. L. Eisenhart, "Discussion of a 2-gap waveguide mount," IEEE Trans. Microwave Theory Tech., vol. 24, pp. 987-990, Dec. 1976. [184] R. L. Eisenhart and P. J. Khan, "Theoretical and experimental analysis of a waveguide mounting structure," IEEE Trans. Microwave Theory Tech., vol. 19, pp. 706 -719, Aug. 1971.

247 [185] J. E. Oswald. Analysis of a millimeter waveguide mounting structure. Cambridge. MA: M.S. Thesis, Dept. of Electrical Engineering and Computer Science. Massachusetts Institute of Technology, 1992. [186] H. Bierman, "Designers strive for low cost MMIC packages," Microuave J., vol. 35. pp. 100-106, Sept. 1992. [187] S. Joshi and B. Becker, "Metal hermetic surface mount packaging," `Microwave J.. vol. 35, pp. 144-145, July 1992. [188] B. Ziegner, "High performance MMIC hermetic packaging," MIicrowave J., vol. 29. pp. 133-139, Nov. 1986. [189] D. Hughes and D. Jackson, "An equivalent circuit for a microwave surface mount package," Microwave J., vol. 34, pp. 134-144, Oct. 1991. [190] R. DeBoo, S. Hinderks, and N. Osbrink, "New surface-mounted package breaks from traditional MIC packaging," Microwave J., vol. 27, pp. 126-134, Mar. 1984. [191] D. Hughes, H. Harris,, and C. Rucker, "An equivalent circuit for a 70 mil microstrip package," Microwave J., vol. 29, pp. 97-102, Aug. 1986. [192] H. R. Malone, "Antenna/MMIC packaging techniques for commercial applications," in IEEE Ant. Propagat. Soc. Int. Symp. Dig., vol. 3, p. 1264, 1992. [193] L. Katehi, "The role of EM modeling in integrated packaging," in IEEE Ant. Propagat. Soc. Int. Symp. Dig., vol. 2, pp. 982-985, 1993. [194] D. S. Wein, "Advanced ceramic packaging for microwave and millimeter wave applications," in IEEE Ant. Propagat. Soc. Int. Symp. Dig., vol. 2, pp. 993-996, 1993. [195] H. Kuno, T. Milford, and J. Wooldridge, "The evolution of MMIC packaging," in IEEE Ant. Propagat. Soc. Int. Symp. Dig., vol. 2, pp. 1013-1016, 1993. [196] M. Ida and T. Nishikawa, "A multilayered package technology for MMICs," in IEEE Ant. Propagat. Soc. Int. Symp. Dig., vol. 2, pp. 1005-1007, 1993. [197] D. Griffin, "A new instrument and technique for diagnosing electromagnetic design problems in microwave module housings and component packages," IEEE Trans. Instrumentat. and Meas., vol. 41, pp. 27-31, Feb. 1992. [198] L. Lu, C. Wu, and J. Litva, "A new 3D-FDTD simulator for modeling electronic interconnections and packaging," in Proc. of Topical Meeting on Electrical Performance and Electronic Packaging (EPEP), pp. 92-95, 1993. [199] L. Dunleavy and P. Katehi, "Shielding effects in microstrip discontinuities," IEEE Trans. Microwave Theory Tech., vol. MTT-36, pp. 1767-1774, Dec. 1988. [200] L. Wu and Y. Chang, "Characterization of the shielding effects on the frequencydependent effective dielectric constant of a waveguide-shielded microstrip using the finite-difference time-domain method," IEEE Trans. Microwave Theory Tech., vol. MTT-39, pp. 1688-1693, Oct. 1991.

248 [201] C. Amrani. NI. Drissi. V. Fouad. and J. Citerne, "Packaging and interconnection mutual coupling effects in planar structures and discontinuities," in IEEE MITT-S Int. Microwave Symp. Dig., pp. 843-846, 1994. [202] P. Mezzanotte, M. Mongiardo. L. Roselli, and R. Sorrentino, "FDTD analysis of high performance MMIC package," in IEEE MTT-S Int. Alicrowave Syrmp. Dig.. pp. 337-340, 1994. [203] G. Strauss and W. Menzel, "A novel concept for mm-wave MMIC interconnects and packaging," in IEEE MTT-S Int. Microwave Symp. Dig., pp. 1141-1144, 1994. [204] R. Faraji-Dana and Y. Chow, "A CAD tool for interference rejection of a microwave circuit making use of response of the circuit packaging," in IEEE MTT-S Int. Microwave Symp. Dig., pp. 1523-1526, 1994. [205] M. Tsuji and H. Shigesawa, "Packaging of printed-circuit lines: A dangerous cause for narrow pulse distortion," IEEE Trans. Microwave Theory Tech., vol. MTT-42, pp. 1784-1790, Sept. 1994. [206] P. Mezzanotte, M. Mongiardo, L. Roselli, R. Sorrentino, and W. Heinrich, "Analysis of packaged microwave integrated circuits by FDTD," IEEE Trans. Microwave Theory Tech., vol. MTT-42, pp. 1796-1801, Sept. 1994. [207] R. Goyal, "Designing a 100 MHz SPARC dual processor using MCM-L packaging technology and microwave design techniques," in IEEE MTT-S Int. Microwave Syrnp. Dig., pp. 1715-1718, 1994. [208] M. Rittweger, M. Werthen, and I. Wolff, "3D FDTD analysis applied to the investigation of the resonant behavior of ceramic feedthrus," in IEEE MTT-S Int. Microwave Symp. Dig., pp. 1719-1722, 1994. [209] T. Shibata, S. Kimura, H. Kimura, Y. Imai, Y. Umeda, and Y. Akazawa, "Design technique for a 60 GHz-bandwidth distributed baseband amplifier IC module," IEEE J. Solid-State Circ., vol. 29, pp. 1537-1544, Dec. 1994. [210] K. Wu, R. Vahldieck, J. L. Fikart, and H. Minkus, "The influence of finite conductor thickness and conductivity on fundamental and higher-order modes in miniature hybrid MIC's (MHMIC's) and MMIC's," IEEE Trans. Microwave Theory Tech., vol. 41, pp. 421-430, March 1993. [211] M. Yu, R. Vahldieck, and J. Huang, "Comparing coax launcher and wafer probe excitation for 10 mil conductor backed CPW with via holes and air-bridges," in IEEE MTT-S Int. Microwave Symp. Dig., pp. 705-708, 1993. [212] D. W. Griffin and A. J. Parfitt, "Electromagnetic design aspects of packages for phased array modules that may incorporate monolithic antenna elements," in IEEE Ant. Propagat. Soc. Int. Symp. Dig., pp. 986-989, 1993. [213] J. Gipprich, L. Dickens, B. Hayes, and F. Sacks, "A compact 8-14 ghz LTCC stripline coupler network for high efficiency power combining with better than 82% combining efficiency," in IEEE MTT-S Int. Microwave Symp. Dig., pp. 1583-1586, 1995.

249 [214] R. P. Colwell, "Latent design faults in the development of multiflow's trace/200." in Proc. 22nd Fault Tolerant Computing Symposium, pp. 468-474, July 1992. [215] H. B. Bakoglu, Circuits, Interconnections and Packaging for IVLSI. VLSI Systems Series. New York: Addison-Wesley, 1990. [216] E. E. Davidson, "Electrical design of a high-speed computer packaging system," IEEE Trans. Components, Packaging, and Alanufacturing Technology: Part A, vol. 6. pp. 272-282, 1983. [217] G. A. Katopis, "Delta-I noise specification for a high-performance computing machine," in Proc. IEEE, pp. 1405-1415, 1985. [218] A. J. Rainal, "Computing inductive noise of chip packages," AT&T Bell Lab. Technical Journal, pp. 177-195, Jan. 1984. [219] R. Senthinathan and J. L. Prince, Simultaneous Switching Noise of CMOS Devices and Systems. Boston: Kluwer Academic Press, 1994. [220] L. W. Nagel, "SPICE2: A computer program to simulate semiconductor circuits," ERL-M520, University of California, Berkeley, Berkeley, CA, May 1975. [221] C. Dike, "Equivalent circuits model subtle traits of advanced CMOS ICs," Electronic Design News, pp. 189-198, April 1988. [222] W. S. Song and L. A. Glasser, "Power distribution techniques for VLSI circuits," IEEE J. Solid-State Circ., pp. 150-156, 1986. [223] W. T. Weeks, "Algorithms for ASTAP-A network analysis program," IEEE Trans. on Circuit Theory, vol. CT-20, pp. 628-634, 1973. [224] D. Shear, "EDN's advanced CMOS logic ground-bounce tests," Electronic Design News, vol. CT-20, pp. 88-114, March 1989. Correction: vol. 34, pp. 35, Sep. 1989.. [225] 0. Pedersen, B.-E. Flaten, and T. Gleditsch, "Modeling of power and ground planes," Circuit World, vol. 13, pp. 30-33, March 1992. [226] S. Maeda, T. Kashiwa, and I. Fukai, "Full wave analysis of propagation characteristics of a through hole using the finite-difference time-domain method," IEEE Trans. Microwave Theory Tech., vol. MTT-39, pp. 2154-2159, Dec 1991. [227] Meta-Software, HSPICE User's Manual, 1992. [228] L. T. Pillage, Asymptotic Waveform Evaluation for Timing Analysis. PhD thesis, Carnegie Mellon University, 1989. [229] H. H. Howe, Stripline circuit design. Dedham, MA: Artech House, 1974. [230] R. Garg and I. J. Bahl, "Microstrip discontinuities," Int. J. Electronics, vol. 45, pp. 81-87, 1978. [231] R. K. Hoffmann, Handbook of microwave integrated circuits. Norwood, MA: Artech House, 1987.

250 [232] D. Nl. Pozar. Alicrowave Engineering. New York: Addison-WXeslev Publislhing Company. 1990. [233] L. P. B. Katehi, "A generalized method for the evaluation of mutual coupling in microstrip arrays," IEEE Ant. Propagat., vol., 35, pp. 125-133. Feb. 1987.