ABSTRACT Title of Thesis: ELETROMAGNETIC MODELING WITH A NEW 3D ALTERNATING-DIRECTION-IMPLICIT (ADI) MAXWELL EQUATION SOLVER Degree Candidate: Xi Shao Degree and year: Master of Science, 2004 Thesis directed by: Professor Neil Goldsman Department of Electrical and Computer Engineering We introduce a time-domain method to simulate the digital signal propagation along on-chip interconnects, aperture radiation, and indoor-communication by solving the Maxwell equation with the Alternating-Direction-Implicit (ADI) method. With this method, we are able to resolve the large scale (i.e. electromagnetic wave propagation) and fine scale (i.e. metal skin depth, substrate current, coating material) structure in the same simulation, and the simulation time step is not limited by the Courant condition. The simulations allow us to calculate in detail parasitic current flow inside the substrate; propagation losses, skin-depth and dispersion of digital signals on non- ideal interconnects; detailed surface current and standing wave pattern in aperture radiation problem; signal power map and propagation delay in complicated in-door communication scenarios. ELECTROMAGNETIC MODELING WITH A NEW 3D ALTERNATING-DIRECTION-IMPLICIT (ADI) MAXWELL EQUATION SOLVER by Xi Shao Thesis submitted to the Faculty of the Graduate School of the University of Maryland in partial fulfillment of the requirements for the degree of Master of Science 2004 Advisory Committee: Professor Neil Goldsman, Chair Dr. Parvez N. Guzdar Professor Omar Ramahi ? Copyright by Xi Shao 2004 ii Dedication To My Parents and My Wife Who Made It All Possible iii Acknowledgements I greatly appreciate the encouragement, guidance and friendship from my advisors, N. Goldsman, P. Guzdar, and O. Ramahi. To Prof. Goldsman, I always remember his help during my most uncertain time and his restless curiosity. To Dr. Guzdar, I thank him for his readiness to answer or find answers for my questions. To Prof. Ramahi, without his advice and pushing, this project would not have gone this far. I also thank Prof. Papadopoulos, Dr. Sharma, Dr. Milikh, and Dr. Sitnov for their encouragements. I want to thank my wife for her patience and letting me do what I want. Thank my parents for asking me to finish whatever I have started. It is always a pleasure to discuss with Prof Goldsman?s students, Baiyun, Bo Yang, Akin, Zeynep, and Gary. iv Table of Contents Chapter 1 ......................................................................................................................1 Introduction..................................................................................................................1 1.1 Challenges for Modern Electromagnetic Modeling.............................................1 1.2 Thesis Structure ...................................................................................................4 Chapter 2 ......................................................................................................................5 Numerical Schemes of ADI Method and Code Validation.......................................5 2.1 Numerical Schemes of ADI Method....................................................................5 2.2 Model Validation ...............................................................................................12 2.2.1 2D Guided Electromagnetic Wave ......................................................12 2.2.2 2D Electromagnetic Wave Scattering..................................................13 2.3 Summary............................................................................................................16 Chapter 3 ....................................................................................................................17 Simulating on-chip Interconnects with ADI-FDTD Method .................................17 3.1 Introduction........................................................................................................17 3.2 Three Fundamental Propagation Modes for MISS Structure ......................18 3.3 Single Metal-strip-Insulator-Semiconductor-Substrate (MISS) Structure...27 3.4 Double Metal-strip-Insulator-Semiconductor-Substrate (MISS) Structure...................................................................................................................34 3.5 Summary............................................................................................................37 Chapter 4 ....................................................................................................................38 Simulating Electromagnetic Field Radiation from Aperture with ADI-FDTD Method ........................................................................................................................38 4.1 Introduction........................................................................................................38 4.2 Aperture Model and Simulation Configuration .................................................39 4.3 Simulation Results .............................................................................................42 4.3.1 General Features of the Leakage Field .......................................................42 4.3.2 Comparison of Metal Surface Current........................................................45 4.3.3 Standing Wave inside the Aperture ............................................................47 4.3.4 Metal Edge Effects......................................................................................49 4.4 Summary............................................................................................................50 Chapter 5 ....................................................................................................................52 Simulating Indoor Communication with ADI-FDTD Method ..............................52 5.1 Introduction........................................................................................................52 5.2 Simulation Configuration...................................................................................53 5.3 Simulation Results .............................................................................................55 5.3.1 Average Power Map ...................................................................................55 5.3.2 Line of Sight (LOS) Approximation Assessment.......................................58 5.4 Summary............................................................................................................60 Chapter 6 ....................................................................................................................61 Summary and Future Work .....................................................................................61 6.1 Summary............................................................................................................61 6.2 Future Work.......................................................................................................62 v Appendix A:................................................................................................................64 Treatment of the Current Term in the ADI-FDTD Scheme..................................64 Appendix B:................................................................................................................71 Guided Wave Propagation in Two Layer Structure ..............................................71 References:..................................................................................................................74 1 Chapter 1 Introduction 1.1 Challenges for Modern Electromagnetic Modeling The advent of radio frequency (RF) electronic devices and denser, larger chip and higher clock rates brings new challenges to modern electromagnetic (EM) modeling. In these applications, the electromagnetic wave frequency is often around one to a hundred GHz, and accurate electromagnetic field details are often desired in the critical regions of small size, i.e. of um or submicron scale. For example, on-chip interconnects often appear in the form of Metal-Insulator- Silicon-Substrate (MISS) structure. Accurate modeling of interconnect effects such as losses, dispersion, and substrate current noise is required to reduce the degraded of the performance of the circuits. But, the insulator layer is of micron scale, and metal skin- current is of submicron scale, which is much smaller compared to the wavelength (around cm scale) of the electromagnetic wave propagating on the chip. How to couple large EM wavelength (mm to cm) scale with fine material structure (of um scale) in the same simulation has become a major barrier in the development of high-speed digital and analog IC?s. Another example is the efficient modeling of electromagnetic wave shielding effectiveness and scattering which requires modeling of the EM field distribution in the close proximity and within natural and synthetic conductive material. At 2 microwave frequencies, the skin depth of highly conductive materials is in the micrometer range and the wave length is in the mm to cm range. How can we resolve these two very different spatial scales in the same simulation? Modeling of these problems requires the solution of the electromagnetic-field problem and often a time-domain Maxwell solver is used. Conventional Maxwell solvers typically use the explicit Finite-Difference-Time-Domain (FDTD) method. These explicit FDTD methods are limited by the Courant condition. Namely, the simulation time step 222 111 1 zyx c t  +  +  < , (1.1) where x, y, z is the discrete grid size in three dimensions, respectively. Table 1.1 l isted several classes of problems and the required simulation time steps [Taflove and Problem Class Required Spatial Resolution Total Simulation Time Simulation Time Step given by Courant?s Condition Required Simulation Steps Low-frequency bio-electromagnetic ~1 mm ~100 ms ~2 ps ~5?10 10 Operation of VLSI digital logic; On-chip Interconnects ~0.1 um ~1 ns (GHz) ~0.3 fs ~3?10 6 Metal skin depth; Polymer coating; Film; Membrane ~0.1 um ~1 ns (GHz) ~0.3 fs ~3?10 6 Indoor Communication ~1 cm ~1 us ~30 ps ~3?10 4 Table 1.1: Important classes of application requiring prohibitive simulation steps due to the Courant condition 3 Hagness, 2000]. As we can see, these are important potential applications of FDTD modeling where the Courant stability bound is much too restrictive. The newly revitalized Alternating-Direction-Implicit Finite-Difference-Time- Domain (ADI-FDTD) method [Zheng et al., 1999; Namiki, 1999] offers a unique advantage in that the time step is not constrained by the smallest space-cell size. The ADI-FDTD [Zheng et al., 1999; Namiki, 1999] has been shown to be unconditionally stable and the simulation time step is limited by the duration of the key temporal features of interest instead of the Courant limit. We briefly summarize the ADI-FDTD method. Maxwell?s equations are discretized with the electric and magnetic fields on different grids. By manipulating Maxwell?s equations, the differential equations are transformed to a system of tri- diagonal algebraic equations. Each matrix of the system corresponds to one specific dimension. The tri-diagonal systems are solved at each time step for the EM fields in 3D. In this thesis, we have developed the ADI-FDTD method and used it to model the Metal-Insulator-Semiconductor-Substrate (MISS) structure; aperture radiation and indoor communication. With this method, we are able to resolve the large scale (i.e. electromagnetic wave propagation) and fine scale (i.e. metal skin depth, substrate current, coating material) structure in the same simulation and the simulation time step is not limited by the Courant condition. The simulations allow us to calculate in detail parasitic current flow inside the substrate; propagation losses, skin-depth and dispersion of digital signals on non-ideal interconnects; detailed surface current and standing wave 4 pattern in aperture radiation problems; signal power map and propagation delay in complicated in-door communication scenarios. 1.2 Thesis Structure The thesis is organized as following. In chapter 2, we introduce numerical schemes of the ADI-FDTD method. Two simulation cases with analytical solutions have been used to validate the newly developed Maxwell solver. In Chapter 3, we first derived the analytical solutions for wave propagation along the MISS structure. The ADI-FDTD code is then applied to study digital signal propagation along on-chip interconnects and cross-talk between two adjacent metal lines. In Chapter 4, we applied the ADI-FDTD method to study the radiation from the bared aperture and coated aperture. In Chapter 5, the ADI-FDTD method is used to model indoor communication problems. The summary and future work is given in Chapter 6. 5 Chapter 2 Numerical Schemes of ADI Method and Code Validation 2.1 Numerical Schemes of ADI Method In the Alternating-Direction-Implicit (ADI) FDTD method [Namiki 1999; Zheng et al. 1999], Maxwell?s equations: EJEDHuB E t B JH t D       === ? = ? = ,, , , (2.1) are discretized on the conventional Yee?s [Yee, 1966] staggered grids with the electric fields ( ),,2/1(, kjix E + , ),2/1,(, kj iy E + , )2/1,,(, +kjiz E ) defined on the grid cell edge center and magnetic fields ( )2/1,2/1,(, ++ kjix B , )2/1,,2/1(, ++ kjiy B , ),2/1,2/1(, kjiz B ++ ) defined on the grid cell face center as shown in Figure 2.1. In principle, if initially condition 0= B  is satisfied, then 0= B  is preserved by Faraday?s law: 0=?  = = E t B t B   . (2.2) This requires the numerical scheme to conserve magnetic flux a priori. The use of Yee?s staggered grid [Yee, 1966] ensures 0= B  throughout the simulation. This can 6 illustrated as follows. With the staggered grids, the magnetic field time integration becomes {} yEE zEE t B kjiykjiz kjiykjiy kjix  = +++ +++ ++ / / )2/1,,(,)2/1,1,(, ),2/1,(,)1,2/1,(, )2/1,2/1,(, (2.3) and analogously for B y and B z . The corresponding discretized Equation 2.2 becomes when Equation (2.3) is substituted in, { } 0 )()( )()( )1,1,2/1(,)1,1,2/1(,)1,,2/1(,)1,,2/1(, ),1,2/1(,),1,2/1(,),,2/1(,),,2/1(, 1 1 1 =  ++ ++=       +       +         +         +      +      = ++++++++++ ++++++ + + +  zyx EEEE EEEE t B yx t B yx t B zx t B zx t B zy t B zysdB t kjixkjixkjixkjix kjixkjixkjixkjix k z k z j y j y i x i x cell    ,(2.4) i,j,k x y z Ez Ex Ey Hy Hx Hz Figure 2.1: Yee?s staggered grid: magnetic field on the grid face center; electric field on the grid edge center. 7 and thus the combined magnetic flux through all 6 cell faces remains unchanged (or constant), during the time integration, as required by Equation 2.2. In the ADI-FDTD method, at each time step, by manipulating Maxwell?s equations, we transform the differential equations to a system of tri-diagonal algebraic equations. Here, we give an example of discretizing the E x component (equations (2.5- 2.10)) during the two alternating steps. In this example, we assume uniform grid spacing in each direction for simplicity. In the actual implementation, the grid spacing is non-uniform in each dimension. This allows us to place enough resolution at places of interest. In step 1, the first half (B z ) of the right hand side in Equation 2.5 and the first half (E x ) of the right hand side in Equation 2.6 are treated as implicit. We substitute Equation 2.6 ( 1 ),2/1,2/1(, + ?+ n kjiz B ) back into Equation 2.5 and obtain the tri-diagonal Equation 2.7. For the other two dimensions (E y and E z ), we perform similar manipulation to form a system of tri-diagonal equations, which can be easily solved with the tri-diagonal matrix solver. The magnetic field 1 ),2/1,2/1(, + ?+ n kjiz B is updated using Equation 2.6 with the newly calculated 1 ),,2/1(, + + n kjix E and 1 ),1,2/1(, + + n kjix E . In the next step, we treat the other half (B y ) implicit in Equation 2.8 and (E x ) implicit in Equation 2.9. We obtain a tri-diagonal system in Equation 2.10 for E x . Similarly, we can obtain the other two tri-diagonal systems for E y and E z and solve them to update the electric field. The magnetic field is updated with equations similar to Equation 2.9. These two steps are alternated thereafter. 8 STEP 1 (for E x and B z component): 1 ),,2/1(, )2/1,,2/1(,)2/1,,2/1(, 1 ),2/1,2/1(, 1 ),2/1,2/1(,),,2/1(, 1 ),,2/1(, 1 1 + + +++ + + + +++ + +       =   n kjix n kjiy n kjiy n kjiz n kjiz n kjix n kjix E z BB u y BB ut EE   (2.5) x EE y EE t BB n kjiy n kjiy n kjix n kjix n kjiz n kjiz      =   + + + + + + + + ),2/1,,(),2/1,1,( 1 ),1,2/1,( 1 ),,2/1,( ),2/1,2/1,( 1 ),2/1,2/1,( (2.6) ) ( 1 )( 1 ) ( 1 1 2 1 1 ),2/1,(,),2/1,1(, ),2/1,(,),2/1,1(, 2 ),2/1,2/1(,),2/1,2/1(, )2/1,,2/1(, )2/1,,2/1(,),,2/1(, 2 2 2 2 2 2 1 ),1,2/1(, 1 ),,2/1(, 1 ),1,2/1(, x EE x EE y t BB y t B B z t Ed y t c t y t b y t a dEc EbEa n kjiy n kjiy n kjiy n kjiy n kjiz n kjiz n kjiy n kjiy n kjix n kjix n kjix n kjix   +      +   +   =   =  +   +=   = =? +?+? + +++ +++ + +++ + ++ + + + + ? ? ? ?  ? ? (2.7) 9 STEP 2 (for E x and B y component): 2 ),,2/1(, 2 )2/1,,2/1(, 2 )2/1,,2/1(, 1 ),2/1,2/1(, 1 ),2/1,2/1(, 2 ),,2/1(, 2 ),,2/1(, 1 1 + + + + + ++ + + + ++ + + + +       =   n kjix n kjiy n kjiy n kjiz n kjiz n kjix n kjix E z BB y BB t EE  ? ?  (2.8) z EE x EE t BB n kjix n kjix n kjiz n kjiz n kjiz n kjiy      =   + + + ++ + + + ++ + ++ + ++ 2 ),,2/1(, 2 )1,,2/1(, 1 )2/1,,(, 1 )2/1,,1(, 2 )2/1,,2/1(, 2 )2/1,,2/1(, (2.9) ) ( 1 )( 1 ) ( 1 1 2 1 1 1 )2/1,,(, 1 )2/1,,1(, 1 )2/1,,(, 1 )2/1,,1(, 2 1 )2/1,,2/1(, 1 )2/1,,2/1(, 1 ),2/1,2/1(, 1 ),2/1,2/1(,),,2/1(, 2 2 2 2 2 2 2 )1,,2/1(, 2 ),,2/1(, 2 )1,,2/1(, x EE x EE z t BB z t B B y t Ed z t c t z t b z t a dEc EbEa n kjiz n kjiz n kjiz n kjiz n kjiy n kjiy n kjiz n kjiz n kjix n kjix n kjix n kjix        +      +=   =  +   +=   = =? +?+? +  + + + + + ++ + + + ++ + + + +++ + ++ + + + + ? ? ? ?  ? ? (2.10) 10 Namiki [1999] and Zheng et al.[1999] show that the 3D-ADI method for solving Maxwell?s equation is unconditionally stable and the simulation time step is not limited by Courant?s condition. As noted by Zheng et al. [1999], the numerical dispersion can be minimized as long as the simulation time step is taken small enough to properly resolve key temporal features of the modeled electromagnetic field waveform, or the period of the highest frequency spectral component of interest. A good rule of thumb for unconditionally stable ADI algorithms is that the simulation time step should be selected to be 1/25 th to 1/20 th of the duration of the key temporal features. Mur?s [Mur, 1981] first order absorption boundary condition (ABC) is applied to simulate free space or an unbounded region. Namely, the wave propagating toward the boundary is convected away with the speed of light, for example, 0 1 ,, = ? x E t E zyzy ? is used along the x direction. ?+? and ??? sign can be chosen according to the wave propagating direction. Since discretizing the Mur?s first order boundary condition only involves two adjacent grids at the boundary along a given normal direction, the implicitness of the ADI scheme is maintained. On the other hand, with Mur?s first order absorption boundary condition, there can still be some residual reflection from the boundary. We have extended the outer boundary far enough to minimize this reflection. Other better absorption boundary conditions such as Higdon?s ABC [Taflove and Hagness, 2000], Perfect Matching Layer (PML) 11 [Berenger 1994] and the Complementary Operator Method (COM) [Ramahi 1997, 2002] are still being developed and coupled to the ADI scheme. The basic idea of Perfect Matching Layer ABC is to terminate the outer boundary of the space lattice in an absorbing material medium. In Berenger [1994]?s split-field approach, the plane waves of arbitrary incidence, polarization, and frequency are matched at the boundary. In PML layer, loss parameters are chosen and are consistent with the dispersionless medium. In the discrete FDTD lattice, spatial scaling of the PML parameters was proposed to reduce the discretization errors at material interfaces. The basic idea in the Complementary Operator Method [Ramahi, 1997] is that first-order reflections of numerical waves from the outer boundary of the FDTD grids can be cancelled. By applying the radiation boundary operators that are complementary to each other, we can obtain two independent solutions to the modeling problem. The averaging of these two solutions cancels the reflection. Later on, Ramahi [2002] developed the Concurrent- Complementary Operator Method (C- COM) to reduce the computational burden introduced in COM for calculating the solution in the entire domain twice. In C-COM (suitable for either 2- or 3-dimension implementation), only the field in the boundary layer needs to be updated independently twice and the main domain is updated within a single FDTD run. 12 2.2 Model Validation 2.2.1 2D Guided Electromagnetic Wave To test the code we applied it to a standard metal skin depth problem with a known analytical solution. We performed a 2D simulation of electromagnetic wave propagation under a metal strip of conductivity = 3.9?10 7 S/m. Fig. 2.2a shows the configuration of the simulation. The domain bottom is bounded with Perfect Electric Conductor (PEC) .The smallest grid size of 0.1 um is placed inside the metal. The grid along Z direction is of uniform size= 150 um. The Courant condition requires t < 0.33?10 -15 sec. Our simulation time step t =2?10 -13 sec. The excitation frequency = 50 GHz, corresponding to wavelength= 6 mm. The skin current Jz inside the metal has the analytical solution )/exp()/cos(  yyzkJz z ?+ , (2.11)  is the skin depth = )/(2 ? and k z is the wave number along the guide. Inside the metal, the wave is damped in the Y direction and grazes along the Z direction. Fig. 2.2b shows the comparison between the simulation and analytical calculation for current Jz inside the metal. Note that the Y axis unit is 0.1um; the Z axis unit is mm. The agreement is excellent. With ADI method, we are able to reveal the grazing wave pattern inside the metal with very large time step compared to the explicit scheme. 13 2.2.2 2D Electromagnetic Wave Scattering Since the (ADI-FDTD) method offers a unique advantage in that the time step is not constrained by the smallest space-cell size (i.e., the Courant limit), we can test it (a) X 10 -9 X 10 -9 (b) Figure 2.2: Model Verification: shows excellent agreement between numerical and analytical result. (a) Geometry. Source at z = 0; f= 50 GHz. (b) pattern of the current Jz inside the metal obtained from the simulation (top) and analytical calculation (bottom). The metal strip conductivity= 3.9?10 7 S/m. Note the Y axis unit is 0.1 um. 14 to simulate the classical 2D EM wave scattering from a highly conductive material. Figure 2.3 shows the simulation configuration. The nominal grid size is 5 mm and a fine resolution of 0.1 um is placed on the metal surface to resolve the metal skin depth. A Gaussian current Jz pulse is excited with bandwidth = 20 GHz. The simulation time step t = 2x10 -13 sec and Courant limit requires t < 0.33?10 -15 sec. Figure 2.4 shows a snap-shot of the electromagnetic field and metal skin current distribution at t = 184.4 ps. We can clearly see that the simulation resolves detailed structure of the induced skin current on the metal surface within 2 um. In order to check whether the distribution of the skin current is accurate, we monitor the electric field at a series of points (0.1 um apart) inside the metal and use Fourier transform to derive the frequency response at these locations. The relative magnitude of the electric field when traveling into the metal for several given frequencies are plotted in Figure 2.5 together with the expected analytical solution ( )/exp( yEz  , )/(2 ? = ). The agreement is good and the ADI scheme Metal 100 Grids (x0.5 mm) Current Source Y X Jz 100 Grids (x0.5 mm) 20 Grids (x0.1 um) Skin Layer Figure 2.3: Top view of the configuration for electromagnetic wave scattering from highly conductive material. 15 successfully resolves the metal skin depth. We note that in Appendix A, we use this example again to show that the treatment of the current term in equations (2.5 and 2.8) fully-implicit or semi-implicit can result in different skin depth. In general, treating the current term in Equation 2.5 fully-implicit resolves the metal skin depth correctly. Figure 2.4: Top 2 panels show the electromagnetic field distribution when the EM wave is scattered the metal. Bottom panel shows the skin current Jz inside the metal surface. Note the horizontal axis in the bottom panel is of unit 0.1 um. 16 2.3 Summary In ADI, the simulation time step and grid size is limited by the signal propagation property to avoid dispersion (20 to 25 grids per wavelength), not the smallest grid size. With an ADI method, at microwave frequency, we were able to resolve fine structure such as metal skin depth and model the electromagnetic field distribution in the close proximity and within the material. This will be further illustrated with several applications of the ADI method in the following chapters. f =10 GHz f =40 GHz Figure 2.5: Relative electric filed vs. positions into the metal surface for 10, 20, and 40 GHz. The solid lines are the analytical solutions and the colored plus signs are the simulation results. 17 Chapter 3 Simulating on-chip Interconnects with ADI-FDTD Method 3.1 Introduction Radio Frequency (RF) effects are a major factor in limiting integrated circuit (IC) performance. The complex IC interconnect structure forms a network of coupled transmission lines. Parasitic coupling between these network elements forms significant barriers in the development of high-speed digital and analog IC?s. Accurate modeling of modern on-chip interconnects (including coupling and losses) usually requires a full-wave solution to Maxwell?s equations. However, such a solution is difficult because the wavelengths of interest are much larger than the fine topological structure of IC?s. Wavelengths of the EM wave at radio frequency are typically on the mm to cm scale, while chip structures are on the micron scale, e.g. metal strip, SiO 2 layer thickness and substrate with finite conductivity. In addition, digital and mixed (broad band) signal applications require analysis in the time domain. Conventional Maxwell solvers typically use the explicit Finite-Difference- Time-Domain (FDTD) method and are limited by the Courant condition, which requires prohibitively small time steps to resolve fine structure on the submicron scale. To overcome this problem, we have applied the Alternating-Direction-Implicit (ADI- 18 FDTD) method [Namiki 1999, Zheng et al. 1999] to solve the Maxwell?s Equation in ICs, and have overcome the Courant limit. In this chapter, we use the ADI-FDTD method to model the digital pulse propagation along Metal-Insulator-Semiconductor- Substrate (MISS) structure with single and double metal strips. The simulations allow us to calculate in detail parasitic current flow inside the substrate; propagation losses, metal skin-depth and dispersion of digital signals on non-ideal on-chip interconnects. 3.2 Three Fundamental Propagation Modes for MISS Structure Guckel et al. [1967] and Hasegawa et al. [1971] investigated the properties of a microstrip line on a Si-SiO 2 system (essentially the MISS structure) and determined three fundamental propagation modes for the MISS structure, namely, Dielectric Quasi-TEM Mode, Skin-Effect Mode, and Slow-Wave Mode. Later, Shibata and Sano [1990] and Groteluschen et al. [1994] simulated the MISS structure and provided details about the properties of these three modes. In this section, we revisit the derivations of the modes in a simple 2D MISS structure and acquaint ourselves about the fundamental propagation properties. Figure 3.1 shows the side view of the MISS structure to be analyzed. To derive analytical solutions, we assume that the structure extends to infinity in the X direction. Also, inside the metal layer and ground plane, a perfect electric conductivity is assumed and the tangential electric field on these boundaries is zero. The wave propagation along the MISS is reduced to the propagation along a two-layered (SiO 2 layer and silicon substrate) structure with the 19 proper boundary condition, i.e. zero Ez at the top and bottom metal. We define b 1 and b 2 as the thickness,  1 and  2 as the relative permittivity, ? 1 and ? 2 as the relative permeability,  1 and  2 as the conductivity of the SiO 2 layer and silicon substrate, respectively. Here,  1 =4.5, ? 1 =1, and  1 =0 (ohm*m) -1 for the SiO 2 layer;  2 =12 and ? 2 =1 for the silicon substrate. The fundamental propagation mode of the waveguide shown in Figure 3.1 is a TM wave and only Ey, Hx, Ez components are involved. Here, ))(exp(,, xktjEHE zxy     . (3.1) We define  1 and  2 (complex number) as the transverse propagation constants (in the Y direction) in the SiO 2 layer and substrate, respectively. Also, we denote  as the Z Y Metal Layer Silicon Substrate SiO 2 Ground Plane b1 b2 Figure 3.1: Side view of the MISS structure. Z is the direction of propagation. To derive analytical solutions, we assume that the structure extends to infinite in the X direction. b1 and b2 are the SiO 2 layer and silicon substrate thickness. 20 longitudinal propagation constant (in Z direction and the same constant for both SiO 2 layer and substrate). In other words, the following relations are satisfied, .;2)1,(i zyii jkjk ===  (3.2) Use the boundary conditions at the top, bottom metal layer and the interface between the SiO 2 layer and substrate, we can derive three eigenvalue equations for  1 ,  2 and  (all are complex) [Guckel et al. 1967]: ' 11 2 0 22 1 ? k=+ (3.3) ' 22 2 0 22 2 ? k=+ (3.4) 0)tanh()tanh( 22 2 2 11 1 1 '' =+ bb       . (3.5) The primed relative permittivities are defined by )/( 0 '  j iii += , i = 1, 2; ck / 000 ? == , where  0 and ? 0 are the permittivity and permeability in the vacuum, respectively. See Appendix B for a detailed derivation of Equations (3.3-3.4). The field components in each layer can be calculated as: z ii ii z i yi e b by EE       = ]sinh[ )](cosh[ 0   z ii ii z i i xi e b by E j H       = ]sinh[ )](cosh[ 0 0 '   z ii ii zzi e b by EE     = ]sinh[ )](sinh[ 0   . (3.6) 21 Here, i = 1, 2 and ?-? sign for i = 1, ?+? sign for i = 2; i =1 represents the field distribution within SiO 2 layer and i =2 represents the field distribution within the substrate layer. ?y? is defined to be 0 at the SiO 2 and substrate interface and positive y is upward. E z0 is the value at z = 0, and y = 0. We can further write ? jj effeff +== ** for the wave propagation constant along the z direction. Here, * eff  and * eff ? are the effective complex permittivity and permeability of the double layer.  is the attenuation factor and  is the propagation factor (or )//( c is the so-called slow wave factor) along the z propagation direction. According to Hasegawa et al. [1971], * eff  and * eff ? can be calculated in the following way once  1 ,  2 and  are given:   =  = ==   ! " # == 2 1 *'''* 1 2 1 0 ' '''* ]tanh[ 11 ]tanh[ 111 i ii i ieffeffeff i ii ii effeffeff b b j b b j   ????    , (3.7) where 21 bbb += , total thickness of the two layer. The characteristic impedance Z 0 can be approximated as * * 0 eff eff a b Z  ? = , where ?a? is the width of the metal strip if the assumption of infinite extension in X is removed. So, the question is how to calculate  1 ,  2 and  from Equations (3.3-3.5). After fixing the geometry of the MISS structure, i.e. b 1 and b 2 are given, we need to 22 calculate  1 ,  2 and  for different substrate conductivity  2 and frequency f. We can rewrite Equations (3.3-3.5) as 0)tanh()tanh( 0 22 2 2 11 1 1 ' 22 2 0 ' 11 2 0 2 2 2 1 '' =+ =+ bb kk       ?? . (3.8) This set of coupled non-linear equations can be solved using the Newton-Raphson method for given initial guesses. Care needs to be taken since both the substrate conductivity and wave frequency can vary by several orders of magnitude. Once  1 ,  2 are solved,  can be obtained using either Equation 3.3 or 3.4. 23 As an example, for a MISS with SiO 2 thickness b 1 = 2 um and substrate thickness b 2 = 200 um, we solved Equation (3.8) and the calculated attenuation factor and slow wave factor are shown in Figures 3.2 and 3.3. 1.5??? f e 4.0??? f  Skin-Effect Mode Slow-Wave Mode Dielectric Quasi-TEM Mode 0.3??? f 0 Figure 3.2: Color-map of log of attenuation factor  (along propagation direction Z) vs. substrate resistivity =1/ 2 and wave frequency for MISS structure with SiO 2 layer thickness = 2 um and substrate thickness = 200um. The three black lines divide the map into 3 regions of fundamental modes. 24 In Figure 3.3, the slow wave factor ph Vcc /)//( = essentially indicates the relative ratio between the speed of light and the wave phase velocity. The larger the slow wave factor, the slower the wave propagates. In Figure 3.3, 2 220 2 2 1 b f ?$  = (characteristic frequency for skin-effect in silicon substrate); 02 2 2 1   $ = e f (dielectric relaxation frequency in silicon substrate). Another two characteristic frequencies can Skin-Effect Mode Dielectric Quasi-TEM Mode Slow-Wave Mode 4.0??? f  1.5??? f e 0.3??? f 0 Figure 3.3: Color-map of slow wave factor )//( c (along propagation direction Z) vs. substrate resistivity =1/ 2 and wave frequency for MISS structure with SiO 2 layer thickness = 2 um and substrate thickness = 200um. The three black lines divide the map into 3 regions of fundamental modes. 25 also be defined: 2 1 01 2 2 1 b b f s   $ = (relaxation frequency of interfacial polarization); 1 11 0 3 2     ! " # +=  fff s (characteristic frequency of slow-wave mode). From Figure 3.3, we can see that there exist three mode regions divided by the three solid black lines. As noted in Hasegawa et al. [1971] and Groteluschen et al. [1994], these three modes can be characterized as follows. 1) Dielectric Quasi-TEM Mode: This mode appears when e ff 5.1% or the product of f and  2 is large and the substrate conductivity is low. The substrate acts like a dielectric at high frequencies. In this mode, the real part of * eff  : & = seff  0 ' and real part of * eff ? : 0 ' ?? = eff , where [ ] 1 2211 /)//(  & += bbb s  . The attenuation factor can be approximated as 00 ?$ & = se f . The propagation factor c s / & = . Both the slow wave factor and the attenuation factor are small. In the practical situation ( 2112 bb  << ), the propagation speed is nearly equal to 2 / c and the dispersion is small. 2) Skin-Effect Mode 26 This mode appears when  ff 4% and  fff s /01.0 2 ) or the product of f and  2 is large and the conductivity  2 is large. Note, 2 220 2 2 1 b f ?$  = . With increasing substrate conductivity, the skin-effect arises in the substrate and the return current is enhanced in the substrate. In this mode, 00 ' seff  = and ) 2 ( 1 10 '  ?? += b b eff , where )/(2 02 ? = is the skin depth of the substrate and )/ ( 110 bb s  = . In most cases, 2 1  < 20. The propagation velocity thus slows down. 27 From Figure 3.2, we can further see that considerable losses appear at high frequencies: in the transition region between the slow-wave mode and the dielectric quasi-TEM mode (mainly due to loss in the substrate transversal electric field) and in the upper part (high frequency) of skin-effect mode (mainly due to longitudinal substrate skin current). On the other hand, Figure 3.3 suggests that the wave propagation speeds in both the dielectric quasi-TEM mode and the skin-effect mode are high. We also note that Equations 3.3-3.5 are not limited in application only to the case b 1 << b 2 . As an example, we can revisit the problem of 2D guided wave under metal strip presented in section 2.2.1. In this case, the top layer is highly conductive metal and the bottom layer is vacuum. The solution to Equations 3.3-3.5 yields the skin current Jz inside the metal )/exp()/cos(  yyzkJz z ?+ , (3.9)  is the skin depth = )/(2 ? and k z is the wave number along the guide. 3.3 Single Metal-strip-Insulator-Semiconductor-Substrate (MISS) Structure In real on-chip scenarios, the metal strip is of finite width and metal edge and skin depth effects need to be considered. To model the detailed propagation of broadband signal and the induced substrate current in such a complex structure, we need to resort to the FDTD Maxwell solver. Here, we have applied the 3D ADI- FDTD code to study the digital pulse propagation along a Metal-Insulator- 28 Semiconductor-Substrate (MISS) structure. Fig. 3.4 shows a cross-section of the interconnect MISS structure we simulated with our ADI code. Table 1 lists the geometry and grid configuration for the MISS structure. Along the Z, the direction of wave propagation, we have 100 grids of uniform size = 25 um. In XY cross-section we have a non-uniform mesh with the finest grid spacing = 0.1 um placed in the metal strip (6 um?1.8 um). Figure 3.4: Cross section of the simulated MISS structure. The XY cross section of the metal strip is 6um x 1.8 um and thickness of the SiO 2 layer is 2 um. Z is the direction of propagation and lumped current flow. Absorption boundary conditions are applied in the vacuum. Region (along Y) a (Vaccum) b (Metal) c (SiO 2 ) D (Substrate) Length 505 um 1.8 um 2 um 500 um # of grids 33 18 10 27 Region (along X) e F f1, f2 G Length 555 um 4 um 1um ea. 555 um # of grids 28 5 10 ea. 28 Table 3.1: Geometry and grid configuration for the MISS structure. Refer to Figure 3.4 for region labels (a-g). 29 The simulation time step is t = 1?10 -13 sec and the Courant condition requires t < 0.33?10 -15 sec. The simulation has 83?89?100 space grid points and is run for 1000 time steps. Mur?s first order absorption boundary condition is applied at the top, front, back and two sides of the structure. PEC boundary condition is applied at the bottom. The excitation of Ex and Ey field in the source plane is a solution of Poisson?s equation and enveloped with a transient profile, e.g. Gaussian or digital function. The simulation time is about 3-4 hours on a single PC (2.4 GHz Pentium 4). Figs. 3.5-3.7 shows simulation results for a fast 1V, 20psec digital pulse of rise-time = 2ps, excited at one end of the interconnect. The metal strip conductivity = 5.8?10 7 S/m, typical for copper. The substrate doping is set to be n = 10 17 /cm 3 , which corresponds to substrate conductivity = 2260 S/m, resistivity  = 4.4?10 -4 ohm?m. The bandwidth of the signal is around 50 GHz. From section 3.1, since the substrate resistivity is low, e f =3385 GHz;  f = 0.448 GHz and  ff *4> , we can say the propagation mode is in the skin-effect region. Fig. 3.5 shows the evolution of the voltage signal at Z = 500, 1000 um along the wave propagating direction and the signal amplitude is lowered and broadened. The results show digital signal losses and dispersion. The higher frequency components of the signal suffer larger damping. These losses occur mainly in the 30 substrate. Figure 3.6a shows the cross section of the electric field Ey (perpendicular to the metal and substrate surface) at Z= 1000 um and t = 50 ps. The Ey field concentrates inside the SiO 2 layer. This corresponds to the skin-depth mode propagation along the MISS structure [Hasegawa et al. 1971; Shibata and Sano 1990]. The electromagnetic field is guided along the channel formed by the metal skin current and substrate current. Figures 3.6b shows the XY-cross section of current Jz inside the metal (along the direction of the signal propagation), at Z = 1000 um and t = 37 ps. We see that due to the skin depth effect, the current concentrates near the metal edges and surfaces and decays into the center of the metal, giving rise to resistive losses. With the ADI-FDTD Figure 3.5: Voltage observed at different Z locations (z=500 um, 1000 um) along the MISS strip. 31 method, we are able to reveal the detailed structure of electric field inside the metal and SiO 2 layer. Figure 3.7a shows the top view of substrate current Jz at 5 um below the SiO 2 layer for t = 37 ps. The blue and red shaded areas correspond to the rising and falling of the signal, showing the spread of the current almost a tenth of a mm away from the interconnect edge, which can obviously lead to parasitic cross talk. Fig. 3.7b shows X (um) X 10 5 SiO 2 metal (a) (b) Figure 3.6: (a) Cross section of Ey (V/m) field at Z= 1mm and t= 50 ps. The electric field Ey concentrates in the SiO 2 layer. (b) Cross section of current Jz (mA/um 2 ) inside the metal at Z=1mm and t = 37 ps. Shows metal skin-depth effect losses. 32 side view giving the depth of the substrate current. The substrate skin depth is around tens micron. The substrate current contributes most to the parasitic interconnect losses. From Figure 3.6b and Figure 3.7, a rough estimation yields the comparable magnitude of traveling and returning current in the metal strip and substrate, respectively. x10 -6 (a) x10 -6 (b) Figure 3.7: (a) Top view of current Jz (uA/um 2 ) in the substrate (5um below SiO 2 layer) at t = 37 ps. Blue and red shaded areas correspond to the rising and falling of the signal. Shows potential interference and coupling in lateral direction. (b) Side view of current Jz (uA/um 2 ) in substrate at X = 0. Shows current penetration to the substrate. Substrate doping n = 10 17 /cm 3 . 33 We also investigate the effects of the substrate doping on the digital signal propagation along the MISS structure. Fig. 3.8 shows the voltage at different Z locations for substrate doping changing from n = 10 16 /cm 3 ( = 226 S/m) and 10 18 /cm 3 (( = 22600 S/m). The voltage is obtained by integrating the electric field between the metal strip and the ground plane. For n = 10 16 /cm 3 , from section 3.1, since the substrate resistivity is low, e f =338 5GHz;  f = 4.48 GHz and  fGHzf *4)50(~ > , we can say the propagation mode is in the skin-effect region. Similarly, for n = 10 18 /cm 3 , e f =33850GHz;  f = 0.0448 GHz and  fGHzf *4)50(~ > , the propagation mode is also in the skin-effect region. From Figure 3.8, the electromagnetic waves, all of which are in the skin-depth mode region [Hasegawa et al. 1971; Shibata and Sano 1990], suffer different dispersion and Figure 3.8: Voltage observed at different Z locations along the MISS strip with substrate doping n1 = 10 18 , and n2 = 10 16 /cm 3 . 34 dissipation for different substrate dopings. Lower substrate doping (lower conductivity or higher resistivity) in this region yields more losses and dispersion, which can be inferred from Figure 3.8. 3.4 Double Metal-strip-Insulator-Semiconductor-Substrate (MISS) Structure We also apply the ADI code to simulate the cross-talk effects between two parallel metal strips on a MISS structure as shown in Fig. 3.9. The passive metal line (left) is grounded. A digital voltage wave is applied on the active line (right). The spacing between the two metal strips is 20 um. The metal strip conductivity = 5.8?10 7 S/m, typical for IC interconnects. The substrate doping is set to be n= 10 17 /cm 3 . Again, along the Z, the direction of wave propagation, we have 100 grids of uniform size = 25 um. In the XY cross-section we have a non-uniform mesh with the finest grid spacing = 0.1um; and the 500 um 500 um 1.8 um 2 um x y z 555 um555 um 6 um 6 um 20 um Lossy Silicon Substrate Active metal line SiO2 VacuumPassive metal line Figure 3.9: Cross section of the simulated double-line MISS structure. The XY cross section of the metal strip is 6um x 1.8 um; the SiO 2 layer is 2 um thick. Substrate doping n= 10 17 /cm 3 . 35 Figure 3.10: Voltage observed at different Z locations (z=500 um, 1000 um) along the active (line 1) and passive (line 2) MISS strips. simulation time step is chosen to be t =1?10 -13 sec and the Courant condition requires t < 0.33?10 -15 sec. Fig. 3.10 shows simulation results for a fast 1V, 20psec digital pulse of rise-time= 2ps, excited at one end of the right-side interconnect. We can see that the voltage is induced on the passive line when the active line experiences a fast voltage change. The figure (solid line) shows the voltage signal at Z= 500, 1000 um and the signal amplitude is lowered and broadened along the active line. Also shown is the magnitude of the induced voltage on the parallel line. We see that coupling gives rise to voltage levels of as much as 0.2V, even though the lines are separated by 20um. To understand the coupling between the metal lines, we plot the substrate current Jz in Fig. 3.11. Figure 3.11a shows the top view of substrate current Jz at 5um below the SiO 2 layer for t = 30 ps. The red and blue shaded areas correspond to the 36 rising and falling of the signal. Metal strips are centered at X = -13 um (passive ) and X= 13 um (active) and extend along Z. We can see that the substrate current clearly spreads to the neighboring interconnect, away from the interconnect edge, which lead to parasitic crosstalk. The passive line shield the further spreading of the substrate (a) (b) Figure 3.11: (a) Top view of current Jz in the substrate (5um below SiO 2 layer) at t = 30 ps. (b) Substrate current Jz in the XY plane. (Note: the units on the color scale in Figs. 3.8a and 3.8b range from 3 x 10 -6 A/um 2 (red) to -3 x 10 -6 A/um 2 (blue)). 37 current in the lateral direction. The induced EM wave propagates in the channel formed by the passive metal strip and the substrate. Fig. 3.11b shows a XY cross-section view giving the depth of the substrate current. The substrate current is concentrated under the active metal strip and spreads to the region under the passive line. The substrate current depth is around tens um. 3.5 Summary In this chapter, we show that the ADI method is efficient in simulating signal propagation along on-chip interconnects. We are able to study both the electromagnetic wave propagation and the detailed metal skin depth and substrate current effects on single MISS and double MISS structure without being limited by the Courant condition. It is worth noting that we assume the substrate conductivity or equivalently doping is constant in depth in our study. In reality, the spatial concentration of the free charge (electrons or holes) can vary when the metal strip is biased. This opens a way to control the signal propagation property by intentionally controlling the doping and biasing voltage. To model this requires incorporating the device level physics to replace the simple relation EJ  = in the Maxwell equation. See Chapter 6 for a discussion on future work. The future challenges will also be to simulate interconnects with spatially varying substrate doping and active load. 38 Chapter 4 Simulating Electromagnetic Field Radiation from Aperture with ADI-FDTD Method 4.1 Introduction Metal enclosures with apertures are typically used as chasses for high-speed computer systems. Apertures present on these enclosures are primarily intended for thermal management (i.e. ventilation) and cables. As the clock frequency increases, the metal wire and interconnects on board can act as antennas or radiators, whose radiation intensity increases with frequency. The electromagnetic field radiation originated from the circuit board can leak undesirably to the environment through the apertures on these enclosures. Since the wavelength of clock frequencies starts to approach the size of the aperture, the leakage of electromagnetic waves through apertures in enclosures has become a critical electromagnetic compatibility (EMC) and electromagnetic interference (EMI) problem. Li and Ramahi [IEEE APS/URSI Symposium, 2002] proposed a novel structure to reduce field leakage from apertures by coating a layer of conducting material (of much smaller conductivity) such as lossy polymer on top of the metallic aperture. The coating works to reduce radiation, especially at frequencies at or close to the natural resonance of the aperture. In this 39 chapter, we applied the ADI-FDTD method that is capable of providing a more accurate definition of the electromagnetic fields within the very thin coating layer surrounding the aperture. With ADI-FDTD method, we can resolve fine structure. In what follows, we will investigate the effect of thin-film material inclusion on the reduction of the field in the exterior of an enclosure containing the aperture. We will also present results showing current distribution on the material surrounding the aperture, and present a discussion on the physical aspects of the aperture radiation phenomenon. 4.2 Aperture Model and Simulation Configuration The simulation configuration for metal aperture with and without coating layer is indicated in Figure 4.1a and 4.1b. The metal plate is of thickness 2mm. A Gaussian current source of bandwidth = 20 GHz is excited at 20 cm to the left of the metal plate. The field source is polarized in X direction for maximum aperture radiation. The observational point is placed at 20 cm to the right of the metal plate. The metal aperture dimension is 2cm (W) ? 2mm (H) ? 2mm (Thick). In Figure 4.1b, we show another configuration: the metal aperture coated with thin layer of conductive material (thickness = 2 mm; width = 5mm) surrounding the region of the aperture on both side. The coating layer has conductivity = 5 mhos, typical for conductive polymer,  r = 4 and u r = 4. The metal plate is extended and terminated into surrounding absorption boundary to isolate the effect of the aperture. We note that the area of the aperture 40 remains unchanged in the scenario with coating, hence, the heat transfer performance of the metal plate is not compromised. Figure 4.2 shows the simulation grid configuration in the YZ plane (perpendicular to the metal plate) together with an amplified view of the grid configuration around the surface of the metal. The largest grid size = 0.75 mm is placed in the vacuum, and the smallest grid size = 0.1 um is placed on both side of the Z Y X Metal Plate 2 mm thick Gaussian current source Aperture Observation point (a) (b) Figure 4.1: Simulation configuration for scenario (a) Bare metal aperture of 2cm(W)?2mm(H)? 2mm(thick). The radiation source and observation point (at 20 mm to the right of the metal) are also indicated on the figure. (b) Aperture with thin layer of coating on both sides. Coating thickness = 2 mm; width = 5mm.The coating conductivity is 5 mhos. 41 metal surface. We used 166?178 ?85 grids in 3D. The simulation time step t = 5?10 - 13 sec and Courant?s limit requires t < 0.33?10 -15 sec. To simulate 2400 steps on a single PC requires about 10 hours. Amplified View Figure 4.2: Top panel: simulation grid configuration in the YZ plane (perpendicular to the metal plate). Bottom panel: amplified view of the grid structure on top of the metal surface. The dense grids are used to resolve the skin depth and reveal metal surface current. 42 4.3 Simulation Results 4.3.1 General Features of the Leakage Field We first present the general features of the leakage field from the simulations with and without the coating layer around the aperture. Here on, we refer to the aperture without coating layer as Scenario A: bare slot; and the aperture with coating layer as Scenario B: with coating. Figure 4.3: The observed electric field Ex (20 mm to the right of the aperture) in Scenario A: Bare Slot and Scenario B: With Coating. 43 Figure 4.3 shows the observed electric field Ex (polarized for maximum aperture radiation) at 20 mm to the right of the metal aperture for Scenario A and B. Although the excitation Gaussian current source decays rapidly, once the electromagnetic field reaches the aperture, the aperture acts as an antenna and radiates at its resonant frequency. This can be seen from the residual electric field oscillation in Figure 4.3. We can also see the amplitude of the transient electric field in the Scenario B (with coating) is significantly smaller than the one in Scenario A (bare slot). This suggests that the coating layer increases the shielding effectiveness. In order to understand the shielding performance of the coating layer at the resonant frequency of the aperture, we transformed the transient electric field signal to the frequency domain by Fourier transformation and divided the amplitude of the excitation electric field at that frequency. The resulting relative electric field in the frequency domain is plotted in Figure 4.4. When the aperture is not loaded, the radiation is maximized at the frequencies at which the aperture becomes resonant. Given the dimension of the aperture 2cm (W) ? 2mm (H) ? 2mm (Thick), the simulation indicates the resonant frequencies are at around 7.0 GHz and 14 GHz. When the lossy layers are placed on both sides of the aperture as in Figure 4.1b, the radiation drops about 6 times at first resonant frequency around 7 GHz and 1.5 times at the second resonant frequency around 14 GHz. There are some trade-offs at the low frequency (< 6GHz) at which the radiation is increased compared to the bare slot 44 scenario. But, this is insignificant compared to the large reduction of the radiation at the aperture resonant frequencies, which is most harmful due to its impulsive nature, and the overall shielding performance over the broad spectrum (from 6 GHz to 16 GHz). Although we can determine the shielding effectiveness of the coating layer from the observed electric field, it is not clear how the efficient shielding is originated. In the following sections, we will investigate the detailed current and electric field Figure 4.4: The observed relative electric field Ex in frequency domain for Scenario A: Bare Slot and Scenario B: With Coating. The electric field at a given frequency is normalized with field amplitude propagates in vacuum. 45 structures around the aperture as revealed with ADI-FDTD method to understand the origin of the reduction of the leakage electric field in the scenario of with coating. 4.3.2 Comparison of Metal Surface Current Figure 4.5 shows snap-shots of electric field and surface current pattern at t = 189 ps for Scenario A (base slot) and Scenario B (with coating), respectively. The top row shows the electromagnetic wave reaches the metal plate and starts to re-radiate from the metal aperture. Here, the metal plate is located at y= 50 mm. The third row shows the detailed surface current Jx on the left side of the metal surface (XZ plane). We can see that the metal surface current surrounding the aperture is largely reduced when the coating layer is present. The resistive material acts as a lossy termination and suppresses the current penetrating into the metal. The fourth row in Figure 4.5 shows the detailed surface current Jx on the right side of the metal surface (XZ plane). Both the spatial distribution and the magnitude of the induced surface current are reduced. On the other hand, these surface currents serve as sources of radiation which propagate to the observational point. The reduction of the induced metal surface current results in the reduction of the leakage radiation. 46 Figure 4.5: Snap-shot of the electric field and current pattern at t = 189 ps for Scenario A (left column) and Scenario B (right column). Top two rows show the electric field Ex in the YZ plane (perpendicular to the metal plate). The second row shows only the Ex on the right YZ plane with smaller color scale. Bottom two rows show the surface current Jx on the Left and right side of the metal surface (XZ plane), respectively. Note, the color scales are different for the third and fourth rows. 47 4.3.3 Standing Wave inside the Aperture In order to understand the resonant nature of the aperture leakage radiation, we investigated the electric field distribution inside the aperture. Figure 4.6 shows the Figure 4.6: Snap-shot of the electric field and current pattern at t = 210 ps for Scenario A (left column) and Scenario B (right column). Top two rows show the electric field Ex in the YZ plane (perpendicular to the metal plate). The second row shows only the Ex on the right YZ plane with smaller color scale. The third row shows the electric field Ex inside the aperture in the YZ plane. The fourth row shows the surface current Jx on the right side of the metal surface (XZ plane). 48 electric field and current pattern at a later time t = 210 ps for Scenario A (left column) and Scenario B (right column). From the top row, we can see that the electromagnetic field is reflected from the metal plate. One interesting phenomenon to note is the electric field inside the aperture remains there. The second row shows radiation starts to emerge from the aperture. In the third row of Figure 4.6, we can see clearly the formation of standing wave inside the aperture region. This explains the remnant radiation from the aperture even after the incident electromagnetic wave is reflected back. The narrow center-fed slot is often referred to as a magnetic dipole. The dominant wave pattern inside the aperture is a symmetrical standing wave of the form |))|(sin( +lk . Here, l2 is the width of the metal aperture, + is the distance from the aperture center along z and ,$ /2=k . 2/2 ,-l is the resonant condition to maximize the radiation. We also note that if the source radiation is polarized along Z direction, the aperture width intersecting the wave front becomes 2 mm and the resonant frequency shifts to around 75 GHz of much smaller wave power (due to smaller component in the source). Another point to note is that the magnitude of the electric field oscillation inside the aperture is much smaller in the Scenario B with coating layer. The fourth row in Figure 4.6 again shows the reduction of both the spatial distribution and the magnitude of the induced surface current on the right side of the metal plate when the aperture is coated with lossy material. 49 4.3.4 Metal Edge Effects As an advantage of using ADI-FDTD method, we are able to show the detailed current pattern around the metal edges. Figure 4.6 shows the snap shots of the electric Figure 4.6: Snap-shot of the electric field and current pattern at t = 313 ps for Scenario A (left column) and Scenario B (right column). Top two rows show the electric field Ex in the YZ plane (perpendicular to the metal plate). The second row shows only the Ex on the right YZ plane with smaller color scale. The bottom two rows show the current Jx at the cross-section (YZ plane) of metal edge (within 1um in Z). Note the bottom two rows differ in color map scale. 50 field and current pattern at t = 313 ps for Scenario A (left column) and Scenario B (right column). Top two rows in Figure 4.6 show the leakage radiation from the aperture. From the second row, we can see that the magnitude of the leakage radiation field (on the right side (Y > 50 mm)) is reduced in Scenario B with coating layer. Interesting aspects of the current Jx at the cross-section (YZ plane) of metal edge (within 1um in Z) are shown in the bottom two rows. The current decays into the center of the metal within 1 um and the current concentrates at the metal edges. With ADI-FDTD, we can resolve the detailed features of the metal skin current. 4.4 Summary In this chapter, we have shown that the multi-grid ADI-FDTD algorithm can model the problem of a radiating object in the presence of an aperture coated with a thin film material. We are able to reveal fine features: metal surface current, standing wave inside the aperture, and metal edge effect. The coating works to reduce radiation, especially at frequencies at or close to the natural resonance of the aperture. In ADI, the simulation time step and grid size are limited by the signal propagation property to avoid dispersion, not the smallest grid size. Mitigation of transmitted radiation sees critical applications in the field of electromagnetic interference, enhancing radiation. It is also important in a wide range of applications from maximizing efficiency of antennas to improving the effectiveness of near field optical microscopes. ADI-FDTD 51 method provides an efficient way of modeling the interaction of electromagnetic wave and thin layered structure and can be applied to these applications. 52 Chapter 5 Simulating Indoor Communication with ADI-FDTD Method 5.1 Introduction With the advent of RF electronic devices and Personal Communication System (PCS), the interest to characterize radio signal propagation inside buildings becomes larger [Rappaport 2002], and so do the difficulties. The indoor radio channel differs greatly from the traditional outdoor radio channel. In the situation of indoor communication, the distances covered are much smaller, and the variability of the environment is much greater for a much smaller range of transmitter-receiver separation distance. The propagation within buildings suffers reflection, diffraction and scattering, and is strongly influenced by specific and variable features such as the layout of the building and the construction material. To accurately model the indoor field propagation requires taking into account of multi-path propagation scenarios. The conventional way of studying of the outdoor radio signal propagation is through ray- tracing. But when the wavelength of radio frequency wave (i.e. ,= 30 cm for 1GHz RF wave) is comparable to the spacing among metal studs, pipes, and furniture inside the building, interference and diffraction can occur. To model the detailed field power at the receiver site and take into account the EM wave interactions with conductive, semi-conductive, dielectric material, or even magnetic material during the propagation, we can use the FDTD method the solve the Maxwell Equations. In 53 particular, the ADI-FDTD method provides a promising way of modeling the interaction between EM waves and fine structure (of size much smaller than the wavelength) inside the building. In this chapter, we will show some examples of modeling the radio frequency wave propagation inside buildings and discuss the effects of metal studs inside the wal,l and wave frequency on the signal propagation. 5.2 Simulation Configuration The simulation configurations for indoor-communication are shown in Fig. 5.1a (wall with metal stud) and 5.1b (wall without metal stud). The two adjacent rooms are of 4m?4.5m each. The wall thickness is 12 cm for inner wall and 20 cm for outer wall. The cross section for the wood door is 90cm?6cm. The cross section for the metal stud is 5cm?8cm. The metal stud spacing is 30 cm. Table 5.1 lists the material property used for the simulation.  0 is the permittivity for the vacuum. The permeability is assumed to be ? 0 everywhere. Material Wall Wood Door Metal Stud Otherwise: empty Permittivity 10  0 42  0  0  0 Conductivity 0.0005 S/m 0 S/m 10 7 S/m 0 S/m Table 5.1: Material Property for the Simulation Configuration In the simulation, the grid resolution is 1.0 cm and 1100?750 grids are used in XY plane. Simulation time step is 1.0?10 -10 sec. The Courant condition requires t < 3.3?10 -11 sec. The transmitter is modeled as a current source and placed at the center of the right room. Mur?s first order absorption boundary condition is applied at four sides 54 to model the signal propagation outside the building. In our simulation, the transmitter radiates in the sinusoidal form with current source Jz  Sin(2$ft) for a period 3T, where f = 1/T. The current source is switched Source Back Front R i g h t L e f t room 1 room 2 Hall Way R i g h t L e f t (a) (b) Figure 5.1: (a) Floor map configuration in XY plane. The dark blue regions are empty space. The light blue regions are the wall. The black spots inside the light blue regions are metal studs. The yellow regions are the wood door; (b) Similar configuration as (a), but without the metal studs inside the wall. 55 off after 3T second and we wait long enough for the signal to propagate throughout the building and decay away. Two sets of simulation are performed: One with the excitation frequency = 433 MHz for floor map 5.1b (wood wall) and floor map 5.1a (wall with metal studs), respectively; another one with the excitation frequency = 2.4 GHz for floor map 5.1b (wood wall) and floor map 5.1a (wall with metal studs), respectively. 5.3 Simulation Results 5.3.1 Average Power Map For 433 MHz RF wave, the wavelength in vacuum is 69.3 cm and inside the wall is 22 cm; for 2.4 GHz, the wavelength in vacuum is 12.5 cm and inside the wall is 4 cm. Figure 5.2a and b show the average electric field power maps for excitation frequency = 433 MHz in building with pure-wood wall and metal stud-filled wall, respectively. Similarly, Figure 5.3a and b shows the average power maps for excitation frequency = 2.4 GHz in the two kinds of building. In the configuration with metal studs, the metal stud spacing is 30 cm. From Figures 5.2 and 5.3, we can see the formation of standing wave pattern. Wall plays an important role in shielding the signal since only refracted waves can penetrate the wall. In general, metal studs cause significant interference pattern. For 2.4 GHz excitation, the power pattern in the source room doesn?t decay as fast as in the 433 MHz excitation scenario. It seems the wall and metal studs shield the waves excited with 2.4 GHz better and the cavity modes have larger power 56 (a) 433MHz with wood wall (b) 433MHz with metal Stud Figure 5.2: Average Signal power (in dB) distribution for transmitted signal with frequency = 433 MHz in (a) building with wood wall; (b) building with metal studs inside the wall. 57 (a) 2.4 GHz with wood wall (b) 2.4 GHz with metal Stud Figure 5.3: Average Signal power (in dB) distribution for transmitted signal with frequency = 2.4 GHz in (a) building with wood wall; (b) building with metal studs inside the wall. 58 inside the room. Leaking of relative power for 433 Mhz excitation is much faster. 5.3.2 Line of Sight (LOS) Approximation Assessment In general, indoor channels can be classified as line of sight (LOS) or obstructed (OBS). We can use the simulation to assess how good the line of sight approximation is in an indoor situation. We placed monitoring points on 4 sides outside the wall: 0.5m to the wall, 1.0m apart. We record the arrival time of the first dip in the received signal at the monitoring point. This is illustrated in Figure 5.4. Here, we assume that the transmitter and receiver have been synchronized and the receiver knows when the transmitter sends out the signal. When the first dip arrives, we use this signal propagation time or the delay time t to estimate the distance Detect First Dip Delay for signal to subside Figure 5.4: Illustration of detecting the first dip of the signal at the receiver side to estimate the distance between the transmitter and receiver. 59 between the transmitter and receiver by using c?t, where c is the speed of light. In this way, we can assess the validity of line of sight approximation in indoor- communication. In Figure 5.5, we plot the estimated distance vs. actual distance with color-coded symbols. The blue line shows the scenario if the estimated distance equals actual distance. We can see the estimation error is different for monitoring at different sides. In general, the shorter the separation, the better is the LOS approximation. At the front and the left sides, the LOS approximation degrades which is due to the multi- (a) Wood Wall, f = 433 MHz (b) Wall with Studs, f = 433 MHz (c) Wood Wall, f = 2.4 GHz (b) Wall with Studs, f = 2.4 GHz Figure 5.5: We plot the estimated distance vs. actual distance with color-coded symbols. The color represents the monitoring points on different sides as indicated in the legend. See Figure 5.1 for the definition of the four sides. The blue line show the scenario if the estimated distance equals actual distance. 60 path propagation delay resulted from the multi-wall in-between the transmitter and receiver. The LOS approximation error with excitation frequency f = 433 MHz is smaller than the simulation with f = 2.4 GHz. This is because the much smaller wavelength suffers larger multi-path reflection. The wall with studs causes more error to the LOS approximation. Here, we further list some extreme error measures. At f = 433 MHz, wood wall: maximum LOS estimation error = 0.4m for distance = 7.6m; wall with metal studs: maximum LOS estimation error = 0.9m for distance = 8.3m. At f = 2.4 GHz: wood wall: maximum LOS estimation error = 1.0m for distance = 8.3m; wall with metal studs: maximum estimation error = 1.4m for distance= 8.3 m. As we can see the LOS estimation error can be minimized by monitoring from optimized location and using magnitude info. (Stronger signal means closer to the source; near side is with less error.) 5.4 Summary In this chapter, we presented simulations of indoor communication using the ADI-FDTD method. As we see, using ADI method, we can not only determine the indoor signal propagation channel, but also retrieve the detailed wave power map in a complicated floor structure. We also assessed the line of sight approximation used in indoor communication. We need to be aware of the limitation of the ADI-FDTD method: we need large nymber of grids to resolve the large scale, i.e. we can only simulate 2D structure in the case of indoor communication. 61 Chapter 6 Summary and Future Work 6.1 Summary In this thesis, we introduced the ADI-FDTD method for electromagnetic modeling. The code is benchmarked with metal scattering and 2D guided wave problems. We applied the ADI-FDTD Maxwell equation solver to simulate 1) digital signal propagation along on-chip interconnects and cross-talk among interconnects; (2) metal aperture radiation and the reduction of the radiation with added polymer coating layer; (3) indoor communication in different building structures and at different frequencies, and assess the validity of the line of sight signal propagation approximation. From these applications, we can see the usefulness of the ADI-FDTD method in resolving the large scale wave propagation wavelength and small scale material structure within the same solution. All of these originate from the facts that the ADI- FDTD method is unconditionally stable, and the simulation time step is not limited by the Courant condition. Of course, caution needs to be taken so that only placing fine grid at the places needed and the simulation time step needs to be much smaller than the fundamental wave propagation period. 62 6.2 Future Work Although ADI-FDTD shows promising aspects in resolving multi-scale electromagnetic problem, there is still much to be done and we listed them as following. 1) Need to couple to a better absorption boundary condition, e.g. PML or COM. 2) Need to be able to insert lumped elements into the ADI method. These lumped elements can be a resistor, capacitor or even a nonlinear device such as a MOSFET which can be approximated as a group of connected lumped elements. 3) We mentioned in Chapter 3 that a better model for the biased silicon substrate should involve device level physics. Initial work by Wang et al. [2002] solves the coupled Maxwell?s equation and drift-diffusion equation in the frequency domain for a 2D MISS cross section. ADI-FDTD provides a promising scheme in which we can solve the problem in the time domain since we can limit the simulation time step to about sub pico-second scale. Here, we assume the device-level semiconductor equations can be solved in sub pico-second scale and the signal propagation time scale is on the pico-second scale, which is reasonable. One starting point might be revisiting the 2D case presented by Hasegawa [1971] which has an 63 analytical solution and modifying the substrate current term as biasing- voltage dependent. These models can largely be applied to delay lines, variable phase shifters, voltage-tunable filters, etc. 4) ADI-FDTD method is still evolving. A higher order ADI method is emerging which can perform better in conservative form and provide higher-resolution simulations. 64 Appendix A: Treatment of the Current Term in the ADI-FDTD Scheme In the ADI-FDTD method we presented in Chapter 2, the current-term ( EJ  = ) in the Ampere?s equation E B t E    ?  ? = (A.1) can be discretized in two ways, i.e. for : . / . 0 1 +       =   + + + + + +++ + + + +++ + + (A.3)) 2 ( (A.2) 1 1 ),,2/1(, 1 ),,2/1(, 1 ),,2/1(, )2/1,,2/1(,)2/1,,2/1(, 1 ),2/1,2/1(, 1 ),2/1,2/1(,),,2/1(, 1 ),,2/1(, n kjix n kjix n kjix n kjiy n kjiy n kjiz n kjiz n kjix n kjix EE E z BB y BB t EE   ? ?  Scheme A.2 treats the current term in a fully implicit way and scheme A.3 treats the current term with Crank-Nicholson method by splitting it into half known and unknown. To assess the performance of these two schemes, we simulate the scattering of electromagnetic wave from the metal surface as illustrated in Figure 2.3 of Chapter 2. The excitation is a point Gaussian current source. The simulation set-up is exactly 65 the same as that in section 2.2.2. The electric field was monitored inside the metal surface every 0.1um, and has been Fourier transformed. In Figure A.1, we plot the (a) Crank-Nicholson Method for Current Term (b) Fully Implicit Method for Current Term Figure A.1: Relative electric filed magnitude vs. positions into the metal surface for 10, 20, and 40 GHz obtained from (a) simulation with Crank-Nicholson method for the current term; (b) simulation with fully implicit method for the current term; The solid lines are the analytical solutions and the colored plus signs are the simulation results. 66 electric field magnitude vs. positions into the metal surface for 10, 20, and 40 GHz obtained from simulations with scheme A.2 and A.3, respectively. While the results of (a) (b) Figure A.2: (a) Transient electric field signal inside vacuum for simulation with (top) Crank-Nicholson method for the current term; (bottom) fully implicit method for the current term; (b) Transient current Jz inside the metal for simulation with (top) Crank-Nicholson method for the current term; (bottom) fully implicit method for the current term. 67 the fully implicit treatment of the current term match the analytical solution ( )/exp( yEz  ; )/(2 ? = ). well, we can clearly see the Crank-Nicholson method shows larger skin depth at given frequencies systematically. This comes as a surprise, since at a first glance, these two methods should not differ that much. Figure A.2a provides a closer look at the transient electric field signal inside vacuum monitored during the simulations with the two schemes. They are exactly the same for the two schemes treating the current term with the Crank- Nicholson and Fully-Implicit method, respectively. But in the top panel of Figure A.2b, the interval between 150-300 ps suggests that, inside the metal, Crank- Nicholson method produces much smaller current (4-5 order of magnitude smaller) at the second step. This can be seen more clearly in Figure A.3 where an interval of transient metal current as shown in the top panel of Figure A.2b is amplified. On the other hand, the results from fully implicit method looks fine. Figure A.3: A amplified view of an interval of transient metal current as shown in top panel of Figure A.2b. 68 In the simulation, the electromagnetic wave is excited by current Jz remotely. Only E z , H x , H y components are excited. In the first step of ADI, we sweep through all the y directions to calculate E z as indicated in Equations (A.4). )( 1 )( 1 )(),( ) ( 1 )( ),()( 2 ),( )( 1 2 1 2 2 1 1 (A.4) )2/1,,2/1(,)2/1,,2/1(, )2/1,2/1,(,)2/1,2/1,(, ),,2/1(,)1,,2/1(, ),,2/1(,)1,,2/1(, 2 )2/1,,(,2 )2/1,,(, )2/1,,(,1 2 2 2 2 2 2 2 1 2 2 2,1 1 )2/1,,1(, 1 )2/1,,(,2,1 1 )2/1,,1(, n kjiy n kjiy n kjiy n kjix n n y n x n kjix n kjix n kjix n kjix n x n y n x n x n kjiz n kjiz n y n x n x n kjiz n kjix n kjiz n kjiz BB x t BB y t B t BBg z EE z EE x t Ef BBgEfEd E t BBg EfEd x t c t x t b t x t b x t a dEc EbEa +++ +++ + +++ + + + + + + + + ++    +    =?  =   +      = ++=   ++=   =  +   +=  +   +=   = =? +?+? ? ? ? ?   ?  ?  ? ?  Here, b 1 and b 2 , d 1 and d 2 refer to Crank-Nicholson method and fully-implicit method, respectively. 69 Since only E z , H x , H y are involved in our simulation, we can say 0)( = n x Ef . Also, inside the metal layer, we have the following relation valid.     t b t bbca  2  2<< 212,1 ; 2 ;, . (A.5) In the Crank-Nicholson method, for the first step, when we calculate B x n+1 and B y n+1 ; only B x n+1 uses the newly updated E z n+1 ; B y n+1 is updated with E z n . Suppose initially Ez n is small, when we calculate the ( ?B) n+1 or ), ( 11 ++ n y n x BBg inside the metal, it is of the order n kjiz E t )2/1,,(, 2 +  2   . The main contribution comes from B x n+1 . On the other hand, in the Fully Implicit method: When we calculate the ( ?B) n+1 inside the metal, it is of the order n kjiz E t )2/1,,(, +  2   . Then in the second step, we sweep through all the x direction to calculate E z . In the second step of Crank-Nicholson method, ) 2 (,; 2 ),( 2 )(),( ) 2 ( 1 )2/1,,(, 11 1 )2/1,,(, 111 1 )2/1,,(, 2 )2/1,,1(, 2 )2/1,,(, 2 )2/1,,1(,         t bcaE t BBg E t EfBBgE EcE t bEa n kjiz n y n x n kjiz n x n y n x n kjiz n kjix n kjiz n kjiz  +<<  2  ++= ?+?  ++? + + ++ + + +++ + + + ++ + + + + (A.6) The solution will be 1 )2/1,,(, 2 )2/1,,(, + + + + << n kjiz n kjiz EE . On the other hand, in the second step of the fully-implicit method, 70 1 )2/1,,(, 11 111 1 )2/1,,(, 2 )2/1,,1(, 2 )2/1,,(, 2 )2/1,,1(, ),( )(),( )( + + ++ +++ + + + ++ + + + +  2 ++= ?+?  ++? n kjiz n y n x n x n y n x n kjiz n kjix n kjiz n kjiz E t BBg EfBBgE EcE t bEa     .(A.7) The solution is 1 )2/1,,(, 2 )2/1,,(, + + + + 2 n kjiz n kjiz EE . We can see the solution of using Crank-Nicholson method to treat the current term has very small electric field value in the metal at the second step. If we substitute this value back to Equation A.4, we will get an electric field solution of normal magnitude. Since 2 )2/1,,(, 1 )2/1,,(, 22 22 ),( + + + + ++  >>  2 n kjiz n kjiz n y n x E t E t BBg     , the solution of (A.4) for Crank-Nicholson method at step n+3 will have 1 )2/1,,(, 3 )2/1,,(, + + + + 2 n kjiz n kjiz EE . This explains the observed oscillatory electric field solution inside the metal layer from using the Crank-Nicholson method to treat the current term. 71 Appendix B: Guided Wave Propagation in Two Layer Structure In this section, we present detailed derivations of the dispersion relation (Equations 3.3-3.5) in a 2D two-layered structure as shown in Figure B.1. We denote subscript ?1? as within SiO 2 layer and ?2? as within substrate layer. Define y = 0 as being at the interface between the two layers and positive y is upwards. The TM propagation mode along Z has only Ey, Hx, Ez components. For the Ampere equation, , 1 E t B EB t E     ? = ? =  ?  , (B.1) we can rewrite it in frequency domain as Z Y Metal Layer Silicon Substrate SiO 2 Ground Plane b1 b2 Figure B.1: Side view of the MISS structure. Z is the direction of propagation. To derive analytical solutions, we assume that the structure extends to infinite in the X direction. b1 and b2 are the SiO 2 and silicon substrate thickness. 72 BkEj ii  ?= ?  1 ' (B.2) EkBj i  ?=  , (B.3) where )/( 0 '  j iii += , i = 1, 2. Substituting (B.2) to (B.3), for the two layers, we can get the first two equations for dispersion relation: ' 11 2 0 22 1 ? k=+ (B.3) ' 22 2 0 22 2 ? k=+ (B.4) Also, Ey, Hx, Ez components can be written as 2,1),exp(,, =? izytjEHE izixiyi  . (B.5) In particular, 2,1),exp())exp()exp(( =+= iztjyByAE iiiizi  (B.6) The boundary conditions on the top and bottom part require 0; 21 | === bybyz E . This suggests that Equation (B.6) can be rewritten as z ii ii zzi e b by EE     = ]sinh[ )](sinh[ 0   , (B.7) where 0z E is the longitudinal electric field at Z =0. Here, we removed the )exp( tj  for simplicity. As noted in Pozar [1998], once zi E is given, the other Hx, Ey components can be solved with . y Ej H z i i x =   0 ' (B.7) 73 y E E z i y =   (B.8) As a result, z ii ii z i i xi e b by E j H       = ]sinh[ )](cosh[ 0 0 '   (B.9) z ii ii zzi e b by EE     = ]sinh[ )](sinh[ 0   . (B.10) Here, i = 1, 2 and ?-? sign for i = 1, ?+? sign for i = 2; i =1 represents the field distribution within SiO 2 layer and i =2 represents the field distribution within substrate layer. Up to now, we have two equations for dispersion relation (Equations B.3 and B.4), we need one more equation to solve for  , 1  , 2  . We can use the boundary condition for Ey at the interface between the two layers, namely + == = 0 ' 2 0 ' 1 || y y y y EE  , (B.11) where )/( 0 '  j iii += . The resultant additional equation is 0)tanh()tanh( 22 2 2 11 1 1 '' =+ bb       . (B.12) 74 References: Alsunaidi, M. A., S. M. S. Imtiaz, and S. M. ElGazaly, Electrmagnetic wave effects on microwave transistors using a full-wave time domain model, IEEE Trans. Microwave Theory Tech., vol. 44, pp.799, 1996. Berenger, J. P., A Perfectly matched layer for the absorption of electromagnetic waves, J. Computational Physics, Vol. 114, pp. 185-200, 1994. Elliott, R. S., Antenna theory and design, IEEE Press Series on Electromagnetic Wave theory, John Wiley & Sons, New Jersy, 2003. Groteluschen, E., L. S. Dutta, and S. Zaage, Quasi-analytical analysis of the broadband properties of multiconductor transmission lines on semiconducting substrates, IEEE Proc. International Conference on WSI, pp. 387, 1994. Guckel, H., P. A. Brennan, I. Palocz, A parallel-plate waveguide approach to micro- miniaturized, planar transmission lines for integrated circuits, IEEE trans. On Microwave Theory and Techniques, vol. 15, pp. 468, 1967. Hasegawa, H., M. Furukawa, and H. Yanai, Properties of microstrip line on Si-SiO2 system, IEEE Trans. Microwave Theory Tech., vol. 19, pp. 869-881, 1971. Hasegawa, H., and A. Seki, Analysis of interconnection delay on very high speed LSI/VLSI chips using an MIS microstrip line model, IEEE Trans. on Microwave Theory and Techniques, vol. 32, pp. 1721, 1984. Kwon, Y. R., V. M. Hietala, and K. S. Champlin, Quasi-TEM analysis of ?Slow- Wave? mode propagation on coplanar microstructure MIS transmission lines. 75 Liou, Jenn-Chorng, and Kei May Lau, Analysis of slow-wave transmission lines on multi-layered semiconductor structures including conductor loss, IEEE Trans. on Microwave Theory and Techniques, vol. 41, 824, 1993. Mur, G., Absorbing boundary conditions for the finite-difference approximation of the time domain electromagnetic field equations, IEEE Trans. Electromagnetic Compatibility, Vol. 23, pp. 377, 1981. Namiki T., IEEE Trans. Microwave Theory Tech., vol. 47, pp. 2003-2007, 1999. Pozar, D. M., Microwave Engineering, John Wiley & Sons, New York, 1998. Ramahi, O. M., The complementary operators method in FDTD simulations, IEEE Antennas and Propagation Magzine, Vol. 39, pp. 33-45, 1997. Ramahi, O. M., The concurrent complementary operators method for FDTD mesh truncation, IEEE Trans. Antennas and Propagation, Vol. 46, pp. 1475-1482, 1998. Rappaport T., Wireless Communications: Principles and Practice, Prentice Hall, NJ, 2002. Shao, X., N. Goldsman, O. Ramahi, P. N. Guzdar, A new method for simulation of on-chip interconnects and substrate currents with 3D Alternating-Direction- Implicit (ADI) Maxwell Solver, IEEE International Conference on Simulation of Semiconductor Processes and Devices, pp. 315, 2003. Shao, X., Neil Goldsman, and Omar M. Ramahi, The Alternating-Direction Implicit Finite-Difference Time-Domain (ADI-FDTD) Method and its Application to Simulation of Scattering from Highly Conductive Material, IEEE International 76 Antennas and Propagation Symposium and USNC/CNC/URSI North American Radio Science Meeting: URSI Digest, 358, Columbus, OH, 2003. Shao, X, Omar M. Ramahi, Lin Li, Baharak Mohajeriravani, and Neil Goldsman, Study of Electromagnetic Field Radiation from Apertures using Alternating- Direction Implicit Finite-Difference Time-Domain (ADI-FDTD) Method, IEEE International Antennas and Propagation Symposium and USNC/CNC/URSI North American Radio Science Meeting: URSI Digest, 630, Columbus, OH, 2003. Shibata T., and E. Sano, Characterization of MIS structure coplanar transmission lines for investigation of signal propagation in integrated circuits, IEEE Trans. Microwave Theory Tech., vol. 38, pp. 881-890, 1990. Taflove, A., and S. C. Hagness, Computational Electrodynamics: the finite-difference time-domain method, Artech House, Boston, 2000. Wang, B Z, Y. J. Wang, Wenhua Yu, Raj Mittra, A hybrid ADI-FDTD subgridding scheme for modeling on-chip interconnects, IEEE Transaction on Advanced Packaging, vol. 34, pp. 528, 2001. Wang, G. F., R. W. Dutton, and C. S. Rafferty, Device-Level Simulation of Wave Propagation along Metal-Insulator-Semiconductor Interconnects, IEEE Trans. On Microwave Theory and Technologies, Vol. 50, 1127, 2002. Yee, K. S., Numerical solution of initial boundary value problems involving Maxwell?s equations in isotropic media, IEEE Trans. Antennas and Propagation, Vol. 14, pp. 302, 1966. 77 Zheng F., Z. Chen, and J. Zhang, A finite-difference time domain method without the Courant Stability conditions, IEEE Microwave Guided Wave Lett., Vol.9, pp. 441-443, 1999.