Determination of Potential Profile In Planar Electronic Structures Using A Semi-Analytical Technique

In this paper, a semi-analytical technique known as the Method of Lines (MoL) with uniform and non-uniform discretization schemes is developed. The aim is to determine static potential profile in planar electronic structures. Even though this method has been known for some time, there has been reported work on its application to planar semiconductor device analysis for voltage profile determination. Since most current electronic devices are manufactured using planar and quasi-planar technology, the proposed algorithm is well suited for device analysis prior to fabrication. Compared with known popular methods such as Finite Difference and Finite Elements methods, the proposed technique is relatively simple, more accurate and unlike other methods, has no convergence problem. In addition to this, its semi-analytical nature, which consists of reducing one computing dimension, allows saving significant memory and computation time. Typical planar electronic structures are considered to demonstrate their suitability for these devices, and the obtained results are presented and discussed.


Introduction
The most popular techniques used for numerical simulations of various electronic devices are the Finite Difference and the Finite Element Methods (FDM and FEM).Although, some other techniques, such as the Boundary Integral and Monte Carlo methods have been used in some cases, the two former methods are the most widely used despite of some of their disadvantages.
The main drawbacks of these methods are long running time, large memory space requirements and conver-____________________________________ Corresponding author's e-mail: hadj@squ.edu.omgence problems.In the present work, an alternative seminumerical approach is used for the simulation of electronic planar devices.The approach is based on a mathematical method known as the Method of Lines (MoL) (Ames, 1977, Sadiku, et al. 2000;Vietzorreck, et al. 2000;Barcz, et al. 2003;Vietzorreck, et al. 2004;Pascher, et al. 2005;Yan, et al. 2005;Gonzalo, et al. 2006;Plaza, et al. 2006 andChen, et al. 2007)  R. Pregla and co-workers for the analysis of planar microwave structures (Schultz, 1980, Pregla and Pascher, 1989, Pregla, 2002).For a given partial differential equation subject to some boundary conditions, all but one of the independent variables are discretized.This allows finding a one-dimensional analytical solution.The semianalytical nature of the MoL makes the computational efforts much less intensive than the above-cited methods applied so far to microwave structures.This method has been used recently for the necessity of reducing long calculation times and therefore to significantly reduce the computation time (Zhou, et al. 2000;Pascher, et al. 2005;Yan, et al. 2005;Gonzalo, et al. 2006;Plaza, et al. 2006 andChen, et al. 2007).The method of lines has been used by several authors for the analysis of microwave structures in both quasi-static and full-wave modes (Schultz, 1980;Pregla and Pascher, 1989;Pregla, 2002;Barcz, et al. 2003;Vietzorreck andPascher, 2004 andYan, et al. 2005).However, the method has not been extended to the voltage profile determination of planar semiconductor devices.In the following sections, the application of the MoL to the analysis of planar structures will be demonstrated based on a two-dimensional Poisson's equation.Since the discretization schemes are important parameters in the efficiency of the algorithm, the uniform and nonuniform schemes will both be presented.Finally, some of the obtained results for basic semiconductor planar structures will be depicted.

Technique of Analysis
Let us introduce the approach of the Method of Lines based on solving the very popular Poisson's equation given by Eq. (1) for the domain shown in Fig. 1.
(1) If we uniformly discretize the variable x, then the functions (x ,y) and f(x ,y) in Eq. ( 1) are transformed into the sets (x i ,y) and f(x i ,y) along the lines x=i.hwhere i=1,2,3…N.
N is the total number of lines within the structure and h is the discretizing size interval as shown in Fig. 1.
This operation transforms Eq. (1) into a system of N differential equations of the form: (2) If we let (x i ,y) =  i (y) and f(x i ,y) = f i (y) then, the i th difference approximation on the discretized variable may be written as: (3.a) and its second derivative may be written as:  2) can be written in the following matrix form: (5) where P is a second order differential operator, which is an NxN tri-diagonal matrix of the form: x where p 1 and p 2 are determined by the left and right boundaries of the structure of Fig. 1.
Since P is a real-symmetric and tri-diagonal matrix, The following two possible conditions on  0 and  N+1 may be met: there exists a nonsingular NxN matrix T such that the matrix  given by: ( 7) is diagonal and where the elements  i are the eigenvalues of P.
T is a matrix of the corresponding eigenvectors and T t is its transpose.
For a non-uniform discretization scheme, the approach presented above has to be slightly modified.
Hence, the second derivative expressed by Eq. (3.b) for the uniform discretization will now be expressed as given below for the non-uniform case: (8) where e i = (h 1 + h i-t ) / 2 denotes the i th interval size between the dotted lines shown in Fig. 2 and h where h is an assumed normalized discretization interval.
Considering the Dirichlet-Dirichlet lateral boundary conditions of the analyzed structure, which means      , Eq. ( 9) takes the following matrix form  18) x y

