JP2950358B2 - Simulation method of semiconductor device - Google Patents

Description

DETAILED DESCRIPTION OF THE INVENTION

[0001]

BACKGROUND OF THE INVENTION 1. Field of the Invention The present invention relates to a method for simulating a semiconductor device capable of accurately interpolating a potential in an entire space and accurately determining an electric field in the entire space.

[0002]

2. Description of the Related Art There is known a method for calculating an avalanche breakdown voltage at a junction of a semiconductor device by performing ionization integral calculation on a current path in the semiconductor device. That is, when the following expression (1) is satisfied in the calculation result of the ionization integral, avalanche breakdown of the junction occurs.

(Equation 1) In addition, in a space of two or more dimensions, an integration path is generally arbitrary, and when performing ionization integration, it is necessary to first determine a path for performing ionization integration. Generally, an electric field direction path is used as an ionization integration path. Incidentally, the ionization integration path is in the form of a continuous line, and its passing point is generally arbitrary. In order to execute the ionization integration, the electric fields at all points on the path must be defined. Further, the passing position of the ionization integration path is not known in advance, and is determined for the first time as a result of a path search following the electric field direction. Therefore, the electric field at an arbitrary position in the entire target space must be defined and clarified. On the other hand, in a general device simulation using a finite difference method using a rectangular grid, a directly obtained calculation result is a potential on the orthogonal grid point, and an electric field is obtained from this potential. According to the simplest method, the electric field calculated from the potential on the grid point is a horizontal electric field component Ex _{i + 1/2, j} along a horizontal grid line, as shown in FIG.
Ex _{i + i / 2, j + 1} (calculated from the potential on the left and right grid points) and the vertical electric field components Ey _{i, j + 1/2,} Ey along the vertical grid line
_{i + 1, j + 1/2} (calculated from potentials on upper and lower grid points).
Ψ in FIG. 1 is a potential on an orthogonal lattice point, and Ex _{i + 1/2, j} , E
y _{i, j + 1/2} is represented by the following equation. Ex _{i + 1/2, j} = − (ψ _{i + 1, j} −ψ _{i, j} ) / Δx _i Ey _{i, j + 1/2} = − (ψ _{i, j + 1} −ψ _{i, j} ) / Δy _j Δx _i = x _{i + 1} −x _i Δy _j = y _{j + 1} −y _j However, the electric field obtained by this method does not represent an electric field at a certain point, so that the calculation result is not practical. Therefore, a method of defining an electric field at the center point (center of gravity) of a rectangular lattice as shown in FIG. 1 has been performed by using linear interpolation of these electric field components. That is, Eo is the electric field at the position of the center of gravity, Eo
x component of Eo and _x, when the Eo _y and y components of _{Eo, Eo x = (Ex i} + 1/2, j + Ex i + 1/2, j + 1) / 2 Eo y = (Ey i, j _+1/2 + Ey _{i + 1, j + 1/2} ) / 2. In the definition of the electric field in the entire target space, the electric field defined for the center of gravity of the grid is regarded as the representative electric field of the entire grid, and the electric field is constant at all positions in the grid equal to the electric field of the previous center of gravity. And With this definition, it is possible to define an electric field in the entire target space, and it is possible to determine an ionization integration path obtained along the electric field direction. At the same time, the electric field on the path is determined, and the ionization integral can be calculated.

[0003]

However, since the conventional electric field definition lacks accuracy, the ionization integration path and the ionization integration value determined based on the definition are extremely low in reliability. That is, the x component of the electric field according to the above definition method is an average value of the x component electric field on the upper and lower boundaries of the lattice. Therefore, as long as the x-component electric fields at the upper and lower boundaries are not equal, the x-component of the electric field according to the conventionally defined method and the x-component electric field on the upper and lower boundaries are different. The electric field component on the lattice boundary is directly obtained from the lattice point potential at both ends, and is the most reliable information. However, in the above-described electric field definition method, this important information is lost on the boundary. In addition, the electric field is obtained by reversing the potential gradient (Ε = −grad 、), and the electric field is a constant electric field in the lattice in the above-described electric field definition method. Therefore, the potential surface in the rectangular lattice is one plane. Become. For each grid,
When the potential surface is set according to the above method, the above-mentioned inconsistency of the electric field components on the grid lines (lattice boundaries)
The potential surface between adjacent lattices is generally discontinuous as shown in FIG. That is, in the above-described electric field definition method, a discontinuous potential definition is inherent, and it is physically extremely unnatural. Therefore, the ionization integral calculation based on the electric field definition including such a contradiction is not reliable.

