Accurate and Efficient
Modeling of Global
Circuit Behavior
in VLSI Layouts
Zhen-Qiu Ning
TR diss
1715
Accurate and Efficient
Modeling of Global
Circuit Behavior
in VLSI Layouts
Accurate and Efficient
Modeling of Global
Circuit Behavior
in VLSI Layouts
Proefschrift
ter verkrijging van de graad van doctor aan de Technische Universiteit Delft, op gezag van de Rector Magnificus, prof. drs. P.A. Schenck, in het openbaar te verdedigen ten overstaan van een commissie aangewezen door het College van Dekanen op donderdag 20 April 1989 te 16.00 uur door Zhen-Qiu Ning geboren te Guangdong electrotechnisch ingenieur
TR diss
1715
v
1
Stellingen behorende bij het proefschrift:
Accurate and Efficient Modeling
of Global Circuit Behavior in VLSI Layouts
door Zhen-Qiu Ning
1. A poly-silicon gate electrode can be considered as a distributed RC line. As the feature sizes are scaled down, the delay induced by this RC time constant will become more dominant in governing device switching speed and circuit delay. A good approach to modeling this problem is the boundary finite element method. 2. The use of the boundary finite element method is only a beginning in the modeling
of global behavior in VLSI, since the exact behavior of such a device is very important in VLSI design and is difficult to model. The next step may lead to the following modeling combinations: transistors with interconnect parasitics, a boundary with bulk finite element method, Laplace with Poisson equation.
3. It is desired to find a global approach to VLSI device modeling, which is usable for design purposes, rather than a single solution to a specific field problem.
4. An IC substrate is neither a perfect conducting plane nor a dielectric layer. The capacitance models which assume the substrate at ground level will overestimate the capacitance to ground and underestimate the couplings. For more accurate results further research should be directed towards the modeling of the substrate. 5. "Fast running" is one of the most important properties required from CAD tools.
The simplicity of models and the efficiency of an algorithm are both important. The combination of matrix approximation with system model reduction will yield outstanding results.
6. Device modeling is the interface between circuit design and fabrication. Much work has been devoted to modeling the dc characteristics of MOS transistors; however, only a few reports discuss short-channel MOSFET capacitances. This
lack of results is partly due to the fact that the transistor capacitances are related to parasitic capacitances.
7. A fast and cheap way of layout verification consists of a complete extraction of layout data and a switch level simulation with defect testing.
8. CMOS is an elegant and promising IC structure since CMOS ICs usually consume less power than other types of digital ICs and are able to operate over a much wider voltage range. Ion implant is a key technology in CMOS fabrications. It brings advantages and conveniences to the fabrications and leaves cumbersome things to the device modeling.
9. Combination of electron beam technology with CAD tools is an approach to the automatic design and fabrication of VLSI without masks. This approach is very efficient for optimum performance in research and development, in the production of photomasks, or in the direct writing of silicon wafers.
10. Language problems consume much time, but fortunately all the people in the world will speak the same language in the Tomorrow of Mankind.
11. Live to an old age and learn to an old age, yet you cannot learn most of what must be learned.
CONTENTS
1. INTRODUCTION 1 1.1 Modeling of MOS Devices 1
1.2 Parasitic Models 4 1.3 Overview 9 2. BOUNDARY VALUE PROBLEMS IN VLSI DEVICE
MODELING 15 2.1 Field Equations 15 2.2 Integral Equation Formulations for MOSFET 18
2.3 Parasitic Capacitance 21 2.4 Interconnection Resistance 28
3. GREEN'S FUNCTION 33 3.1 Definition 33 3.2 The Bounded Multilevel Dielectric Problem 34
3.3 The Unbounded Multilevel Dielectric Problem 39 4. GALERKIN BOUNDARY FINITE ELEMENTS 51
4.1 Element and Local Shape Function 51 4.2 Solution at the Element Level 54 4.3 Element Assembly and System Solution 57
4.4 Evaluation of Green's Function Integrals 63 4.5 Determination of the Number of Terms Required for Green's
Function 69 4.6 Results and Comparisons 70
5. POINT FITTING BOUNDARY HNITE ELEMENTS 77
5.1 Model Reduction 77 5.2 Solution at the Element Level 79
5.3 Solution of the Connected Problem 81 5.4 Short Circuit Capacitance Matrix 85 5.5 Numerical Considerations 86
6.1 The method of Triangular Decomposition 97
6.2 Approximate Inverse 98 6.3 Results and Comparisons 114 7. ON THE MODELING OF A SHORT-CHANNEL MOSFET BELOW
THRESHOLD 119 7.1 Analytical Solution of Poisson Equation 119
7.2 Boundary Conditions 123
7.3 Discussion 126 8. PARASITIC CAPACITANCES AND THEIR LINEAR
APPROXIMATION 129 8.1 Parallel Conductors 129
8.2 Corners 135 8.3 Crossing Strips 136 8.4 Combination of Corner and Crossing Strips 138
9. DISCUSSION 141 9.1 Complexity and Other Issues 141
9.2 Application to 3D Extraction 142 9.3 More about Interconnection Resistances 143
9.4 Hybrid Finite Elements 143
10. APPENDICES 153 10.1 Appendix 3.1: Solution of Equation (3.8) 153
10.2 Appendix 3.2: Fourier Integral Evaluation 157 10.3 Appendix 4.1: Evaluation of Singular Integrals 162
10.4 Appendix 4.2: Derivation of (4.41) 166 10.5 Appendix 5.1: Evaluation of Green's Function Integral 167
SUMMARY 173 ACKNOWLEDGMENTS 177
BIOGRAPHY 179
-SAMENVATTING 181
-1. INTRODUCTION
1.1 Modeling of MOS Devices
The dependence between geometry and electric field in scaled MOSFET's has generated a great interest in device modeling. A variety of modeling techniques have been presented in the literature [1,2,3,4,5,6]. These concentrate on the behavior of transistors, characterize the device physics and represent the current and capacitance of the short-channel MOSFET.
For short-channel devices, the bulk charge region can no longer be considered as a rectangle lying under the gate, because the depleted charges are shared between the gate and source/drain depletion regions near the source/drain diffusion edges. This results in the well-known short-channel effect. This short-channel effect can be qualitatively explained by a two-dimensional charge sharing concept first proposed by Yau [3]. This concept reduces the two-dimensional problem to one dimension and results in a simple analytical formula to determine the threshold of short-channel MOSFET's [3,7]. Many analytical models are based on this concept [8,9,7,10] . However, this concept applies only to uniform doping cases and is somewhat arbitrary, which limits the accuracy and the universality of these models as the channel becomes shorter.
The MOS devices in VLSI are made in a very nonuniform substrate. The nonuniform doping profile consists generally of (1) a moderately heavy but shallow implant for threshold control (gate implant); (2) a deeper implant for punch-through suppression (field doping); (3) a lightly doped substrate for minimization of capacitance.
For the nonuniform doping, there are two methods which can generally be used in CAD models at present. One is that the implanted profiles are approximated by a rectangular box, the other is that the profiles are described by Gaussian or exponential functions [11,12]. But the models of MOSFET for circuit simulation are invariably based on the box approximation because it is only for this particular case that a simple analytical solution of Poisson's equation can be obtained [1,7,13].
2 INTRODUCTION Ah CO. NB --(a) 0 x, w NAM, NB (b) *i NAW NB «- w -» 1 1 (0 * ƒ
Figure 1.1. Rectangular box approximation of the implanted profiles, xj denotes the
depth of ion implanted region, w the depletion depth, (a) xj <£ W, (b) x\~2>w and (c) xi ■ ■ w.
It has been found that the CAD short-channel MOSFET model with an average constant doping approximation in the channel region results in an acceptable accuracy for circuit simulation if the depth of the ion implanted region is small compared to the total depletion depth in the substrate bias range or if the depletion depth is always confined to a small fraction of the implantation layer. The approximation results in serious error if the depletion region depth is comparable to the implant depth. For example, the body
effect factor y (= y2EsÖ^Beff/Cox) varies with respect to the voltage between source and
body (VBS), since the doping under the channel is not uniform . The effective substrate concentration Nseff may be determined by the implant for low substrate bias VBs and by
the background level NR for larger VBS- However up to the present, most circuit simulators, including SPICE2, do not allow for a nonconstant value of y [14]. This prevents these simulators from being used in the trend towards using higher voltage bootstrap models and from meeting the requirements of lower diffusion line capacitance. To solve this problem, a dynamic average method [15,16] has been presented, in which an effective doping density is determined by a "dynamic average doping transformation" for nonuniform implant devices and the threshold voltage is then calculated by using that effective density according to the charge sharing concept. Since estimating the charges correctly is much more important than perfectly modeling their distribution for circuit simulation, this model can approximate the threshold voltage well. The shortcomings of this method are that the doping transformation procedure must be done for every MOSFET of different channel length fabricated by the same channel implant and that the extracted flatband voltage and effective doping concentration lose their physical meanings. This shortcoming becomes more serious when the extracted parameters are used to calculate the surface potential for applications in the subthreshold current simulation and to determine the capacitances related to the substrate doping concentration [17].
Recently another approach based on an analytical solution of the two-dimensional Poisson equation with approximate boundary conditions has emerged [4,18,19,20,21]. This approach is very promising, since it is free of the drawbacks of the charge sharing. However, the accuracy of the analytical solution is strongly dependent on the simplified assumptions used in the derivation. In existing models there are two sets of boundaries corresponding to two solution domains. The first set contains the boundaries of the silicon layer under the gate given in [19,20] , the second set those of the oxide-silicon composite under the gate [4,21]. For the first set of boundaries, an infinite junction depth and a constant surface potential along the surface channel are assumed in [19]. An improvement is made in [20] by replacing the constant surface potential with the normal electric displacement, but the infinite junction depth is still assumed. These unrealistic assumptions are easily removed by using the boundaries of the second set. However, in the second set the different dielectric permittivities for the oxide and the semiconductor
4 INTRODUCTION
cause difficulties in the mathematical derivation. Thus the Si-Si02 system is always
assumed to be a uniform dielectric [21,4], and then a transformation is made to compensate for the influence. As we know, the transformation using the scaled £s,7eox to
deal with the different dielectric permittivities for oxide and semiconductor in [21] becomes less accurate as the channel is made shorter. An improvement by the image charge method of this transformation is proposed in [4] but, unfortunately, a flat depletion layer is assumed for their derivation. This assumption is valid only for long channel devices. As the channel length is reduced, this flat part vanishes.
In this thesis, we present an analytical Green's function technique to solve the two-dimensional Poisson equation in which the Si—Si02 system is considered as a multi-dielectric composite [22]. The doping can be approximated in terms of a piecewise linear distribution or described as a nonuniform distribution function. The solution that we will show is then more straightforward and more flexible. It can be applied to any MOSFET with either a uniform or a nonuniform doping profile if the doping function is specified.
Based on the analytical solution, we can develop a CAD model which deals with the second-order phenomena accurately and maintains the circuit simulation efficiency.
1.2 Parasitic Models
In the past, VLSI modeling techniques have concentrated on local effects, especially the behavior of transistors [2,21,23,24]. However, with the decrease in feature size and the increase in chip area, global effects are becoming more prominent, with interconnect capacitance playing a major role.
The scaling theory predicts that a reduction of circuit dimensions results in a proportionate reduction in the capacitance that the circuit has to drive. Optimally, the insulator thickness and conductor width and height scale must be reduced at the same rate [25]. Since this has not occurred, capacitances will have higher values than desirable.
The reduced metal thickness causes an unacceptable increase in the resistance of signal lines and power busses. In order to maintain sufficiently low current densities, the height of the conductors is often increased [26]. The aims of reduced line resistance and capacitance are mutually exclusive and will remain so unless there is a significant change
in technology. Thus, the two-dimensional fringing fields and the capacitances between neighboring lines become more important. If these effects are neglected and the parallel-plate formula is used for estimating the capacitances, a very large error may result [27,28]. Some results and comparisons are shown in Tables 1.1, 1.2 and 1.3.
TABLE 1.1. Ground capacitances for single conductor (/ = 4(i, w = 4\i, t = \\i and
h=l[i) (see Fig.5.9)
SPIDER (in fF) 1.695
p-p formula (in fF) 0.552
TABLE 1.2. Ground capacitances for single conductor (/ = 100|x, w = 6(i, t = 1^.,
h = In and d = 3n) (C in fF) (see Fig.4.10)
Dierking and Bastian 30.84 ICPC 30.91 SPIDER 29.38 p-p formula 20.72
TABLE 1.3. Comparison of capacitances for crossing strips (/ = 25n., w = 5|i,
tx =0.5|i, t2 = lM-, hi =0.8(1 and h2=2\i) (see Fig.5.8)
C (in fF) C n Cn C22 Cn Ruehli andB. 8.4 2.2 4.3 2.2 SPIDER 8.307 2.252 4.276 2.252 p-p formula 5.395 1.233 2.158 1.233
Table 1.1 and Table 1.2 are for a single conductor, the first with protecting material over the conductor and the second without. In Table 1.1 the result given by SPIDER [29] is more than 3 times that of the result obtained from the parallel plate formula. In Table 1.2 the result obtained with the parallel plate formula is about two thirds of that from the three other methods [30,31]. Table 1.3 is for crossing strips; the result with the p-p
6 INTRODUCTION
formula is only about half of the result of SPIDER or from Ruehli and Brennan [27].
The increasing chip size and number of circuits per chip cause in addition an increase in the ratio of average interconnection length to minimum feature size. The signal propagation delay introduced by the interconnections may be comparable to or longer than the active device delay. This can be the case for wiring lengths as short as lmm, even with A\im minimum feature size [32]. Consequently, the capacitance of interconnections becomes a dominant effect in governing device switching speed and circuit performance [32,13,33], when the minimum feature size is scaled down to the submicrometer level and the chip size is increased. Two- and three-dimensional models are needed to accurately predict the capacitances.
The modeling of this capacitance encounters the "exterior problem" of Laplace's equation. That is. the electric field is not restricted to a small domain in the immediate vicinity of the conductors but instead extends far beyond that vicinity. The natural mathematical model for this configuration has an infinite domain.
In recent years, methods to deal with this exterior problem have been investigated. One of them is to use Green's function for conductors in an infinite, uniform or layered, dielectric medium [34,35,36,37,38]. This approach maps, using Green's theorem, the exterior problem into a finite boundary problem on the surface of the conductors. Finite elements are used to discretize the resulting integral equation into a finite number of linear equations, which are then solved to obtain the capacitance coefficients. Since the elements are only assigned over the finite boundaries of the exterior problem (surface of the conductors), we call this approach a boundary finite element method to distinguish it from the (traditional) bulk finite element method.
Another way to handle this exterior problem is to use finite differences or finite elements [39,40]. Since these numerical methods can manage only a finite number of difference equations, the medium must be truncated finitely and artificial boundary conditions along the surface of truncation must be imposed. The artificial boundary conditions are iteratively improved on during the search for a solution. Thus this method is more prone to errors induced by the represention of the boundaries. Recently, a new finite-difference procedure presented in [41] has provided a means around this difficulty, in which equivalent terminating capacitors duplicate the behavior of the infinite grids below and above a strip which contains all of the conductors. This leads to an exact determination
of the artificial boundary conditions on the upper and lower surfaces, but the strip still remains an infinite domain because of its unending extension horizontally. Thus two truncating surfaces somewhere to the left and to the right are required, resulting in an approximation. In [41] the potentials along those surfaces are taken to be zero.
Conformal mappings can also be used to obtain analytical expressions [42,43,44], but the method is feasible only for simple geometries [41].
The work in this area can be classified according to the dimensions treated: two or three. Most of the investigations are on the two-dimensional problem [34,35,39,30,41]. For the three-dimensional problems, the best method is to use boundary finite elements [38]. This is due mostly to the fact that the integral-equation mesh need only be placed over the surface of the conductor itself, and not over the entire half-space. In this way, the problem becomes tractable.
Based on different weight functions by which the subspace for projections is spanned, the boundary finite element method corresponds to different weighted residual procedures [45,46]. In these procedures, the commonly adopted methods are the method of subarea [38], the Galerkin method [47] and the point fitting method [37]. Let us partition the surface of the conductors into N elements S j , S2, •••> $N ar|d let W^(x) denote a set of
weight functions. The method of subareas corresponds to the procedure in which the weight function W^x) is equal to one when xsS^ and equal to zero when x outside S*. In the Galerkin method, the weight function W^ is chosen the same as the shape function which denotes the charge density of S^. And for the point fitting method, the weight function W^ = 8(x - x^), i.e. in the surface of the conductors we choose A' points and set functionals equal to the potentials at the points.
For the boundary finite element method, the two-dimensional techniques have proved to be useful for certain simple multiconductor systems in which infinite parallel conductors or planar multiconductors are considered. The three-dimensional modeling techniques are more powerful, but they require larger computer memories and more computer time. In a two-dimensional analysis, the unknowns are surface charges (on the boundary of a cross-section), and are therefore one dimensional. For example, approximately 30 unknowns lead to a good solution for an average problem on a small computer since the matrix requires only 900 words of storage. But the unknowns representing the surface charge will be two-dimensional when the boundary finite elements are applied to a
8 INTRODUCTION
three-dimensional capacitance problem. If the above example is extended to this case, matrices larger than 1000 x 1000 result, requiring excessive storage for a direct matrix inversion. This represents a serious limitation of the three-dimensional techniques [38]. From among the variety of methods shown above, we adopt the boundary finite element method for the three-dimensional problem. It addresses effectively a major problem, the need to reduce the memory requirement and the computing time.
For the three-dimensional problem, the subarea method has been presented in [38], in which the modeling of the charge distribution consists of constant surface charges on a rectangular grid on the conductors' boundaries. The problem with [38] is that a large set of functions is needed to approximate the charge distribution. This is equivalent to requiring more memory and more computing time. For the two-dimensional case, the use of nonconstant charge distributions was considered in [37,48]. It leads to more efficient computations. By replacing the constant surface charge distribution with a one-dimension interpolation, the methods presented in this thesis will be sufficient for the three-dimensional case.
We have written two programs (ICPC and SPIDER), one based on the Galerkin method [31], the other on the point fitting method [49,29]. The Galerkin method results in a positive and symmetric elastance matrix for the finite element equations. It can be inverted by an optimal inversion method with a low order of complexity [50,51]. This leads to accurate and efficient computations.
In the point fitting method, we use a linear distribution of charge concentrated on a web of edges located on the surface of the conductors, and modified into surface charge densities near the nodes. As a consequence, the order of Green's function integrations is reduced, whereas most of the integrals can be evaluated analytically a priori. This reduces the computing time and allows much larger layouts to be dealt with. The program described here permits any geometry of conductors in a stratified medium. It is easy to calculate capacitances of not just parallel conductors, but also of complicated situations like crossing strips, corners, contacts and their combinations. Because there is no context dependent limitation on the usage of assembly elements, this method is ideally suited for extraction purposes (no feature recognition needed). The SPIDER program has been incorporated in an accurate and efficient extractor for submicron VLSI circuits named SPACE.
1.3 Overview
This thesis is organized as follows: Chapter 2 is devoted to the general theory of boundary value problems in VLSI device modeling, including the modeling of transistors, parasitic capacitances and interconnect resistances. Some relations of numerical and analytical importance are established. An appropriate set of Green's function for the Si -Si01 composite is developed in Chapter 3. For the bounded problems, the derivation of Green's function is based on the separation of variable technique and for the unbounded problems, we use the Fourier integral method. The boundary finite element method is presented in Chapters 4 and 5. Chapter 4 corresponds to the Galerkin method and Chapter 5 to the point fitting method. Discretization of the charge distribution and the system solutions are described in detail. Results and comparisons are also given to illustrate the accuracy and efficiency of the methods. The approximate Schur inversion method is studied in Chapter 6. Chapter 7 is devoted to the analytical solution of the Poisson equation in a nonuniform doping short-channel MOSFET. Based on the analytical solution, transistor parameters can be extracted automatically. The parasitic capacitance extract and the linear approximation of the resulting capacitances are presented in Chapter 8. From the accurate models presented in the previous chapters, approximate heuristics for frequently occurring conductor geometries are derived. This can be used to extract models more quickly, at the cost of accuracy. Finally, in Chapter 9 we discuss and pay attention to future developments, especially to the hybrid finite element method which consists of solving a finite element problem in a localized region containing the conductors and the associated low-doping silicon layer in combination with the integral-equation method to obtain the fields outside that region. This hybrid method may lead to a good solution not only for the parasitic capacitance but also for the transistor capacitances which are eclipsed by parasitic capacitances associated with the drain and source junction and by interconnection lines.
References
1. W.L. Engl, H.K. Dirks, and B. Meinerzhagen, "Device Modeling," Proc. of The IEEE 71(1) pp. 10-32 (Jan. 1983).
2. S. Selberherr, A. Schutz, and H. W. Potzl, "MINIMOS-A Two-Dimensional MOS Transistor Analyzer," IEEE Transactions on Electron Devices ED-27(8) pp. 1540-1550 (August 1980 ).
10 INTRODUCTION
3. L.D. Yau, "A Simple Theory to Predict the Threshold Voltage of Short-Channel IGFETs," Solid-State Electronics 17 pp. 1059-1063 (1974).
4. P.S. Lin and C.Y. Wu, " A New Approach to Analytically Solving the Two-Dimensional Poisson's Equation and Its Application in Short-Channel MOSFET Modeling," IEEE
Transaction on Electron Devices ED-34(9) pp. 1947-1956 (Sept. 1987).
5. J.E. Meyer, "MOS models and circuit simulation," RCA review 32 pp. 42-63 (March 1971).
6. D.E. Ward and R.W. Dutton, "A Charge-Oriented Model for MOS Transistor Capacitances," 1EEEJ. Solid-State Circiuts SC-13(5) pp. 703-707 (Oct. 1978).
7. L.M. Dang, "A Simple Current Model for Short-Channel IGFET and Its Application to Circuit Simulation," IEEE]. Solid-State Circuits SC-14(2) pp. 358-367 (April 1979). 8. G. Merckel, "A Simple Model of the Threshold Voltage of Short and Narrow Channel
MOSFET's," Solid-State Electronics 23 pp. 1207-1213 (1980).
9. G.W. Taylor, "The effects of Two-Dimensional Charge Sharing on the Above-Threshold Characteristics of Short-Channel IGFETs," Solid-State Electronics 22 pp. 701-717 (1979). 10. P.K. Chatterjee, P. Yang, and H. Shichijo, "Modelling of small MOS devices and device
limits," IEEPROC. 130, pt. 1(3) pp. 105-125 (June 1983).
11. R.R. Troutman, "Ion-Implanted Threshold Tailoring for Insulated Gate Field-Effect Transistors," IEEE Transactions on Electron Devices ED-24{3) pp. 182-192 (March 1977).
12. G. Doucet and F. Van De Wiele, "Threshold Voltage of Nonuniformly Doped MOS Structures," Solid-State Electronics 16 pp. 417-423 (1973).
13. P. Yang and P.K. Chatterjee, "SPICE Modeling for Small Geometry MOSFET Circuits,"
IEEE Transactions on CAD of integrated circuits and systems CAD-1(4) pp. 169-182
(October 1982).
14. L.A. Glasser and D.W. Dobberpuhl, The Design and Analysis of VLSI Circuits. 1985. 15. P.K. Chatterjee, J.E. Leiss, and G.W. Taylor, "A Dynamic Average Model for the Body
Electron Devices ED-28(5) pp. 606-607 (May 1981).
16. P. Ratnam and C. Andre T. Salama, "A New Approach to the Modeling of Nonuniformly Doped Short-Channel MOSFET's," IEEE Transactions on Electron Devices ED-31(9) pp. 1289-1298 (Sept. 1984).
17. C.Y. Wu, G.S. Huang , H.H. Chen , F.C. Tseng , and C.T. Shin , "An Accurate and Analytic Threshold-Voltage Model for Small-Geometry MOSFETS with Single-Channel Ion Implantation in VLSI," Solid-State Electronics 28(12) pp. 1263-1269 (1985).
18. T. Toyabe and S. Asai, "Analytical Models of Threshold Voltage and Breakdown Voltage of Short-Channel MOSFET's Derived from Two-Dimensional Analysis," IEEE
Transaction on Electron Devices ED-26(4) pp. 453-461 (April. 1979).
19. K.N. Ratnakumar and J.D. Meindl, "Short-Channel MOST Threshold Voltage Model,"
1EEEJ. Solid-State Circuits SC-17(5) pp. 937-948 (Oct. 1982).
20. D.R. Poole and D.L. Kwong, "Two-Dimensional Analytical Modeling of Threshold Voltage of Short-Channel MOSFET's ," IEEE Electron Device Letters EDL-5(ll)pp. 443-446 (Nov. 1984).
21. J.R. Pflester, J.D. Shott, and J.D. Meindl, "Performance Limits of CMOS ULSI , " IEEE
Transaction on Electron Devices ED-32(2)(Feb. 1985).
22. Z.Q. Ning, P. Dewilde, and F.L. Neerhoff, "A New Approach to Analytically Solving the 2D Poisson Equation in MOSFET," Voortgang IOP-IC Modelleringsprojekten , Delft Univ. of Techn. The Netherlands , (Sept. 8,1988).
23. S.E. Laux, "Accuracy of an effective channel length/external resistance extraction algorithm for MOSFET's," IEEE Transactions on Electron Devices ED-31(9)pp. 1245-1251 (Sept. 1984).
24. B.J. Sheu, W.J. Hsu, and P.K. Ko, "An MOS Transistor Charge Model for VLSI Design,"
IEEE Trans. CAD of Integrated Circuits and Systems 7(4) pp. 520-527 (April, 1988).
25. H.B. Bakoglu and J.D. Meindl, "Optimal Interconnection Circuits for VLSI , " IEEE
Trans. Electron Devices ED-32(5)(May 1985).
26. Network Theory Section , "Cooperative Development of An Integrated, Hierarchical and Multiview VLSI-Design System with Distributed Management on Workstations," Internal
12 INTRODUCTION
Report, Delft Univ. of Technology , (1985).
27. A.E. Ruehli and P.A. Brennan, "Accurate Metallization Capacitances for Integrated Circuits and Packages," IEEE J.Solid-State Circuits SC-8 pp. 298-290 (Aug. 1973). 28. Z.Q. Ning, "On the Parasitic Capacitances of VLSI Interconnections ," The ninth CAVE
Workshop, (May 24-27,1987).
29. Z.Q. Ning and P. Dewilde, "SPIDER: Capacitance Modelling for VLSI Interconnections,"
IEEE Trans. CAD of Integrated Circuits and Systems 7(12) pp. 1221-1228 (December,
1988).
30. W.H. Dierking and J.D. Bastian, "VLSI Parasitic Capacitance Determination by Flux , Tubes," IEEE Circuits and Systems Magazine, pp. 11-18 (March 1982).
31. Z.Q. Ning, P.M. Dewilde, and F.L. Neerhoff, "Capacitance Coefficients for VLSI Multilevel Metallization Lines," IEEE Transactions on Electron Devices ED-34(3) pp. 644-649 (March 1987).
32. J. Rubinstein, P. Penfield, JR., and M.A. Horowitz, "Signal Delay in RC Tree Networks,"
IEEE Trans. Computer-Aided Design CAD-2(3) pp. 202-211 (1983).
33. S.P. McCormick, "EXCL: A Circuit Extractor for IC Designs ," 21st Design Automation
Conference Proceedings, (1984).
34. D.W. Kammler, "Calculation of Characteristic Admittances and Coupling Coefficients for Strip Transmission Lines," IEEE Trans. Microwave Theory Tech. MTT-16 (11) pp. 925-927 (Nov. 1968).
35. W.T. Weeks, "Calculation of Coefficients of Capacitance of Multiconductor Transmission Lines in the Presence of a Dielectric Interface," IEEE Trans. Microwave Theory Tech. MTT-18(l)(Jan. 1970).
36. P.D. Patel, "Calculation of Capacitance Coefficients for a System of Irregular Finite Conductors on a Dielectric Sheet," IEEE Trans. Microwave Theory Tech. MTT 19 (11) pp. 862-869 (Nov. 1971).
37. P. Balaban, "Calculation of the Capacitance Coefficients of Planar Conductors on a Dielectric Surface," IEEE Trans. Circuit Theory CT-20 pp. 725-731 (Nov. 1973).
38. A.E. Ruehli and P.A. Brennan, "Efficient Capacitance Calculations for Three-Dimensional Multiconductor Systems," IEEE Trans. Microwave Theory Tech. MTT(Feb. 1973). 39. CD. Tayler, G.N. Elkhouri, and T.E. Wade, "On the Parasitic Capacitances of Multilevel
Parallel Metallization Lines,'' IEEE Trans. Electron Devices ED-32(1 l)(Nov. 1985 ). 40. P.E. Cottrell, E.M. Buturla, and D.R. Thomas, "Multi-Dimensional Simulation of VLSI
Wiring Capacitance," IEDM Tech. Dig., pp. 548-551 (1982).
41. A.H. Zemanian , "A Finite-Difference Procedure for the Exterior Problem Inherent in Capacitance Computations for VLSI Interconnections," IEEE Trans. Electron Devices 35(7) pp. 985-992 (July, 1988).
42. W.H. Chang, "Analytical IC Metal-Line Capacitance Formulas," IEEE Trans. Microwave
Theory Tech. MTT-24 pp. 608-611 (Sept. 1976).
43. M.I. Elmasry, "Capacitance Calculations in MOSFET VLSI," IEEE Electron Device Lett.
EDL-3 pp. 6-7 (1982).
44. T. Sakurai and K. Tamaru, "Simple Formulas for Two- and Three-Dimensional Capacitances," IEEE Trans. Electron Devices ED-30 pp. 183-185 (Feb. 1983).
45. S.G. Mikhlin and K.L. Smolitskiy , Approximate Methods for Solution of Differential and
Integral Equations, American Elsevier Publishing Company Inc., New York (1967).
46. O.C. Zienkiewicz and K. Morgan, Finite Elements and Approximation, John Wiley & Sons (1983).
47. P. Silvester and R.L. Ferrari, Finite Element for Electrical Engineers, Cambridge University Press (1983).
48. P. Benedek, "Capacitances of a Planar Multiconductor Configuration on a Dielectric Substrate by a Mixed Order Finite-Element Method," IEEE Trans. Circuits and Systems CAS-23(5) pp. 279-284 (May 1976).
49. Z.Q. Ning and P. Dewilde, "An Efficient Modelling Technique for Computing the Parasitic Capacitances in VLSI Circuits ," IEEE ISCAS Proceedings 2 pp. 1131-1134 (June 7-9,1988).
14 INTRODUCTION
50. P. Dewilde and Ed.F. Deprettere, "Approximative Inversion of Positive Matrices with Applications to Modelling," NATO ASl Series F34 pp. 211-237 (1987).
51. P. Dewilde , "New Algebraic Methods for Modelling Large Scale Integrated Circuits ,"
2. BOUNDARY VALUE PROBLEMS
IN VLSI DEVICE MODELING
The physical effects of importance in large VLSI circuits are determined by finding approximate solutions to field equations, in particular the Poisson equation and the Laplace equation. In each case, our objective is to determine the potential distribution function <1> in the function of charges in space and at boundaries and then to determine the parameters of transistors and interconnections. When an analytical method is used, O is obtained in closed form as a function which is valid in the solution domain of the problem. When a numerical method is used, O is only determined at all the nodes of the discretized domain of the problem. In both cases, O must satisfy the prescribed boundary conditions.
It is instructive to review the important equations of electrostatics and to illustrate their physical meaning with a simple example. The formulation of the closed analytical solution in the case of the Poisson equation (transistor modeling) is then derived from Green's theorem. In the case of numerical computations, the characteristics of finite elements (including the boundary finite elements and the traditional bulk finite elements) for the Laplace equation are derived from the Green's function method and the method of calculus of variations, respectively. The former corresponds to the unbounded problem (parasitic capacitance problem), the latter to the bounded problem (resistance problem).
2.1 Field Equations
Let V be a volume enclosed by a surface S, D the dielectric flux density (or electric displacement) at a point on the surface, D = e£ with E the electrical field and e the permittivity of the medium, n be the outwardly directed unit normal to the surface at the point and ds be an element of surface area.
For a continuous charge density p(r) in the volume V, Gauss's dielectric flux theorem can be written as an integral form [1,2]
JD .ndS=4njp(.r)dV.
16 BOUNDARY VALUE PROBLEMS
Equation (2.1) can also be written in differential form by using the divergence theorem, which states that for any well-behaved vector field D(r) defined within a volume V surrounded by the closed surface S, the relation
\D .ndS=jV.DdV (2.2) s v
holds between the volume integral of the divergence of D and the surface integral of the outwardly directed normal component of D.
Substituting equation (2.1) into (2.2) yields
J(V.D-p)rfV = 0 (2.3) v
for an arbitrary volume V. It follows that the integrand in (2.3) must be equal to zero and we obtain the differential form of Gauss's law of electrostatics
V . D = p. (2.4)
In addition, we know from Coulomb's law that the electric field can be derived from the scalar potential O by a gradient operation. That is
E = -V <D. (2.5) Equation (2.5) is equivalent to the statement that curl E vanishes [1]. It is known that a
vector field can be specified almost completely if its divergence and curl are given everywhere in space. The behavior of the electrostatic field can thus be described by differential equations (2.4) and (2.5).
Combining equation (2.4) and (2.5) into one partial differential equation for the single function O yields
-V . (e VcD) = p, (2.6a)
VO . Ve + eV24> = p. (2.6b)
If the dielectric is homogeneous (e = cons.), the differential equation (2.6b) reduces to
• V2«D = --£-, (2.7)
e
where V2 denotes the Laplacian, V2 = —r- + —— + —r- for a three-dimensional dx2 dy2 dz7
space. This equation is called the Poisson equation. If there is no charge density in regions of space, the Poisson equation will reduce to the Laplace equation
V2O = 0. (2.8)
In order to illustrate the physical phenomena to be governed by these equations in a device, we consider a purely resistive two-dimensional conductor interconnection problem. A differential element of the conductor is shown in Fig. 2.1 , in which ix and iy
denote the current densities (defined as current/length), I(x,y) denotes the source term (defined as current/area). Since the current flow satisfies the continuity conditions, i.e. a conductor is not able to accumulate electric charges, we obtain
dix dL
- r i + - ^ = / ( x , y ) . (2.9) dx dy
From Ohm's law, the current densities are related to the potential O by the equations
5— = h (2.10a) r ox
and
r dy ■
where r denotes the sheet resistance of the conductor (defined as ohm per square). If the sheet resistance is constant, substitution of (2.10) into (2.9) then yields
18 BOUNDARY VALUE PROBLEMS
32<5 32<D
dx7
dy'= -rt(x,y).
(2.11)This is a two-dimensional Poisson equation. When there is no source (i.e. I(x,y) = 0), equation (2.11) reduces to a two-dimensional Laplace equation. The interconnection resistance problem becomes a boundary problem of the Laplace equation. Given a set of boundary conditions, one may solve the internal potential distribution on the conductors exactly or approximately. Ay \ 'x 1 * 1 1 t l(x,y) i h Ax
Figure 2.1. A differential element of a two-dimensional conductor of interconnections
2.2 Integral Equation Formulations for MOSFET
Consider an 5/—SiO 2 composite under the gate of an n-channel MOS transistor as shown in Fig. 2.2. Domain Q surrounded by boundary C consists of silicon oxide, depletion layer and silicon. tox denotes the thickness of gate oxide, d(x) the depletion depth.
Neglecting an inversion charge and assuming total ionization of impurities in the depletion layer, this domain can be characterized by the Poisson equation with boundary
conditions.
Let the two-dimensional Laplacian V2 = d2 + 32 , A^ denote substrate doping
concentration, £5; the permittivity of silicon , c(y) the distribution of impurity. For the two-dimensional problem, the potential can be described as
V2<&{x,y) = « 0 P(x,y) £5; 0 qNAc(y) for 0 < y < tox forto x<y <d(x) fory > d(x), (2.12a)
with boundary condition
®(x>y) - ^c< for (*>y)on boundary C. (2.12b)
Let Niose denote implanted dose {cm 2), ƒ (y) implant doping distribution function, the
distribution of impurity c(y) can be written as
c(y) = 1
™dose
~N7 f(y)+i
for uniform doping
20 BOUNDARY VALUE PROBLEMS
gate
source drain
Figure 2.2. The cross-section of an n-MOS transistor
The problem represented by (2.12) is an interior Dirichlet problem. The potential <ï> can be solved by means of Green's theorem. For the two-dimensional domain Q, Green's theorem takes the form [1,3]
ƒ ƒ (G V20> - <t>V2G)dxdy = \iG~ - Q>^-)dL,
a dn dn
(2.14)
where G can be any arbitrary scalar field, n denotes the exterior normal on the boundaries.
If the scalar field G is defined as
V2G=--S(x-x',y-y'), (2.15)
and satisfies the boundary condition
then substitutions of (2.12), (2.15) and (2.16) into (2.14) yield the general solution of the two-dimensional Poisson's equation with the Dirichlet boundary conditions for an observation point (x',y')
<S>(x',y') = z(x\y')(\\ &^-G dxdy - j&c^-dL}. (2.17)
a Si c ön
The scalar field G(x',y';x,y) is called the Green's function for this situation. It represents the electrical potential at a point (x',y') due to a unit line charge density at point (x,y). e(jc',y') denotes the dielectric permittivity at observation point (x',y'). If the observation point belongs to the gate oxide, z(x',y') = eox. If the observation point is in
the silicon layer, E(x',y') = e^,.
Equation (2.17) converts the Poisson differential equation into an integral equation. Since in this case the Green's function G can be evaluated (see next chapter), given a set of boundary conditions Oc, the potential O in the Si-SiC>2 composite can be determined
analytically. The transistor parameters, for example the subthreshold current and threshold voltage, can then be derived from the resulting potential distributions.
2.3 Parasitic Capacitance
To determine the parasitic capacitance of interconnections in a VLSI circuit one must determine a relationship between the conductor voltages and the charge on the conductor surfaces . For a three conductor problem as shown in Fig. 2.3(a), this relationship will be expressed as
Qi =C110>1 + C1 2( 0 , - 02) + C i3( 01 -<D3) (2.18a)
02=C2l(«I>2-01)-|-C22<I>2 + C23(<I>2-(I>3) (2.18b)
Ö3=C3l(<I>3-<D1) + C32(<I>3-<I>2) + C33<I>3, (2-18c)
where Qi denotes the total charge on conductor i", Oy the potential of conductor j , and C,j
the capacitance coefficients [4]. If i * j , Cy is a coupling (or mutual) capacitance, if ' = J' Cu 1S a ground capacitance between the conductor i and the reference. The
22 BOUNDARY VALUE PROBLEMS
(b)
Figure 2.3. (a) A three conductor sample, (b) Equivalent circuit.
For an M-conductor system, the relationship is denned in more general form by
M
Qi = Ca*,- + 2Cy(*,- - Oj) (i = 1,2,....,M), (2.19)
where the M x M capacitance matrix is the node (or conductor) admittance capacitance matrix. Due to energy conservation, the capacitance matrix will be symmetric [2] (i.e. Cjj = Cp for all numbers i and j) and for a complete characterization of the given M-conductor system only M (M+l)/2 capacitances have to be considered.
In our case, the capacitance of interconnections is assumed to depend only on the geometry of the conductors and dielectric interfaces. (If silicon substrate is considered as a nonlinear semiconductor, i.e. can not be considered as a perfect conductor, the assumption is not valid. This investigation is reserved for future studies, see Chapter 9.) The equation (2.19) can be rewritten in the form
M
where Csij denotes a so-called short circuit capacitance. The short circuit capacitance of
a conductor is related to the total charge on the conductor when it is maintained at unit potential, all other conductors being held at zero potential.
By comparing (2.19) with (2.20), the admittance capacitances are obtained from the short circuit capacitances as
N
Cii=IlCsij i=l,2,..JV, (2.21)
C ■ - - C ■■ i*j- (2.22)
From (2.21), (2.22) and (2.20), the parasitic capacitance problem to be considered reduces to the determination of the charge Q; on each conductor in the function of the potentials <fy.
Consider the configuration of interest to integrated circuits, which is often accurately described as a stratified medium with a conducting ground plane and a sandwich of dielectric layers in which the conductors "float", as shown in Fig. 2.4.
y» 0 x
™
M
e2 SiliconFigure 2.4. Multiconductors in the Si-Si02 composite
Since the dielectric is homogeneous (e, = cons.) in each layer and without space charge , in finding the relationship between the conductor voltage and charge one encounters the "exterior problem" of the Laplace equation. The electric field is not restricted to a small
24 BOUNDARY VALUE PROBLEMS
domain and instead extends infinitely. For our situation (see Fig. 2.4), the conductors are located in the semi-infinite half-space above the silicon.
Figure 2.5. A conductor is located in the semi-infinite half-space.
Let V, denote the volume of the j'th conductor, V the domain outside V; (see Fig.2.5). The domain V,- is surrounded by the closed surface 5,-. The domain V is bounded by two surfaces, one is the surface of the conductor (S,-); the other at the infinite conducting ground plane and at infinity, say C. The ith conductor has constant potential <I>,. Since the total energy stored in the field is finite, we expect the corresponding potential O to vanish at infinity. The potential in the conducting ground plane is chosen as a reference, let it be zero. Then the "exterior problem" can be written as
V2<D(p) = 0 pinV, (2.23a)
with boundary conditions
C>(p) = 0, p on ground plane or at °° (2.23b)
<ï>(p) = d>;, pon Si. (2.23c)
Because of the relatively low frequency range in which integrated circuits are operated, the electric field inside conductors can be neglected. Hence, the charge is essentially concentrated on the surface of the conductor. Let Oi(q) denote the surface charge density on the conductor surface S;, n the outwardly directed unit normal to the surface S; at point q. Another boundary condition is then derived from Gauss's law [2]
^l{q) = - ^ 3 l lS r (2.24)
e dn
The problem described by equation (2.23) and (2.24) is an exterior Neumann problem in three dimensions. The potential distribution can be solved by many methods as shown in Chapter 1. Here we only discuss the boundary finite element method. The (hybrid) finite element method is left for Chapter 9 .
Let p denote the observation point, q the integration variable, n' the outwardly directed unit normal to the surface C at point q. Since the boundary surface area C is infinite, the Neumann boundary condition becomes [1,3]
^ 1 = 0, qonC. (2.25)
an
If we define an appropriate Green's function similarly to (2.15) and (2.16),
V2G=--8(p-q) (2.26a)
e with
G = 0 at boundary C, (2.26b) Green's theorem applied to domain V yields [3,5] (see Fig.2.5)
26 BOUNDARY VALUE PROBLEMS
e j [ - < ! M ^ + ^ ] 4 ( r t =
3n
* ( p )
i f p e V y<D, ifp e Si. n i f P e vi
(2.27)
Similarly, Green's theorem applied to domain V,- results
eJ[G<p,,)^-*fo)^k£]*fo)
dn *<P> if p e V,-2 0 0>(p) if p s 5,. (2.28) Up e VWhen we take in equation (2.28)
_ Jl if P e V,-^ ' ~ 10 if p e elsewhere, equation (2.28) becomes St dG(p,q) dn ds(q)--1 if p e Vt \ ifpeSi. 0 if P e V (2.29) Hence
, ƒ JSiMl „ , „ , . .
if p e S,-if p ê V;- (2.30)<D(p) ifp e V
<D, if/? e Si. (2.31) 0 if p é V
Thus if we substitute equation (2.24) into equation (2.31), we obtain
®(p)=JG(p,q) Oi(q)dS (q) p s V or p e S,-. (2.32) s,
This is a formal solution of (2.23) in the form of the integral. The potential becomes the Green's function integral for all charges on the surface of the conductor. It transfers the exterior problem into an "interior" problem in a finite domain (surface of the conductor). Since the Green's function is the potential due to a unit charge and there is no space charge in the domain V, the integral represents the total of the contributions from all surface charges.
For an M-conductor system, the potential of the jth conductor is then written as
*;<P>= £ ƒ Gip, q)Oi(q)dSi(q) j=l,2,....,M, (2.33) 1 = 1 S,
where G;(<7) denotes the charge density on the j'th conductor surface 5, .
To determine the relationship between potential and charge for a known potential, we have to invert equation (2.33). This can be done numerically by the boundary finite element method, since the Green's function can be found (see next chapter). In the boundary finite element method, the surface of conductors is first divided into a number of elements and a shape function is assumed for the solution of charge density within and on the boundary of the element. Numerical element equations are then built up from equation (2.33). The elements are assembled and a system of equations is obtained representing the given conductor system. Finally, the equations are solved by a standard method of matrix inversion, for example the method of triangular decomposition, or by an optimal inversion method from which a reduced charge-potential relation is derived (see Chapter 6).
-zJG(p,q)
s,
dO(.q) dn ds(q)-.
28 BOUNDARY VALUE PROBLEMS
2.4 Interconnection Resistance
For VLSI verification, the resistance of metal interconnections can nearly always be neglected, only polysilicon and diffusion resistances may be important. Let Q. be a two-dimensional resistive region surrounded by a boundary C (Fig.2.6). Some parts of the boundary consist of highly conducting contact regions. Our purpose is to determine the resistance network between all the contacts.
boundary C / / / / contact
n
Figure 2.6. A polygon with contacts (terminals).
Based on the principles of current conservation and continuity (see Section 2.1), the calculation of potential distribution in the domain £2 becomes a mixed boundary value problem for the Laplace equation, one part is the Dirichlet boundaries at the contact regions, the other is a vanishing normal derivative of the potential at the remaining boundaries. That is
V2O(p) = 0, pinQ.; (2.34a)
with boundary conditions
<&(p) = <t>i, p on contact i, i = 1,2,...
-r— = 0, p on remaining boundaries on
(2.34b)
The potential governed by equation (2.34) can be found numerically by a finite difference [6] and a finite element [7,8] method. The finite difference method is well suited for orthogonal geometries and regions of uniform electrical properties, but certain difficulties arise when the shapes involved have nonorthogonal sides and the region is made up of irregular geometries with different square resistance and different boundary conditions. The finite element method, unlike the finite difference method, can easily handle complicated geometries with variable material properties and mixed boundary conditions. It also meshes well with extraction. We only consider here the finite element method.
In the finite element method, the domain Q is first divided into a number of elements, usually triangles, and a shape function is assumed for the solution function within and on the boundaries of the elements. Let JV denote the total number of elements, fyt the local
shape functions and 0(p^(-) the coefficients to be determined. The potential is
approximated in terms of triangular shape function as
*(p)4l%)/ü- (2-35)
j=\ i=\Numerical element equations can then be derived by using the well-known principle of minimum potential energy, since the minimization method is both simpler and mathematically more appealing [9] for this problem. Solving the coefficients will yield an approximation to the potential distribution. A resistor network is finally determined by the relationship between currents and the voltages.
It can be proved that the minimum energy principle is mathematically equivalent to the Laplace equation in the sense that a potential distribution which satisfies the latter equation will also minimize the energy, and vice versa.
Within the domain Q, the energy is given by
W(&) = — e ƒ I V<D 12dQ, (2.36)
Consider a perturbation O + x¥, where *F is equal to zero at every boundary point. The
30 BOUNDARY VALUE PROBLEMS
W (&+V) = j E ƒ I V(OH->F) 12dQ 2 Ü
= W(®) + WQ¥) + e ƒ VG>V¥dQ. (2.37)
From Green's first identity , we have
JV<DWdQ = J > - | ^ - < i s - J > V2Od£2 = 0, (2.38)
since V20 = 0 in £2 and *P = 0 on the boundary. Substitution of (2.38) into (2.37) yields
W (0H-T) = W (O) + W OF). (2.39) Since for any meaningful perturbation WQV) > 0, W(O) is indeed the minimum value of
energy. Minimization of the energy associated with the potential <I>(p) given in Equation (2.35) can then determine the coefficients, and implicitly determine an approximation to the potential distribution.
References
1. J.D. Jackson, Classical Electrodynamics, John Wiley & Sons, Inc. (1975).
2. E. Weber, Electromagnetic Fields Theory and Applications, New York: Wiley (1957). 3. I. Stakgold, Boundary Value Problems of Mathematical Physics, New York: Macmillan
(1968).
4. P. Balaban, "Calculation of the Capacitance Coefficients of Planar Conductors on a Dielectric Surface," IEEE Trans. Circuit Theory CT-20 pp. 725-731 (Nov. 1973). 5. A.C. van der Woerd, J.P.M, van Lammeren, R.J.H. Janse, and R.H. van Beynhem,
"Calculation of the Resistance Value of Laser-Trimmable Planar Resistors in an Interactive Mask-Layout Design System," IEEE Journal of Solid-State Circuits SC-19(4)(August, 1984).
6. S.P. McCormick, "EXCL: A Circuit Extractor for IC Designs ," 21st Design Automation Conference Proceedings, (1984).
7. CM. Sakkas, "Potential Distribution and Multi-Terminal DC Resistance Computations for LSI Technology,'' IBM. J. Res. Develop. 23(6) pp. 640-651 (Nov. 1979).
8. P. Kazil and P. Dewilde , "A Simple and Fast Method for Obtaining Resistance of VLSI Interconnect," IEEEICCD Proceedings, pp. 342-345 (Oct. 6-9,1986).
9. P. Dewilde , "New Algebraic Methods for Modelling Large Scale Integrated Circuits ,"
33
3. GREEN'S FUNCTION
For the limited purposes of this chapter we look mainly at the Green's function for a multilevel stratified dielectric medium (an approximation for VLSI). We consider first the two- and three-dimensional problem in a homogeneous medium and give a review of the definition of the Green's function and its physical meaning. The Green's function for the Dirichlet problem inside a rectangular box in the multilevel stratified dielectric medium as shown in Section 2.2 is then directly sought through the separation of variable technique [1]. Finally, the problem of an appropriate Green's function for the unbounded multilevel dielectric regions is treated by the Fourier integral method [2,3] in which the Green's function is expressed in terms of cylindrical coordinates and satisfies the boundary conditions as in the conventional boundary value problem. Using this method, Patel [4] has obtained an appropriate Green's function for the case when the ground plane is at infinity. The case in which the ground plane (silicon) is at the dielectric surface is the subject here. Compared with the images method [2,5], the Fourier integral method is less cumbersome and its results are in a more practical form [2].
3.1 Definition
Let V be an open, bounded region in three-dimensional space and let its boundaries be denoted by S. The Green's function for the Laplacian in V is the solution G (p,q) of the boundary value problem
V2G = - - 8 ( p - o ) , p and q in V; G=0iorponS. (3.1)
e
The boundary value problem (3.1) has a simple physical interpretation in electrostatics. We can view G as the electrostatic potential at p due to a unit charge at q when the boundary potential vanishes, which is the case, for instance, when S is a grounded shell or at infinity. For a two-dimensional problem which arises as a degenerate version of a three-dimensional problem when the applied source and the geometric configuration are independent of one of the Cartesian coordinates, the Green's function G may be thought of as the potential due to the unit line charge density of an infinite line source parallel to an infinite surrounding surface at zero potential.
3.2 The Bounded Multilevel Dielectric Problem
Consider first a two-dimensional configuration as shown in Fig. 3.1, which is the SiC>2—Si composite under the gate of a MOS transistor. The domain £2j consists of silicon oxide, the domain J2j °f silicon, Q is the union of Qi and Q.2, and C denotes the
boundary of £2. L and b are the size of the domain Q. and the interface between domain £2i and Q2 is at y = tox. The Green's function G(p,q) is simply the potential at point
p(x,y) due to the unit line charge density at point q(x',y') when the potential on boundary C vanishes. Since the potential and the normal component of the dielectric flux density are required to be continuous on the interface, the Green's function for domain Q can be written as V2 G=-—5(x El x',y-y') E2 8(x-x',y-y') if q e Q\ ifqe Q2 (3.2a)
with boundary conditions:
G{p,q) = 0 on boundary C, G1(p,q) = G2(p,q) at y = tox, (3.2b) (3.2c) and dGi dG2 e1-r— = e2-^— aty = tox, ay ay
where the Laplacian V2 = dx + 32 for the two-dimensional problems.
3.2. The Bounded Multilevel Dielectric Problem 35 (0,0) (L,0) (0, tox) (0,b) y* E2 ?(*',ƒ) J><x,y) Q2(5/) boundary C
Figure 3.1. A rectangular box under the gate for an n-MOSFET
For the finite interval in the x direction with G (x,y; x',y') = 0 at x = 0, x = L for all y, we try a solution of the Green's function as a Fourier series of the form:
G(x,y;x',y')= £ Am(y, x',y')sin(amx) (3.3a)
with
mn
o.m = —, m = l,2,. (3.3b)
Substitution of equation (3.3) into (3.2a) yields
- d2Am m=\ dyz
■a2mAm) sin (amx) = 6(x-x'; y - / ) .
£1,2
(3.4)
functions \ \L/2 if m = n ƒ sin (amx) sin (anx) dx = -L if m * „ (3 5 a) x=0 I and L
| -5(JC -A:'; y -y') sin (anx) dx = -sin (anx') 8(y - ƒ ) , (3.5b) x=0
we obtain an ordinary differential equation for An if we multiply both sides of (3.4) by
sin (anx) and integrate over L :
(3.6) where An Let d2An dy2 = An(y,y').
*-ï
-a2nAn = - 2 Z\,1L Fnsin(anx'),equation (3.6) may be rewritten
d2Fn dy2 a2Fn = -as 1 6 0 El.2 (3.7) 5 ( j - / ) . (3.8)
Fn(y,.y') is determined by solving the ordinary differential equation (3.8) (see Appendix
3.1). Substitution of Fn into equation (3.7) gives An. Let n equal to all m and substitute
(3.7) into (3.3a), we finally arrive at
G(x,y -y,y') = j - £ sin (amx) sin (anx') Fm(y,y'), (3.9) L m=l
3.2. The Bounded Multilevel Dielectric Problem 37
where am = mn/L, m is integer and
F
-
( y y ) =W . y ) if^Pz.
(3-
10)F)*}' and F<? are
F^=B^sinh(amy) ifO<y<y' (3.11a)
F W = f i £ W ( am( y - ro ; c) )
+ //£W/i(am()>-fo;c)) ify'<y<tox (3.11b)
F W = B L3 )« ^ ( am( ö - j ) ) iff0x<y<ö (3-llc)
and
where
F$=AWsinh(.amy), ifO<y<tox, (3.12a)
F^=A^sinh{am{y-t0X))
+ E$cosh(am(y-tox)), titox<y<y', (3.12b)
Fk2)=/lL3)«>>Mam(6-)0) ify'<y<b, (3.12c)
m Ei sinh0iCosh(-02) + E2Cosh8i sinh(-82)
*£> = — . ! —tL, (3.13a)
Htf
Eisinh9isinh(am/) (3.13c)and
B 0)^ £isinh(gm/) Ai
Ai = amE] fE] sinh9] cosh (amtox) + E2cosh9! sinh (amtox)},
,., E2sinh93 i ( i ) — _ t^m — E$-A2 ' E] sinh93Cosh(amfox) &2 ' E2sinh93sinh(amfOT)
(3. Ei cosh (amfox)sinh92 + E2Sinh (amrox)cosh92
= _
4(J)
-&2 -otmE2/'E2«n^(cxmfox)cosh9i + EiCo^/!(amfox)sinh9iJ
am(b-tox) Omiy'-tax) am(b-y') (3.13d) (3.13e) (3.14a) (3.14b) (3.14c) (3.14d) (3.14e) (3.15)
If E! = E2 = E , then the domain Q is uniform dielectric medium. Let tox = 0 for
q(x',y')e&2 (or tox = biorq(x',y')&£l\). We obtain from (3.14)
AO) - 4 ( 2 )
(2) _ sinh(gm(fr-y') Eccmsinh(amfc) '
(3.16a)
3.2. The Bounded Multilevel Dielectric Problem 39
4
2 )=0 (3.16c)
and
,,, sinh(am)>')
eamsinh(amfe)
Substitution of (3.16) into (3.12) results in
,2) s i n h ( c w ) s i n h ( am( 6 - / ) ,
^m = Tv;—rr , y < Y (3.17a) e amsinh(ccmft)
and
,~ sinh(am/) sinh(am(è -y)
E amsinh(amo)
These two formulas are the same as given for the uniform dielectrics [6]. The Green's function for the SiO 2-Si composite as shown in (3.9) is then simply reduced to the same form as for the uniform dielectrics.
3.3 The Unbounded Multilevel Dielectric Problem
It follows from the linearity of equation (3.1) that in the open domain V, G is the sum of the potential due to the unit source at q in free space and the potential due to the charge induced on the boundary S [7]. It can be written as
G(p,q) = E(p,q) + u(p,q), (3.18) where E(p, q) denotes the free-space solution and u(p, q) is a harmonic function in V.
They satisfy the following equations, respectively,
V2£(p, q) =•-£-! 5(p - q), for all p and q (3.19)
V2u(p,q) = 0, for/? in V. (3.20)
For the situation given in (3.1), the boundary condition is G =E+u=0 on S. The problem of finding Gip,q) is thus reduced to that of finding the harmonic function u(p, q) in the domain V, when certain special boundary values on 5 are assigned.
Consider a multilevel stratified dielectric medium as shown in Fig. 3.2, in which the domain V[ may be silicon oxide, Vi silicon nitride and V3 air or protecting material. Each domain is supposed as a uniform dielectric (e = constant).
V3, (£3)
Let
Figure 3.2. A multilevel stratified medium with a unit source.
P = [ U - * T + ( y - y T ]A.
Since the physical situation is symmetrical about the z axis (see Fig. 3.2), the Laplace equation in cylindrical coordinates can be written as
„2 32« 1 du d2u n
V2u = —T + — — + —T = 0. dp2 P dp dz2
3.3. The Unbounded Multilevel Dielectric Problem 41
The Laplace equation is known to have a unique solution. Hence the cylindrical-symmetric solution, if it exists, will be the only one. We try a solution of (3.21) of the type
u=h(z)f(p). (3-22) Substitution of (3.22) into (3.21) yields
, . , d2f 1 , , . df , , , , d2h n h (z) —-f + — h (z) -f- + f (p) —— = 0. dp2 p dp dz2 That is 1 ,r d2f 1 df 1 d2h f(p) dp2 P dp h{z) dz2
Since the left side of the equation is a function of p and the right side is a function of z, they must be equal to a constant. Let it be equal to -m2 for convenience . Hence
variables separate into two ordinary equations
d2h 2U n — T - - m h = 0 , dz2 (3.23) and J ^ + ±-C + m2/ = 0. (3.24) dp2 P dp
Equation (3.23) leads to two independent basic solutions. One increases in magnitude as z increases; the other decreases in magnitude . Let
Z\ =Z -z',
h(z) = e"u' (3.25a)
and
h(z) = e~mZl. (3.25b)
The radial equation (3.24) can be put in a standard form by the change of variable Pi = mp. Then it becomes
dpi Pi dpi
This is a Bessel equation with zero order, and its solutions are called Bessel functions of zero order. It also leads to two independent basic solutions. One is a Bessel function of the first kind and zero order; the other a Bessel function of the second kind and zero order (since the order is an integer). For any domain as shown in Fig. 3.2, the solution of ƒ (p) should be satisfied for the range of p from zero to infinite. Since the Bessel function of the second kind and zero order becomes infinite at p - 0, only the Bessel function of the first kind and zero order is acceptable from the physical point of view. So
f(p) = Jo(mp), (3.26) where ƒ o denotes the Bessel function of first kind and zero order.
Hence for each m, substitution of (3.25) and (3.26) into (3.22) yields two independent basic solutions for equation (2.21)
u=emZlJ0(mp) (3.27a)
and
u = e " " ' /0( « P ) - (3-27b)
Any linear combination of these basic solutions will provide a general solution to equation (2.21). If we introduce arbitrary function H\(m) and H2(in), a general solution
3.3. The Unbounded Multilevel Dielectric Problem 43
can be written as
u(p,g)= ƒ [Hx(m) emz' m=0
+ H2(m)e~mZl]J0(mp)dm. (3.28)
The unknown functions #i(/n) and H2(m) can be determined by the boundary conditions (see below).
The free-space solution E in equation (3.18) satisfies the Laplace equation with boundary condition E = 0 at infinity. That is
1 (3.29)
47te-y z\ + p2
It can also be expanded as a Fourier integral in terms of the linear superposition of the basic solution (3.27). Since E = 0 at z = °°, only the basic solution which decreases in magnitude for I z 11 —>°° is valid. So we have
E(p,q)= ƒ H{m) e m'Z l' J0(mp)dm. (3.30) m=0
For the free-space field, one may take [8,9]
H(m)=-^-. 4rtE Hence
E(p,q) = Y~ ƒ e mhi* J0(mp)dm. (3.31)
4KE m=0
Finding the Green's function for the composite as shown in Fig. 3.2 now becomes relatively easy. If q (x',y',z') in the domain Vt, it is clear that the Green's function in the