Figure 2. Non-equidistant discretization pattern with the method of lines
The difference Eq. (3.a) for the non-uniform case is firstly normalized by multiplying both of its sides by h h h i and to get the relation: Similarly, the second order operator given by Eq. ( 8) is normalized by multiplying both of its sides by the factor i h e / h to get the following expression: and D t is the transpose of D.
Inserting the expression of  x given by Eq. ( 10) into Eq.( 16) gives: then Eq. ( 19) may be written as: Inserting the expression of  xx from ( 19) into the system given by Eq. ( 2) gives the following matrix form (23) where f is a vector of functions obtained from discretizing f (x,y).Since D xx is a tri-diagonal matrix, then the system given by Eq. ( 24) needs to be decoupled before being solved.Furthermore, since D xx is real-symmetric, there exists an orthogonal matrix T such that the elements of the diagonal matrix  are the eigenvalues of P.
T in this case, is also the matrix of the corresponding eigenvectors and T t its transpose.Using the above result, the system given by Eq. ( 24) is decoupled by pre-multiplying it by T t and obtaining the following system: The general solution of each equation of the system (31) takes the form: where the coefficients a i are determined by substituting V i of Eq. (31) by the above expression.
In the end, the vector of potentials  is obtained from V using the following expression: (34) which expresses analytically the solutions along the discretization lines as a function of y.
Based on the above algorithm, a computer program has been developed for the simulation of semiconductor devices to determine the potential distribution within planar structures.It includes the following steps: 1. Enter the structure geometry and the dimension of each layer, the type of boundary conditions, the electric charge profile as well as the discretization interval h. 2. Perform discretization by finding the total number of lines and construct the second order operator P. 3. Determine the eigenvalues of the matrix P, the elements of  and the corresponding matrix of eigenvectors T. 4. Find in the transformed domain the space charge density and homogeneous horizontal surface boundary vectors by pre-multiplying by T t . 5. Compute the constants A i , B i by considering the field and potential interface continuity conditions as well as the surface boundary conditions.6. Perform the inverse transformation to find the analytic solution in the original domain along each discretization line.
Multiplying both sides of Eq. ( 23) by r e 1 yields:   2 then Eq. ( 28) is further transformed into a set of N ordinary differential equations as where A i and B i are constants depending on the horizontal side boundaries.The specific solution V P i is expressed as a linear combination of the functions F i and its linearly independent derivatives