Accordingly, an object of the present invention is to provide a method of simulating a semiconductor device that can determine a more accurate ionization integration path by a novel electric field definition method that solves the above problem.

[0005]

According to the present invention, a method of simulating a semiconductor device using grid points scattered on a grid having polygonal grid cells in a two-dimensional space includes a grid point of one grid cell forming the grid. The process of obtaining each potential (ψ _a , ψ _b , ψ _c , ψ _d ) on (A, B, C, D), and the process of obtaining the potential (ψ _g ) of the center of gravity (G) of one grid cell , One grid cell (A, B, C, D) and the potential (ψ _a , ψ _b , ψ _c , ψ _d , ψ _g ) of the center of gravity (G) and the potential of the center of gravity (ψ _g ) As a common vertex and a plurality of triangular potential planes adjacent to each other (ψ _a ψ _b ψ _g , ψ _b ψ _c ψ _g , ψ
_c ψ _d ψ _g and ψ _a ψ _d ψ _g ) to interpolate the potential in one grid cell. In addition, a plurality of triangulation potential planes (ψ _a ψ _b ψ _g , ψ _b ψ _c ψ _g , ψ _c ψ _d ψ _g ,
The electric field in one grid cell is obtained from the slope of ψ _a ψ _d ψ _g ).

[0006]

The electric field component on the grid boundary coincides with the value directly obtained from the grid point potentials at both ends, so that the reliability is higher than in the conventional method, and the potential plane between adjacent grids is continuous. As a result, the ionization integration path can be determined accurately, and a more accurate breakdown voltage can be calculated.

[0007]

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS An embodiment of a semiconductor device simulation method according to the present invention will be described below with reference to FIGS.
In the electric field definition method according to the present embodiment, first, four points A (i, j + 1), B (i + 1, j + 1), C (C + 1) are arranged on an xy plane composed of an x-axis and a y-axis orthogonal to each other as in the conventional example.
(I + 1, j), D (i, j) potentials ψ _a , ψ _b ,
ψ _c and ψ _d are obtained. Here, the AB line segment and the CD line segment are parallel to the x-axis, and the AD line segment and the BC line segment are parallel to the y-axis. As shown in FIG. 3, the potentials ψ _a , ψ _b , ψ _c , ψ _d are xy
Take on the ψ axis perpendicular to the plane. Also, 4 on the orthogonal grid point
The center of gravity G (i of a square ABCD consisting of points A, B, C, D
Similarly, the potential ψ _{g of} +1/2, j + ／) is also on the ψ axis. Here, the potential ψ _g is an average of four grid point potentials, and ψ _g
= (Ψ _a , ψ _b , ψ _c , ψ _d ) / 4. Next, the potentials in the four grid cells (in the square ABCD) are represented by ψ _a , ψ _b ,
Interpolation is performed using a potential plane composed of ψ _c , ψ _d , and ψ _g . In the present embodiment, interpolation is performed using four triangular potential planes. That is, _a first triangular potential plane (1) composed of ψ _a , ψ _b , ψ _g , _a second triangular potential plane (2) composed of ψ _b , ψ _c , ψ _g , ψ _c , ψ _d , ψ _g The potential in the grid cell is interpolated and defined by _a third triangular potential plane (3) composed of the following, and a fourth triangular potential plane (4) composed of ψ _a , ψ _d and ψ _g . FIG. 3 shows the interpolation potential planes (1), (2), (3), and (4). By this potential definition, an electric field can be obtained in the entire space by calculating the gradient of each planar potential surface constituting the space, as shown in FIG. As understood from FIGS. 3 and 4, in this embodiment, the potential of the grid boundary is interpolated on the potential plane including the grid boundary as one side, and an electric field is obtained based on the interpolated potential. Ie AB,
The electric field components on the CD, AD, and BC line segments match the values directly obtained from the grid point potentials at both ends of each line segment, and provide highly reliable information unlike the conventional example. Since there is no inconsistency in the electric field components on the lattice boundaries, the potential planes of adjacent lattices are continuous. The electric field is obtained by reversing the sign of the potential gradient (−grad−). Therefore, the integration path is determined by searching the maximum inclination line direction of the potential surface in two directions, that is, the downstream side (electric field vector direction) and the upstream side (counter electric field vector direction) with the search start point as a starting point. . According to the above-described potential definition method, since a continuous potential can be defined for all regions, a reliable and accurate integration path can be obtained according to the definition. Further, in each triangular region, the potential surface is a plane, and the electric field vector is also constant in that region. Therefore, the integration path in each triangular area is a straight line (line segment), and the entire integration path is a continuous line segment connecting them. Incidentally, in the route search process following the maximum inclination line direction, a somewhat special situation may occur such as reaching the semiconductor boundary, reaching the valley on the downstream side, reaching the ridgeline on the upstream side, and the like. However, in any case, measures should be taken based on the principle of following the maximum slope of the potential surface under such circumstances.

Next, a specific electric field calculation method will be described. As described above, in the triangular region obtained by dividing the lattice into four parts, the electric field is constant with the potential surface as a plane. As in the present embodiment, when each grid point is a rectangular grid having grid lines parallel to the coordinate axes, an electric field can be easily obtained for both components. That is, in FIG. 5, the electric field in the lateral direction E _x, E _y the vertical electric _{field, ψ 1, ψ 2, ψ} 3, the potential on [psi ₄ grid points, [Delta] x, small distance between the potential Δy Then, the following equation is obtained.

(Equation 2) In addition, even when the grid is not such a rectangular grid, the electric field can be calculated as follows. That is, when two sides of the triangular potential plane are represented by the vectors shown in FIG. 6, the vectors indicating the normal direction of the plane are as follows.

(Equation 3) The side vectors a and b of the triangular potential surface have their starting points at the same position, and the vector b is within 180 degrees counterclockwise from the vector a when viewed from above the potential axis. The component of this normal vector

(Equation 4) Then, the electric field is expressed by the following equation.

(Equation 5) Each component of the normal vector is obtained by the following equation (FIG. 6). _{a x = x 2 -x 1,} a y = y 2 -y 1, a ψ = ψ 2 -ψ 1 (11) b x = x 3 -x 1, b y = y 3 -y 1, b ψ = ψ ₃ −ψ ₁ (12)

(Equation 6)

Next, a method of calculating the ionization integral based on the obtained ionization integral path and the electric field on the path will be described.
As shown in FIG. 7, on the ionization integration path, l is set in the direction from the N layer side to the P layer side, and the end point on the N layer side of the path is set to l = 0. Assuming that the extension of the ionization integration path is L, the end point on the P layer side is 1 = L. At this time, the ionization integral is as follows.

(Equation 7) Hereinafter, Equation (14) is used as a basic equation for ionization integral counting.
The ionization rates α _p and α _n are functions of the electric field strength.
A path is a set of continuous line segments, and the electric field is a constant value on each of these line segments. The l coordinate of each end point of these continuous line segments is denoted by l ₁ , l ₂ ,. . . , L _N (l ₁ = 0, l _N =
L), since the ionization rate is a function of the electric field, it is constant here on one line segment, and the ionization rates on the line segment between l _i and l _{i + 1} are α _{p (i)} and α _{n ( i)} . These are expressed by equation (14)
When applied to

(Equation 8) here,

(Equation 9) The coordinates of the end points of the route constituent line directly obtained are xy coordinates. The xy coordinates for _{l i (x i, y i} ) _{When, l i + 1 -l i =} {(x i + 1 -x i) 2 + (y i + 1 -y i) 2} 1 / ² (17). In this way, it is possible to write down the sum of ionization integrals (Σ) without introducing a new approximation method such as numerical integration.

Next, a method of obtaining a breakdown voltage from the ionization integral calculation will be described. In the case of a simulator that solves only the Poisson equation, the breakdown voltage is determined by repeating calculations with different trial voltages applied and searching for an applied voltage at which the ionization integral value becomes 1. That is, i = f (υ) (18) i: ionization integral value υ: applied voltage In a function of: i = i _BV (i _BV : ionization integral value at the time of junction breakdown, usually 1) It is equivalent to solving for υ. Here, the function f represents a numerical calculation part until a Poisson equation is solved for a given applied voltage to obtain an ionization integral value. In this embodiment, when solving this equation, the Newton method (Newton-Raphson method) is used, and the above equation i = f is used as a function to which the Newton method is applied.
Without using (υ) as it is, logi = g (υ) (g (υ) = log
[f (υ)]). Since the ionization integral value changes exponentially with respect to the applied voltage, the change in the logarithm of the ionization integral value with respect to the applied voltage becomes considerably linear. For this reason, by adopting the logarithm of the ionization integral value, it is possible to avoid an extreme increase in the next trial application voltage and ionization integral value for any trial application voltage, and eliminate the risk of divergence due to this problem. At the same time, since the Newton method is applied to a relationship close to a straight line, the approximation of the trial applied voltage during the repetition calculation is very high, and the number of repetitions required for convergence is also small.
It can be reduced more reliably than the case where the ionization integral value of the exponential change itself is used, and the calculation can be performed at extremely high speed. From the above, the actual solution of this part is as follows. logi = g (υ) (19) i: ionization integral value υ: applied voltage In Expression (19), assuming that i = i _BV (i _BV : ionization integral value at the time of junction breakdown), υ at this time is solved. When solving the equation h (x) = 0 by the Newton method, the derivative coefficient h ′ (x _n ) is required to determine the (n + 1) th trial value ( _xn : the nth trial value). However, the part corresponding to h (x) is not an analytic function, but includes a whole series of numerical calculation processes until a Poisson equation is solved in two dimensions and an ionization integral value is obtained. x _n ) cannot be obtained. Therefore, it is conceivable to obtain h ′ (x _n ) by a numerical calculation method. Basically, Δx is set to a small amount, and h ′ (x _n ) is calculated as h ′ (x _n ) ≒ (h (x _n + Δx) −h (x _n )).
/ Δx. This method makes it possible in principle to use the Newton method in this embodiment. However, in this method, to obtain a derivative of x _n a point x _n
And calculation of a function for two points x _n + Δx,
In calculating this function, it is necessary to solve a two-dimensional Poisson equation and perform a series of enormous calculations for calculating the ionization integral value. As a result, there is a problem that the calculation time is doubled by the differential coefficient calculator.