Case Studies
In this section, the developed algorithm is applied to a set of typical planar structures used in current electronics technology.The objective of the simulation is to determine the potential profile throughout the structures having different physical, geometrical and boundary conditions.Using a uniform discretization scheme with an interval size h = 1 m results in a total of 30 intervals, the following lateral boundary conditions have been considered: Neumann-Neumann The results showing the obtained potential profile contours are depicted in Fig. 4.
With reference to this figure, one can observe that the change in the potential profile is strongly related to the boundary conditions.This is in accordance with most practical situations where a boundary of a semiconductor The effect of boundary conditions on the potential profile has been first considered by finding the potential profile throughout a two-layer planar structure as shown in Fig. 3.
After this, two-layer symmetrical structures shown in Fig. 5 have been simulated.Note that in both structures the layers near to the groundsides have a constant charge density , whereas the charge profile throughout the second layers varies linearly with y as    A(y-H 1 ) + B. The dimensions and the physical parameters for both structures are:  1 = 3.5 x 10 11 C / cm 3 ; A = -10 13 C / cm 4 B = -10 11 C / cm 3 ; = 11.0   ; L = 30 m  mH 2 = 15 m; and V a = 20 V.For the structure 5-a, 1 1 = 1 2 = 10 m and the boundary conditions are Dirichlet-Dirichlet; whereas for the structure 5-b, 1 1 = 1 2 = 6 m and the boundary conditions are Neumann-Neumann.
The analysis is carried out with 63 uniform intervals for the device 5-a and 65 intervals for the device 5-b.
Figure 6 shows the simulation results of the obtained potential profiles for both structures 5-a and 5-b.Due to symmetry, only half of these structures need to be simulated by considering the Neumann condition at the axis of symmetry.It has been verified that the analytical results are the same for all symmetric lines with respect to the axis of symmetry.Hence, for symmetrical devices the numerical effort can be significantly reduced by reducing the simulated space and therefore, reducing the total number of lines by a factor of two.
A two-dimensional metal oxide semiconductor (MOS) capacitor has been also simulated using the developed algorithm.The structure consists of a metallic gate, a Ptype silicon substrate and an oxide film layer as shown in Fig. 7.Because of the presence of a layer of silicon dioxide layer (which is a dielectric material) between the gate and the substrate, this device exhibits the properties of a capacitor.It is assumed in this analysis that the metal and the semiconductor work functions are equal, and that there is no electric field either within the oxide or at the oxidesilicon interface.Hence, the electric field is supposed to be zero whenever the applied voltage is zero.
Since the two-dimensional model of the considered device is symmetric, the investigation is limited to the left part of the device as illustrated in Fig. 7.When a positive voltage is applied at the electrode, free majority carriers (holes) are repelled by the induced electric field.Equilibrium condition requires that the potential is zero at the self-adjusting boundaries of the resulting fully-depleted layer and elsewhere in the neutral region.The lateral boundary conditions of the silicon and oxide layers are of Neumann type.
The discretization pattern is formed by 10 concentrat-ed equidistant lines crossing the curved area and 31 largely spaced lines elsewhere.The space charge layer is subdivided into 10 sub-layers to perform a stair-step approximation on the curved boundary.The dimensions used are: L = 35 m; l = 17 m; 0.1 < X o < 5.0 m; 30 < X s < 50 m; 10 14 < N A < 10 18 cm -3 ; 10 < V G < 100 (V G < V breakdown ).The program has been executed for different X o and N A and the results obtained are shown in Fig. 8.The results shown in this figure agree with both the previously stated boundaries and, the physical requirements.An important observation is that considering only a part of the neutral region can still save a lot of numerical efforts.This is because the potential within the whole neutral region is zero.This last point is particularly important when the left lateral side is far from the space charge region boundary, and consequently, the number of discretization lines can also be reduced.

Conclusions
The Method of Lines (MoL) is an efficient numerical tool for solving partial differential equations.Its strength lies in its assured convergence and its semi-analytical nature for both uniform and non-uniform discretization schemes.As it appears from the analysis done in this work, the MoL can be applied to electronic devices to find static potential profiles for different structures.A direct application of the MoL is to determine the breakdown voltages of power electronic devices needed for design prior to manufacturing.The developed algorithm allows significant saving in memory storage and computation time with assured convergence.This facilitates its extension to more sophisticated electronic structures as well as its integration into current computer-aided design tools.

Figure 1 .
Figure 1.Equidistant discretization pattern using the method of lines the matrix D of dimension (N+1xN) is a first order difference operator given by:

Figure 3 .
Figure 3. Cross sectional view of a two-layer planar semiconductor structure.(Dimensions are given in text)

Figure 6 .
Figure 5. Cross sectional view of two-layer semiconductor symmetrical structures

Figure 8 .
Figure 7. Two-dimensional model of the MOS capacitor used in the present analysis