[0011] To solve the above problems, in the present _{embodiment, h (x n) '(} x n) h as an approximation of' ≒ (h (x _n)
−h (x _n−1 )) / (x _n −x _n−1 ) is used. According to this equation, a sufficient degree of approximation is shown at the stage of approaching convergence, and the function of this embodiment is close to a straight line even at the beginning of repetition, and is a sufficiently practical approximate value. Note that the trial applied voltage for the first time gives a withstand voltage that is expected by common sense. In addition, since the second trial applied voltage cannot be determined by the above-described method, it is given by expectation. However, in order to reduce the number of times of convergence, it is desirable to apply a voltage as close as possible to a breakdown voltage that is a solution. The above-described relationship between log (ionization integral value) and applied voltage is close to a straight line. In addition, as a result of examining the calculation results
It has been found that there is some correlation between the slope and the breakdown voltage. The relationship between the slope and the breakdown voltage in the vicinity of the breakdown voltage in several structures is shown by a straight line (broken line) in FIG. Using this correlation, the second trial applied voltage is determined as an approximate value. The relationship of a straight line as a representative of the correlation is represented by a = q (υ). Where a is the log at the time of surrender
(Ionization integral value) −Slope of applied voltage, υ is voltage. Now, the first trial applied voltage is υ ₁ , the ionization integral value for υ ₁ is i ₁ , the ionization integral value at the time of junction breakdown is i _BV, and the second trial applied voltage υ ₂ is determined as follows. a ₁ = q (υ ₁ ) (20)

(Equation 10) As can be understood from FIG. 9 showing this relationship, υ ₂ is log
(Ionization integral value) - the inclination relationship between applied voltage corresponds to the breakdown voltage of assuming a linear relationship a _1, a sufficiently valid approximation. Three FLRs (Field Limiting Ring)
When this embodiment was actually applied to a PN junction semiconductor device having the following formula, the calculation result of the breakdown voltage was 1632 volts, and the error from the measured value was 2% or less. FIG. 10 shows the potential distribution at the time of breakdown together with the ionization integration path where breakdown occurs.
FIG. 11 shows the electric field intensity distribution (main junction and the first FLR part enlarged).

According to this embodiment, the following effects can be obtained. (1) Since the electric field component on the grid boundary coincides with the value directly obtained from the grid point potential at both ends, reliability is higher than in the conventional method, and the potential plane between adjacent grids is continuous. As a result, an accurate ionization integration path can be determined, and a more accurate breakdown voltage can be calculated. (2) Unless the initial applied voltage deviates significantly from common sense,
Convergence is reached by repeating the number of trials required for convergence within 4 to 5 times, and the calculation time can be reduced.

[0013]

As described above, according to the present invention, the potential interpolation for the entire space and the electric field for the entire space can be accurately determined.

[Brief description of the drawings]

Fig. 1 Graphic showing simple electric field definition and electric field definition at center of gravity

FIG. 2 is a diagram showing discontinuities in adjacent lattice potential surfaces caused by a conventional method.

FIG. 3 is a diagram showing a potential surface

FIG. 4 is a diagram showing continuity of a potential surface according to the present embodiment.

FIG. 5 is a diagram showing a method of calculating an electric field in a triangular area of a rectangular grid parallel to a coordinate axis.

FIG. 6 is a diagram showing a normal vector of a triangular potential plane.

FIG. 7 is a diagram showing an ionization integration path.

FIG. 8 is a graph showing a relationship between a slope and a breakdown voltage near a breakdown voltage.

FIG. 9 is a graph showing a relationship between an ionization integral value and an applied voltage.

FIG. 10 shows a potential distribution at the time of breakdown.

FIG. 11 is a diagram showing an electric field intensity distribution.

[Explanation of symbols]

A, B, C, D... Orthogonal grid points, ψ _a , ψ _b , ψ _c , ψ _d
The potential at the lattice point, ψ _a ψ _b ψ _g , ψ _b ψ _c ψ _g , ψ _c ψ _d ψ _g ,
ψ _a ψ _d ψ _g

Claims (2)

(57) [Claims] 1. A method of simulating a semiconductor device using grid points scattered on a grid having polygonal grid cells in a two-dimensional space, wherein the grid points (A, B,
C, D) obtaining the potentials (ψ _a , ψ _b , ψ _c , ψ _d ), obtaining the potential (ψ g) of the center of gravity (G) of the one grid cell, and obtaining the one grid. Grid points (A, B, C,
D) and the potential of the center of gravity (G) (ψ _a , ψ _b , ψ _c , ψ _d ,
ψ _g ) and a plurality of triangular potential planes (ψ _a ψ _b ψ _g , ψ _b ψ _c ψ _g , ψ _c ψ _d ψ _g , adjacent to each other with the potential of the center of gravity (の_g ) as a common vertex ψ _a ψ _d一_g ) to interpolate the potential in the one grid cell. Wherein further, said plurality of triangular potential plane _{(ψ a ψ b ψ g,} ψ b ψ c ψ g, ψ c ψ d ψ g, ψ a ψ d ψ g) grating inclination from the one square The method for simulating a semiconductor device according to claim 1, wherein an electric field in the inside is obtained.