JP3001917B2 - Method for analyzing operation of semiconductor device and apparatus used therefor - Google Patents

Description

DETAILED DESCRIPTION OF THE INVENTION [Object of the Invention] (Field of Industrial Application) The present invention relates to a semiconductor device and modeling of a predetermined physical phenomenon.

(Prior Art) A basic equation used in semiconductor device modeling is usually expressed in the following form. In the following equation:
τ represents physical time.

∂p / ∂τ = (− 1 / q) (div J _p ) + GU (1) ∂n / ∂τ = (1 / q) (div J _n ) + GU (2) J _p = −qD _p (grad p) −qμ _p p (gradψ) (3) J _n = qD _n (grad n) −qμ _n n (gradψ) (4) div (gradψ) = (− q / ε) ( N _d −N _a + p−n) (5) These equations are commonly used in the conventional method and the method of the present invention. Equations (1) and (2) represent the continuity equation for holes and electrons, respectively. G, U
Represents generation and disappearance of carriers, respectively. Equation (3),
(4) specifically describes the hole current density and the electron current density, respectively. The hole current density _Jp and the electron current density
J _n is composed of a diffusion current component (first term) proportional to the carrier density gradients grad p and grad n, and an algebraic sum of a potential gradient gardψ and a drift component (second term) proportional to each carrier density. Represents Equation (5) is a Poisson equation representing the relationship between the potential ψ and the space charge density.

Hereinafter, the conventional solution will be specifically described. The first step in finding the numerical solutions of equations (1)-(5) is to rewrite these equations from a continuous form (differential equation) to a discrete system (difference equation). At this time, it is necessary to divide the device space to be analyzed into rectangular meshes. FIG. 9 conceptually shows mesh division for a two-dimensional space model. The vertical and horizontal division point intervals h _x , h _y are defined as shown in FIG. 9, and the distances h _x ′, h between the respective intermediate points are defined.
_y ′ is defined by the following equation.

_{h x '(M) = (} 1/2) [h x (M'-1) + h
_x (M ')] (6-1) _hy ' (N) = (1/2) [ _hy (N'-1) + h
_y (N ′)] (6-2) From equations (3), (6-1), and (6-2), in a discrete system, the divergence div J _p of the hole current can be expressed by the following approximate expression. .

div J _p = ∂J _px / ∂x + ∂J _py / ∂y [J _px (M ′) − J _px (M′−1)] / h _x ′ (M) + [J _py (N ′) − J _py (N′−1)] / h _y ′ (N) (7) Similar rewriting is possible for the electron current.
Since it is obvious, the description is omitted.

From the above, the result of rewriting the entire formulas (1) and (2) is as follows.

(1 / q) ｛[J _px (M ′) − J _px (M′−1)] / h _x ′ (M) + [J _py (N ′) − J _py (N′−1)] / h _y ′ (N)｝ = − U (M, N) (8-1) (1 / q) ｛[J _nx (M ′) − J _nx (M′−1)] / h _x ′ (M) + [ _Jny (N ') _-Jny (N'-1)] / _hy ' (N) = U (M, N) (8-2) where (8-1), (8-2) In the expression ()), the carrier generation term G is not required for normal device analysis, and is omitted for simplification. Also, for the sake of simplicity, it is assumed below that the calculation of a DC steady-state solution is performed, and the physical time differential terms ∂p / ∂τ and ∂n / ∂τ
Is zero.

However, it should be pointed out that even in the case of non-stationary problems where these physical time derivative terms are non-zero, the calculations can be performed without substantial difference.

Similarly, the Poisson equation (5) is expressed as follows.

[1 / h _x '(M)] ｛[ψ (M + 1, N) -ψ (M, N)] / h
_x (M ')-[ψ (M, N) -ψ (M-1, N)] / h _x (M'-
1)｝ + [1 / _hy '(N)] ｛[ψ (M, N + 1) -ψ (M, N)] / h
_y (N ')-[ψ (M, N) -ψ (M, N-1)] / _hy (N'-
1)｝ = − (q / ε) [Γ (M, N) + p (M, N) −n (M, N)] (9) However, for simplicity, the difference between the donor and acceptor impurity concentrations is put a Γ = N _d -N _a.

Now, for the following development, the carrier recombination term U
(M, N) is specifically given in the form of Shockley-Read-Hall. That, U (M, N) = [p (M, N) · n (M, N) + n i 2 (M, N)]
/ ｛Τ _n [p (M, N) + n _i (M, N)] + τ _p [n (M, N)
+ N _i (M, N)]｝ (10) Next, the current equations (3) and (4) are discretized. At this time, as is conventionally known, an approximate form of Scharfetter-Gummel is adopted in order to secure stability in numerical calculation. As a result, the x and y components of the hole and electron current densities are given as follows.

J _px (M ′) = [q / h _x (M ′)] × [λ _px1 (M ′) p (MN) + λ _px2 (M ′) p
(M + 1, N)] (11-1) _Jpy (N ') = [q / _hy (N')] x [? _Py1 (N ') p (MN) +? _Py2 (N') p
(M, N + 1)] ... (11-2) J nx (M ') = [q / h x (M')] × [λ nx1 (M ') n (MN) + λ nx2 (M') n
(M + 1, N)] (11-3) J _ny (N ′) = [q / _hy (N ′)] × [λ _ny1 (N ′) n (MN) + λ _ny2 (N ′) n
(M, N + 1)] (11-4) where λ in the above equation is given by the following equation.

λ _px1 (M ′) = [μ _p (M ′) / θ (M ′)] × [β (M ′) / (1-e− ^{β (M ′)} )] λ _px2 (M ′) = [μ _p (M ′) / θ (M ′)] × [β (M ′) / (1- _eβ ( ^{M ′)} )] λ _py1 (N ′) = [μ _p (N ′) / θ (N ′ )] × [β (N ′) / (1-e− ^{β (N ′)} )] λ _py2 (N ′) = [μ _p (N ′) / θ (N ′)] × [β (N ′) / (1− _eβ ^{(N ′)} )] λ _nx1 (M ′) = [μ _n (M ′) / μ _p (M ′)] × λ _px2 (M ′) λ _nx2 (M ′) = [μ _n (M ′) / μ _p (M ′)] × λ _px1 (M ′) λ _ny1 (N ′) = [μ _n (N ′) / μ _p (N ′)] × λ _py2 (N ′) λ by _{ny2 (N ') = [μ} n (N') / μ p (N ')] × λ py1 (N') ... (12) or more, write to the discrete system of the original equation (1) to (5) The replacement has been completed. That is, the equation we should solve is
Equations (8-1) and (8-2) obtained by substituting equations (11-1) to (11-4) for current, auxiliary relational equation (12) and equation (10) for recombination, (9), which consist of three equations with three unknowns p, n, ψ. However, looking at the expressions (8-1) and (8-2), since these include nonlinear relational expressions for the unknown quantities, it is difficult to find a solution as it is. Therefore, next, these equations are Taylor-expanded for unknown quantities, and a second-order or higher-order term is ignored to derive a linear equation for the changes Δp, δn, δψ of the unknown quantities.

First, it is assumed that each unknown quantity is composed of a trial value (represented by a superscript zero) and a sum of changes as follows.

p (M, N) = p ⁰ (M, N) + δp (M, N) n (M, N) = n ⁰ (M, N) + δn (M, N) ψ (M, N) = ψ ⁰ ( (M, N) + δψ (M, N) (13) If the current equations (11-1) to (11-4) are Taylor-expanded and second-order or higher terms are ignored, the following result is obtained.

J _px (M ') J _px ⁰ (M') + [∂J _px ⁰ (M ') / ∂p (M, N)] × δp (M,
N) + [∂J _px ⁰ (M ′) / ∂p (M + 1, N)] × δp
(M + 1, N) + [∂J _px ⁰ (M ′) / ∂ψ (M, N)] × δψ (M,
N) + [{J _px ⁰ (M ') / {(M + 1, N)] × δ}
(M + 1, N) (14-1) J _py (N ') J _py ⁰ (N') + [∂J _py ⁰ (N ') / ∂p (M, N)] × δp (M,
N) + [∂J _py ⁰ (N ′) / ∂p (M, N + 1)] × δp
(M, N + 1) + [{J _py ⁰ (N ') / {(M, N)] × δ} (M,
N) + [{J _py ⁰ (N ') / {(M, N + 1)] × δ}
(M, N + 1)… (14-2) J _nx (M ′) J _nx ⁰ (M ′) + [∂J _nx ⁰ (M ′) / ∂n (M, N)] × δn (M,
^{_{N) + [∂J nx 0 (}} M ') / ∂n (M + 1, N)] × δn
(M + 1, N) + [{J _nx ⁰ (M ') / {(M, N)] × δ} (M,
N) + [{J _nx ⁰ (M ') / {(M + 1, N)] × δ}
(M + 1, N) (14-3) J _ny (N ') J _ny ⁰ (N') + [∂J _ny ⁰ (N ') / ∂n (M, N)] × δn (M,
N) + [∂J _ny ⁰ (N ′) / ∂n (M, N + 1)] × δn
(M, N + 1) + [{J _ny ⁰ (N ') / {(M, N)] × δ} (M,
N) + [{J _ny ⁰ (N ') / {(M, N + 1)] × δ}
(M, N + 1)… (14-4) where the quantities with superscript zero are unknown trial values p ⁰ , n ⁰ ,
対 応 It shall correspond to ⁰ .

Since the recombination term U (M, N) is nonlinear in addition to the current, the same result is obtained by performing the same calculation.

U (M, N) = U 0 (M, N) + U p 0 (M, N) × δp (M, N) + U n 0 (M, N) δn (M, N) provided that U _p (M, N ) = ∂U (M, N) / ∂p (M, N) = [n (M, N) −τ _n U (M, N)] / ｛τ _n × [p
_{(M, N) + n i} (M, N)] + τ p [n (M, N) + n i (M, N)]} U n (M, N) = ∂U (M, N) / ∂p ( M, N) = [p (M, N) −τ _p U (M, N)] / ｛τ _n × [p
(M, N) + n _i (M, N)] + τ _p [n (M, N) + n _i (M, N)]｝ (15) where (13), (14-1), (14) -2), (14-
Substituting 3), (14-4) and (15) into (8-1), (8-2) and (9), and rearranging the terms, the following matrix / vector equation is obtained.

_{A T (M, N) Θ} T (M, N-1) + B T (M, N) Θ T (M, N) + C T (M, N) Θ T (M, N + 1) + D T (M, N _{) Θ T (M-1,} N) + E T (M, N) Θ T (M + 1, N) = F T (M, N), ... (16) where In addition, when each matrix and vector element are specifically written down,
It looks like this:

a ₁₁ = − [1 / qh _y ′ (N)] × J _py ⁰ (N−1) / Mp (M, N−1) a ₁₂ = 0 a ₁₃ = − [1 / qh _y ′ (N ^{_{)] × ∂J py 0 (N}} -1) / ∂ψ (M, N-1) a 21 = 0 a 22 = - [1 / qh y '(N)] × ∂J ny 0 (N-1) / ∂n (M, N-1 ) a 23 = - [1 / qh y '(N)] × ∂J ny 0 (N-1) / ∂ψ (M, N-1) a 31 = 0 a 32 _{= 0 a 33 = γ 1 (} M, N) = 1 / [h y '(N) h y (N-1)] b 11 = [1 / h x' (M)] × [∂J px 0 ( M) / ∂p (M, N) −∂J _px ⁰ (M−1) / ∂p (M, N)] + [1 / _hy ′ (N)] × [∂J _py ⁰ (N) / _{∂p (M, N) -∂J py} 0 (N-1) / ∂p (M, N)] + ∂U 0 (M, N) / ∂p (M, N) b 12 = ∂U 0 ( M, N) / ∂n (M, N) b ₁₃ = [1 / qh _x '(M)] × [∂J _px ⁰ (M) / ∂ψ (M, N) −∂J _px ⁰ (M− 1) / ∂ψ (M, N )] + [1 / qh x '(M)] × [∂J py 0 (N) / ∂ψ (M, N) ^{_{∂J py 0 (N-1)}} / ∂ψ (M, N)] b 21 = -∂U 0 (M, N) / ∂p (M, N) b 22 = [1 / qh x '(M) ^{_{] × [∂J nx 0 (M}} ) / ∂n (M, N) -∂J nx 0 (M-1) / ∂n (M, N)] + [1 / qh y '(N)] × [ ∂J _ny ⁰ (M) / ∂n (M, N) −∂J _ny ⁰ (N−1) / ∂n (M, N)] −∂U ⁰ (M, N) / ∂n (M, N ) B ₂₃ = [1 / qh _x '(M)] × [∂J _ny ⁰ (N) / ∂ψ (M, N) −∂J _nx ⁰ (M−1) / ∂ψ (M, N)] _{+ [1 / qh y '(} N)] × [∂J ny 0 (N) / ∂ψ (M, N) -∂J ny 0 (N-1) / ∂ψ (M, N)] b 31 = q / ε, b ₃₂ = −q / ε b ₃₃ = − [1 / h _x ′ (M)] [1 / h _x (M) + 1 / h _x (M−1)] − [1 / _hy '(N) + 1 / h y (N-1)] c 11 = [1 / qh y' (N)] × [∂J py 0 (N) / ∂p (M, N + 1)] c 12 = 0 c ₁₃ = [1 / qh _y '(N)] x [∂J _py ⁰ (N) / ∂ψ (M, N + 1)] c ₂₁ = 0 c ₂₂ = [1 / qh _y '(N)] × [∂J ny 0 (N) / ∂n (M, N + 1)] c 23 = [1 / qh y' (N)] × [∂J ny 0 (N) / ∂ψ (M , N + 1)] c ₃₁ = 0, c ₃₂ = 0 c ₃₃ = 1 / _hy ′ (N) _hy (N) d ₁₁ = − [1 / qh _x ′ (M)] × [∂J _px ⁰ ( M-1) / ∂p (M -1, N)] d 12 = 0 d 13 = - [1 / qh x '(M)] × [∂J px 0 (M-1) / ∂ψ (M- 1, N) d ₂₁ = 0 d ₂₂ = − [1 / qh _x ′ (M)] × [∂J _nx ⁰ (M−1) / ∂n (M−1, N) d ₂₃ = − [1 / qh _x ′ (M)] × [∂J _nx ⁰ (M−1) / ∂ψ (M−1, N) d ₃₁ = 0, d ₃₂ = 0 d ₃₃ = [1 / h _x ′ (M) h _{x (M-1) e 11} - [1 / qh x '(M)] × [∂J px 0 (M) / ∂p (M + 1, N)] e 12 = 0 e 13 = [1 / qh x' (M)] × [∂J _px ⁰ (M) / ∂ψ (M + 1, N)] e ₂₁ = 0 e ₂₂ = [1 / qh _x ′ (M)] × [∂J _nx ⁰ (M) / ∂ n (M + 1, N)] e ₂₃ = [1 / qh _x '(M)] × [ ∂J _nx ⁰ (M) / ∂ψ (M + 1, N)] e ₃₁ = 0, e ₃₂ = 0 e ₃₃ = 1 / h _x ′ (M) h _x (M) f ₁ = −U ⁰ (M, ^{_{n) + U p 0 (M}} , n) p 0 (M, n) + U n 0 (M, n) n 0 (M, n) - [1 / qh x '(M)] × [∂J p 0 ( M−1) / ∂ψ (M−1, N)] × [ψ ⁰ (M−1, N) −ψ ⁰ (M, N)] + [1 / qh _x ′ (M)] × [∂J ^{_{p 0 (M) / ∂ψ (}} M, N)] × [ψ 0 (M, N) -ψ 0 (M + 1, N)] - [1 / qh y '(N)] × [∂J p 0 ( N-1) / ∂ψ (M , N-1)] × [ψ 0 (M, N-1) -ψ 0 (M, N)] + [1 / qh y '(N)] × [∂J ^{_{p 0 (N) / ∂ψ (}} M, N)] × [ψ 0 (M, N) -ψ 0 (M, N + 1)] f 2 = -U 0 (M, N) -U p 0 (M, N) p ⁰ (M, N) −U _n ⁰ (M, N) n ⁰ (M, N) − [1 / qh _x ′ (M)] × [∂J _n ⁰ (M−1) / ∂ψ (M−1, N)] × [ψ ⁰ (M−1, N) −ψ ⁰ (M, N)] + [1 / qh _x ′ (M)] × [∂J _n ⁰ (M) / ∂ ψ (M, N ^{] × [ψ 0 (M,} N) -ψ 0 (M + 1, N)] - [1 / qh y '(N)] × [∂J n 0 (N-1) / ∂ψ (M, N-1 ^{) × [ψ 0 (M,} n-1) -ψ 0 (M, n)] + [1 / qh y '(n)] × [∂J n 0 (n) / ∂ψ (M, n)] × [ψ ⁰ (M, N) −ψ ⁰ (M, N + 1)] f ₃ = − (q / ε) Γ (M, N) (18) The nonlinear equation (8-1) discretized as described above ,
As a result of linearizing (8-2) and (9), it was found that a matrix / vector equation (16) was derived.

Since equation (16) is equivalent to one of a large number of division points arranged two-dimensionally, in order to obtain a solution for the entire device plane, a combination of equation (16) for all division points is used. I have to solve it. Assuming that the division point numbers corresponding to the entire device are in the range of 1 ≦ M ≦ K, 1 ≦ N ≦ L, the result of rewriting equation (16) as a whole is
As shown in Figure 10. However, for simplicity, A _T (M,
N) was rewritten as A _MN , _{Ｔ T} (M, N) as Θ _MN , and F _T (M, N) as F _MN .

In Figure 10, each of the matrix such as B ₁₁ is because having a dimension of (3 × 3), the dimensions of the whole of the matrix becomes (3KL) × (3KL). Assuming that K = 30, L = 40 is given as a standard device model, the dimension of the entire matrix is 3,600 × 3,600 = 12.960,000 = 1.
The size is 296 × 10 ⁷ .

There are various methods for actually solving the matrix / vector equation in FIG. 10, and the most commonly used method is Gaussian elimination. Today, with the development of computers, high-speed calculation of large-capacity memories has become possible. Therefore, the problem of the magnitude of the above numerical example can be solved by using an erasure method. Other methods include various iterative methods, which are effective when the matrix size is very large.

The matrix solution method is widely applied not only to device modeling but also to other physical and engineering problems, and is a well-known calculation process. Therefore, it is assumed that the solution is obtained with a predetermined calculation time. That is, it is assumed that the solution of equation (16) is obtained for the entire space of the given device. This means that the nonlinear equations (8-1), (8-
It is inevitable that a solution was obtained by linearly approximating 2) and (9). However, when this solution is substituted into a nonlinear equation, both sides are generally not equal due to a nonlinear effect. Then, [p (M, N), n
(M, N), ψ (M, N), 1 ≦ M ≦ K, 1 ≦ N ≦ L]
⁰ (M, N), n ⁰ (M, N), ψ ⁰ (M, N), 1 ≦ M ≦ K, 1 ≦ N ≦
L] and execute the same calculation process as described above again. In the same manner, if the iterative calculation is performed the required number of times, the nonlinear equations (8-1) and (8-
A solution that completely satisfies 2) and (9) is obtained.

(Problems to be Solved by the Invention) In the matrix solution described above, since the performance of today's computers has been improved, it has become possible to increase the memory capacity and perform high-speed calculations.
This can be applied to fairly large problems. However, the problem to be solved becomes larger with the improvement of the performance of the computer. Therefore, speeding up the matrix solution is still one of the biggest problems in performing the modeling.

In order to solve this problem, the present invention provides a method and an apparatus for solving a predetermined equation without applying a matrix solution, analyzing the operation of a semiconductor device, and analyzing a predetermined physical phenomenon, without applying a matrix solution. Is provided.

[Constitution of the Invention] (Means for Solving the Problems) The first method for analyzing the operation of a semiconductor device according to the present invention is a method for analyzing a system of simultaneous equations composed of electron and hole transport equations and Poisson equation in semiconductor device modeling. In the method of analyzing using the system, the system of simultaneous equations f _p = 0, f _n = 0, f _ψ = 0 is converted into an artificial time differential term dp / dt, dn /
dt, d [phi] / dt and sensitivity coefficients λ _p, λ _n, below which includes a lambda _[psi (a
1), (b1), ( c1) replacing the _{formula, dp / dt = -λ p f} p ... (a1) dn / dt = λ n f n ... (b1) dψ / dt = λ ψ f ψ ... (c1 The splitting point arrangement (M, N) of the semiconductor device is determined, and the sensitivity coefficients λ _p , λ _n , λ _を are also given appropriate spatial position dependency that matches the physical characteristics of the semiconductor device. Equations (a1), (b1), and (c1) are replaced by the following equations (d), (e),
(F), dp (M, N) / dt = −λ _p (M, N) f _p (M, N) (d) dn (M, N) / dt = λ _n (M, N) f _n (M, N)… (e) dM (M, N) / dt = λ _ψ (M, N) f _ψ (M, N)… (f) where f _p , f on the right side of each equation is _n, f _[psi defined as _{follows, f p (M, n)} = (1 / q) {[J px (M) -J px (M-1)] / h x '(M)} + (1 / q) {[J py (n) -J py (n-1)] / h y '(M)} + U (M, n) ... (22-1) f n (M, n) = ( 1 / q) ｛[J _nx (M) −J _nx (M−1)] / h _x ′ (M) + (1 / q) ｛[J _ny (N) −J _ny (N−1)] / h _y ′ (M)｝ + U (M, N) (22-2) f _ψ (M, N) = [1 / h _x ′ (M)] × ｛[ψ (M + 1, N) −ψ (M , N)] / h _x (M) − [{(M, N) −ψ (M−1, N)] / h _x (M−1)} + [1 / _hy ′ (N)] × ｛ [Ψ (M, N + 1) −ψ (M, N)] / _hy (N) − [ψ (M, N) −ψ (M, N−1)] / _hy ( N−1)｝ + (q / ε) [Γ (M, N) + p (M, N) −n (M, N)] (22-3) The above (d), (e), (f) A solution of the system of simultaneous equations is obtained by time-integrating the equation.

However, in the above description, for simplicity, the original function
The DC steady-state problem in which f _p and f _n do not include the derivative with respect to the physical time τ is taken up. In the case where these problems are included, in equations (8-1) and (8-2), ∂p / ∂τ and ∂n /
The respective differential equation ∂τ [p (M, N) -p 0 (M, N)] / ∂τ [n (M, N) -n 0 (M, N)] are replaced in the form of / ∂τ It should be pointed out that the solution of the system of simultaneous equations can be obtained without departing from the principle of the present calculation method.

A second operation analysis method of a semiconductor device according to the present invention is a method of analyzing a Poisson equation representing a relationship between an electric potential and a space charge in a semiconductor device modeling by a computer.
Replaced by the following (a2) expression containing ∂t and sensitivity coefficient _{λ, ∂ψ / ∂t = λf ψ} ... (a2) dividing point arrangement of semiconductor devices (M, N) and determines, said in sensitivity coefficient lambda Given the appropriate spatial position dependency that matches the physical characteristics of the semiconductor device, the above equation (a2) is converted into the following equation (b2), and 、 (M, N) / ∂t = λ (M, N) _fψ (M, N) (b2) The solution of the Poisson equation is obtained by time-integrating the above equation (b2).

A third method for analyzing the operation of a semiconductor device according to the present invention is a method for analyzing an equivalent circuit constant determination problem of a semiconductor device based on a least squares method by a computer, wherein the equivalent circuit constant determination problem includes a time differential term dq / dt and a sensitivity. The following equation (a3) including the coefficient λ is substituted, and dq / dt = −λ (Aq−b) (a3) By giving a value suitable for the physical characteristic of each equivalent circuit constant to the sensitivity coefficient λ, The solution of the problem of determining the equivalent circuit constant is obtained by time-integrating the equation (a3).

First analysis device of the present invention, in a semiconductor device modeling, electron, hole transport equation and the system of simultaneous equations consisting of the Poisson equation time derivative term and differential equations and a sensitivity coefficient dp / dt = -λ _p f _p , dn / dt = λ _n f _n , dψ / dt =
An apparatus for analyzing replaced with lambda _[psi f _[psi, the device includes a single or plurality of m logic cells, the m amplifiers the bound to each logic cell having an adjustable amplification gain lambda _i or therewith It is characterized by comprising an equivalent multiplier and detecting means provided in a network for each of the logic cells and detecting that the transient phenomenon in the device has reached a steady state.

Second analyzing device of the present invention, there is provided an apparatus for analyzing by replacing the differential equation _dψ / dt = λf ψ Poisson equation representing the relationship between the potential and space charge in the semiconductor device modeling, the apparatus comprising one or more m Logic cells, m amplifiers or multipliers equivalent to the amplifiers having an adjustable amplification gain λ _i coupled to each of the logic cells, and provided in a network for each of the logic cells, And detecting means for detecting that the transient phenomenon has become a steady state.

The third analysis apparatus of the present invention solves the problem of determining the equivalent circuit constant of a semiconductor device based on the least squares method using a differential equation dq / dt =
An apparatus for analyzing replaced with -λ (Aq-b), the apparatus, the m amplifiers having one or a plurality of m logic cells, the bound adjustable amplification gain lambda _i to each logic cell Alternatively, it is characterized by comprising a multiplier equivalent to this, and detecting means provided in a network form for each of the logic cells and detecting that the transient phenomenon in the device has become a steady state.

(Operation) In the method of the present invention, a method of solving a system of simultaneous equations composed of electron and hole transport equations and Poisson equation in semiconductor device modeling, a method of solving Poisson equation representing a relationship between potential and space charge in semiconductor device modeling, In solving the equivalent circuit constant determination problem of a semiconductor device based on the square method, each of the above equations is replaced with a differential equation having a time differential term and a sensitivity coefficient λ, and the sensitivity coefficient λ
Is set in accordance with the physical characteristics of the equation to be solved, so that the matrix solution can be executed at high speed even if the target equation is further enlarged.

The device of the present invention regards the independent variable t of the following equation dx / dt = λf (x) (a4) as time, and the unknown quantity x is represented by the voltage of a predetermined terminal of the electric network satisfying the equation (a4). Shall be If the transient state of the network progresses with time and a steady state is realized, the entire apparatus is configured so as to give the solution of equation (a4).

In general, in the digital calculation based on the conventional method, in the calculation of the solution of the matrix equation Ax = b, when the elimination method is executed, the so-called inverse matrix calculation or practically written as x = A ⁻¹ b is used. Equivalent calculations are involved, and performing this requires a great deal of calculation time.

In contrast, in equation (a4) used in the present invention,
Does not include calculation of the inverse matrix. Simply the function value x ^{o at} some current time
Using the matrix A, the multiplication Ax ^o of the matrix A and the subtraction for subtracting the constant vector b from the matrix A, and further, for the obtained Ax ^o −b, the time step length δt from the left The value obtained by multiplying the coefficient matrix λ is equal to the magnitude of the time differential determined by the connection of the network. This calculation process, as it were, consists only of calculations in the forward direction.

The advantage of forward computation is especially pronounced when matrix A contains a large number of zero elements. In this case, multiply Ax
^{This is because} only the non-zero elements of A need to be considered when executing ^o . There is one condition for the present calculation method to give a true solution of the matrix equation Ax = b. That is, when the integral calculation on the time axis of the differential equation (a4) is executed for a sufficiently long time, the unknown quantity x reaches a steady value. If you write this in the formula, Becomes Assuming that this condition is satisfied, the left side of the differential equation (a4) becomes zero after a sufficient time has elapsed. Therefore, since the right side must be zero,
The obtained x satisfies Ax = b.

In the above-described calculation execution process, whether or not x can reach a steady value is determined by the suitability of selecting the diagonal coefficient row example λ as the sensitivity coefficient. In general, λ is too small, since the time differential quantity dx / dt too small, the change from the initial trial value x = x _o is not Okoshira little, therefore stationary solution it takes a lot of time to obtain. On the other hand, if λ is too large, the amount of time differentiation is too large and the differential equation (a4)
The behavior of the entire system of equations (or network) represented by Looking at an example of actual numerical calculation by digital simulation, if λ is too large, the change of unknown quantity x will be too large, and after a short period of time, the change in the opposite direction to correct it will be corrected. Happens.

If this is repeated with the progress of time, it turns out that a vibration phenomenon on the time axis eventually occurs, and the unknown quantity x cannot reach a steady value forever. It is also known that if the unknown quantity change is further large, a divergence phenomenon such as the amplitude of the vibration becoming infinitely large or the unknown quantity value becoming infinitely large in one direction (either positive or negative) is caused. ing.

To summarize the above, in order to properly operate the apparatus of the present invention by the calculation method of the present invention, the magnitude of the diagonal coefficient matrix λ is determined to an appropriate value so that the change of the unknown quantity of a moderate magnitude is obtained. It is important to make it happen.

More specifically, since λ is a diagonal coefficient matrix, n elements (x ₁ , x ₂ ,.
Note that _n coefficients (λ ₁ , λ ₂ ,... λ _n ) correspond to x _n ). For this reason, the calculation method and apparatus of the present invention have a degree of freedom in which n coefficients can be individually adjusted for n unknown quantity elements, and each of them can be determined to an optimum value. In this regard, it is much more efficient than the iterative method conventionally known as a matrix solution (in which a so-called relaxation coefficient is introduced, but the number is usually one for the whole matrix equation), The amount can be brought to a steady-state value.

Next, the optimization method of λ will be discussed. λ is a diagonal matrix as described above, and has n elements λ ₁ , λ ₂ ,... λ _n . That is, λ = diag (λ ₁ , λ ₂ ,... Λ _n ) (19). Here, in equation (2), the unknowns x _i , i = 1,2,
.., N, there must be a maximum absolute value of the time differential amount | dx _i / dt | _max so that the entire device system can operate stably. Therefore, the component λ _i of λ is The value must not be exceeded. If this is expressed by the formula, Becomes Since the magnitude of the right side of equation (20) depends on the nature of the problem, it seems difficult to discuss it more specifically within the scope of general theory. According to matrix algebra, this kind of decision condition is usually associated with the eigenvalues of the matrix, but solving a real problem and finding the eigenvalues of the matrix are problems that are almost the same and almost as difficult. Therefore, it is rare that practically effective determination conditions are given. Therefore, when actually assembled device of the present invention, the amplifier associated with each of the n logic cells, as will be described later is provided, each of which configured to have an amplification gain lambda _i. Each amplifier shall be designed to have a sufficiently large amplification gain to sufficiently address the individual practical problem, and shall include a knob or controller for gain adjustment.

Hereinafter, the present invention will be described with reference to the drawings.

FIG. 1 is a flowchart for explaining an operation analysis method of a semiconductor device according to a first embodiment of the present invention.

In the present invention, equations (8-1), (8-2), (9)
To find the solution of, replace them with the following form, including the time derivative term.

dp (M, N) / dt = −λ _p (M, N) f _p (M, N) (21-1) dn (M, N) / dt = λ _n (M, N) f _n (M , N)… (21-2) dψ (M, N) / dt = λ _ψ (M, N) f _ψ (M, N) (21-3) where f _p , f _n , f _ψ is defined as follows.

f _p (M, N) = (1 / q) ｛[J _px (M) −J _px (M−1)] / h _x ′ (M)｝ + (1 / q) ｛[J _py (N) −J _py (N−1)] / h _y ′ (M)｝ + U (M, N) (22-1) f _n (M, N) = (1 / q) ｛[J _nx (M) − _{J nx (M-1)]} / h x '(M) + (1 / q) {[J ny (N) -J ny (N-1)] / h y' (M)} -U (M, _{N) ... (22-2) f ψ} (M, N) = [1 / h x '(M)] × {[ψ (M + 1, N) -ψ (M, N)] / h x (M) - [Ψ (M, N) −ψ (M−1, N)] / h _x (M−1)｝ + [1 / _hy ′ (N)] × ｛[ψ (M, N + 1) −ψ (M , N)] / _hy (N)-[{(M, N)-{(M, N-1)] / _hy (N-1)} + (q / [epsilon]) [{(M, N) + P (M, N) -n (M, N)] (22-3) Returning to (21-1) to (21-3) again, the time differential term on the left side of each equation is approximated by difference. At this time, a finite number of division points are provided on the time axis, and the value of the old time is represented by a superscript numeral 0 and the value of the new time is represented by a superscript numeral 1. As a result, the following equations are obtained.

p ¹ (M, N) = p ⁰ (M, N) −δtλ _p (M, N) f _p ⁰ (M, N) (23-1) n ¹ (M, N) = n ⁰ (M, N n) + δtλn (M, n ) f n 0 (M, n) ... (23-2) ψ 1 (M, n) = ψ 0 (M, n) + δtλ ψ (M, n) f ψ 0 (M, N) (23-3) where 1 ≦ M ≦ K and 1 ≦ N ≦ L.

Hereinafter, description will be made with reference to the flowchart of FIG.
In this method, when starting the calculation, in step a1, the structure and impurity data of the semiconductor device are input to the memory of the computer in step a1. Next, in step b1, the arrangement of the dividing points of the semiconductor device is determined, and the bias voltage V is input in step c1. Then the basic variables p (M, N), n (M, N), ψ (M, N) initial trial value p ⁰ of ^{(M, N), n 0} (M, N), ψ 0 (M, N ) Is given in step d1, and the initial time for starting time integration in step e1
defining the t ^1.

Next, in step f1, equations (23-1) to (23-3)
The values on the right side of each expression are calculated for all points 1 ≦ M ≦ K, 1 ≦
For N ≦ L, these are calculated as p ¹ (M, N), n ¹ (M, N),
ψ ¹ (M, N); 1 ≦ M ≦ K, 1 ≦ N ≦ L. This means that the integration on the time axis has advanced by one step.

Next, in step g1, the absolute value of the variation of each variable | p
¹ (M, N) −p ⁰ (M, N) |, | n ¹ (M, N) −n ⁰ (M, N) |, | ψ ¹
Find (M, N) −０ ⁰ (M, N) | for all points and check if all of them are less than a predetermined error limit.

If at least one point in the device space has a change larger than the error limit, in step h1, the correction values p ¹ (M, N), n ¹ (M, N), ψ ¹ obtained earlier. (M, N) are regarded as new trial values, and these are p ⁰ (M, N), n ⁰ (M, N), ^{０ 0}
(M, N) placed equal to, further time t ¹ as t ^0, performing integration calculation on further time axis in the same manner as described above.

If the absolute value of each variable change falls below the error limit at all points, the calculation is terminated assuming that a steady solution has been obtained. The resulting solution is (21-1)-(21
In -3), since the time differential term on the left side of each equation becomes zero, f _p (M, N) = 0, f _n (M, N) = 0, f
_ψ (M, N) = 0; 1 ≦ M ≦ K, 1 ≦ N ≦ L holds. This means that the equations (8-1), (8-2) and (9) have been solved.

The above is the outline of the present embodiment. The features in comparison with the conventional method are that the present embodiment does not include the process of linearizing the nonlinear equation, and the equation (16) shown in the description of the conventional method and FIG. There is no need to calculate the solution of the large-sized matrix / vector equations, and therefore, without performing complicated calculations using a large memory capacity, simply by using equations (23-1) to (23-3). It is only necessary to determine the correction value of the basic variate by finding the value on the right side of the equation, and the coefficients λ _p (M, N) and λ _n (M, N) as adjustment factors
And λ _ψ (M, N) _-hereinafter simply referred to as sensitivity coefficients--can take different values for each point in the device space, so by giving them appropriate values that match the physical properties That the convergence of the solution can be achieved efficiently,
The process of calculating the solutions of the equations (21-1) to (21-3) in this embodiment is nothing but a numerical integration calculation on the time axis. It is of course possible to perform this calculation by digital calculation using a large-sized computer as usual, but in addition, by applying an electronic circuit, the time response (transient phenomenon) of the circuit can be obtained.
Once has reached a steady state, it is possible to construct a device whose state corresponds to the solution of a given equation. This point will be described below as an apparatus of the present invention.

As described above, the method of the first embodiment is completely different from the conventional method based on the matrix solution, and has a possibility of opening up a new field in device modeling.

Next, an operation analysis method of a semiconductor device according to a second embodiment of the present invention will be described with reference to FIG.

In device modeling, for example, when the p-side of the pn diode is kept at zero potential and the n-type is kept at positive potential V, a depletion layer is formed near the interface of the pn junction. At this time, if the potential V is equal to or lower than the breakdown voltage of the diode, the current flowing through the diode is small enough to be ignored. In such a case,
Porin equation (5) among the basic equations (1) to (5)
It is known that the distribution of the potential ψ of the entire device can be accurately obtained by solving only the above (for example, Kurata, “Numerical A
nalysis for Semiconductor Devices, Heath, 1982, No.8
Chapters or Kurata, "Behavior Theory of Bipolar Transistors", Kindai Kagaku-sha, Sections 3.2 and 5.1, 1980).
At this time, the hole density and the electron density p, n are respectively given by the following equations.

Substituting ^{_{p = n i e -θψ (24-1}} ) n = n i e θ (ψ-V) (24-2) them in equation (5), the following equation when a two-dimensional space is obtained .

^{(∂ 2 ψ / ∂x 2)} + (∂ 2 ψ / ∂y 2) = - (q / ε) [Γ (x, y) + n i e -θψ -n i e θ (ψ-V)] ... (25) Converting this equation to a degenerate system is performed in the same manner as above. The specific shape is (24) on the right side of equation (9).
-1) and (24-2). That is, [1 / h _x '(M)] × ｛[ψ (M + 1, N) −ψ (M, N)] / h _x (M) − [ψ (M, N) −ψ (M−1, N)] / h _x (M−1)｝ + [1 / _hy ′ (N)] × ｛[ψ (M, N + 1) −ψ (M, N)] / h _y (N) − [ψ ( M, N) − {(M, N−1)] / _hy (N−1)} = − ρ (M, N) / ε (26) where ρ (M, N) = q [Γ (M , n) is ^{_{+ n i e -θψ (M,}} n) -n i e θ (ψ (M, n) -V)] ... (27). Solving this equation by conventional methods is
Following the same procedure as detailed above, it results in calculating the solution of the matrix vector equation and performing iterative calculations.

In contrast, the method of the present invention starts with the following equation including the time derivative term and the sensitivity coefficient λ.

∂ψ (M, N) / ∂t = λ (M, N) × << [1 / h _x ′ (M)] × ｛[ψ (M + 1, N) −ψ (M, N)] / h _x ( M ′) − [ψ (M, N) −ψ (M−1, N)] / h _x (M′−1)} + [1 / _hy ′ (N)] × ｛[ψ (M, N + 1 ) − {(M, N)] / _hy (N ′) − [{(M, N) − {(M, N−1)] / _hy (M′−1)} + ρ (M, N) / Ε >> (28) The procedure for integrating this equation along the time axis is executed in the same manner as already described in the description of the first embodiment. FIG. 2 shows a flowchart for this purpose. That is,
In this method, when starting the calculation, in step a2, the structure and impurity data of the semiconductor device are input to the memory of the computer. Next, in step b2, the division point arrangement of the semiconductor device is determined, and the bias voltage V is input in step c2. Next, an initial trial value of the potential ψ ⁰ is given in step d2, and
defining the initial time t ¹ for the time integration start in e2.

In step f2, to calculate the [psi ^1. [psi ¹ is given by the following equation.

ψ ¹ (M, N) = ψ ⁰ (MN) + δt · λ (M, N) × << [1 / h _x ′ (M)] × ｛[ψ ⁰ (M + 1, N) −ψ ⁰ (M, N )] / H _x (M ′) − [ψ ⁰ (M, N) −ψ ⁰ (M−1, N)] / h _x (M−
1)｝ + [1 / _hy ′ (N)] × ｛[ψ ⁰ (M, N + 1)-０ ⁰ (M, N)] / _hy (N ')-[ψ ⁰ (M, N)- ^{０ 0} (M, N−1)] / h _y (M′−
1)｝ + ρ (M, N) / ε >> (29) According to equation (29), the value on the right side is changed to all points 1 ≦ M ≦ K,
≦ N ≦ L, and these are calculated as ^{１ 1} (M, N); 1 ≦ M ≦
Let K, 1 ≦ N ≦ L. This means that the integration on the time axis has advanced by one step.

Next, in step g2, the absolute value of the change of the variable | ψ
^{1 (M, N) -ψ 0} (M, N) | look for all points, all of which checks whether less than a predetermined error limit.

If at least one point in the device space has a change larger than the error limit, in step h2, the previously determined correction value ^{１ 1} (M, N) is regarded as a new trial value, and these are set to ψ. ^{0 (M,} N) placed equal to, further time t ¹ as t ^0, performing integration calculation on further time axis in the same manner as described above.

If the absolute value of the variate change at all points falls below the error limit, the calculation is terminated assuming that a steady-state solution has been obtained. The resulting solution is such that in (21-3), the time derivative term on the left side of the equation becomes zero, so that _fψ (M, N) = 0; 1 ≦ M ≦ K, ≦ N ≦ L holds. This means that equation (9) or equation (2
It is nothing but the solution of 6).

Note that when performing integration on the time axis in the above,
Absolute value of change of unknown variate ψ (M, N) | ψ ¹ (M, N) −ψ ⁰
If (M, N) | is too small, the stability of the calculation process is good, but 時間 (M, N) cannot easily reach the true solution even after the lapse of time t. On the other hand, if the absolute value of the above change is too large, the calculation process becomes unstable, and oscillation or divergence of ψ (M, N) occurs with the passage of time t, so that a true solution cannot be obtained. Absent.

The sensitivity coefficient λ (M, N) is an adjustment factor for appropriately causing such a change in the variable. In particular, in the case of the Poisson equation of a semiconductor device as in this example and the two-dimensional space model shown in equation (29), the following equation is used as λ (M, N) satisfying the above condition. It is known to be good. That is, first, the following quantities are defined in relation to equation (29).

_{H = 1 / [h x (} M'-1) h x '(M)] + 1 / [h x (M') h x '(M)] + 1 / [h y (N'-1) h y' (N)] + 1 / [ h y (N ') h y' (N)] ... (30-1) α = (qθn i / ε) × [exp {-θψ 0 (M, N)} + exp {θ ({ ⁰ (M, N) -V)}] (30-2) λ (M, N) · δt ≦ 2 / (2H + α) (30-3)

Hereinafter, the meaning of Expression (30-3) will be qualitatively described. In general, in a region having a high impurity concentration inside a semiconductor device, one of the electron and hole densities n and p is substantially equal to the impurity concentration, and the space charge neutral condition is maintained. In this case α is very large compared to H and 2 of the molecule, so λ
(M, N) · δt is very small. As a result, the variate change |
ψ ¹ (M, N) −ψ ⁰ (M, N) | is also very small. On the other hand, in a region where the impurity concentration is relatively low and a depletion layer is easily formed, α is equal to H or negligibly smaller than H. As a result, λ (M, N) · δt
Is relatively large and therefore the variate change is also relatively large.

For the above reasons, if a value that satisfies the space charge neutral condition everywhere in the device is given to the variable ψ (M, N) as an initial condition for starting the calculation, In the low-concentration region, a large variable change occurs, which promotes the formation of a depletion layer, and the variable ψ (M, N) converges toward a true solution.

Next, an apparatus for executing the method of the second embodiment will be described with reference to FIG. To assemble such a device, an amplifier associated with each of the N _T logic cells is configured to have an amplification gain λ _i . Here, _NT is the total number of mesh division points when the original equation is replaced with a discrete system based on an approximate solution method such as the difference method or the finite element method in order to solve the original equation belonging to the continuous system using a computer. It is. Each amplifier is designed to have a large amplification gain enough to cope with individual practical problems, and includes a knob or controller for gain adjustment. However, the amplification gain may be smaller than 1 due to a problem.
In this case, the amplifier is actually an attenuator, but such a case is also included. Alternatively, as an apparatus equivalent to an amplifier, a multiplier for obtaining a product of an input amount and an amplification gain may be used, and this case is also included.

Reference numerals 1 to 12 are logic cells, which have a non-linear input / output response. In this embodiment, _NT is 12. Amplifiers a1, a2,... A12 are provided at the input terminals of each cell. Also, the node between each amplifier and each cell is
There are provided resonance circuits b1, b2,... B12 each including a resistor and a capacitor. Each amplifier has a variable gain corresponding to λ. A current is input from each of the input current sources c1, c2,... C12 to each cell, and the fact that this input current has reached a steady state in the apparatus of the present invention is provided in a network manner for each cell. , 1212. The detection result is addressed by the address unit 100, converted into a digital signal by the D / A converter 101, and sent to the circuit 103. The circuit 103 includes a differentiator, a multiplier, and an adder, executes a predetermined operation, and outputs measurement data to the circuit 104.

In the above embodiment, the total number of mesh division points is equal to _NT .
Although an example using _NT cells has been shown, if a sweeper is used, the analyzer of the present invention can be configured using m cells equal to or smaller than _NT , and the same operation as in the above embodiment can be performed. it can.

Next, with reference to FIGS. 4A, 4B, and 5, the problem of determining the equivalent circuit constant of the semiconductor device by the least square method according to the third embodiment of the present invention will be described. Now, four y parameters (all of which are complex admittances) y ₁₁ , y ₁₂ , y ₂₁ , y
₂₂ are measured at N frequencies f _k (k = 1, 2,..., N), and this is measured by an equivalent circuit (the first circuit) including _M equivalent circuit constants q ₁ , q ₂ ,. If fitting is attempted using (see Figures 4A and 4B), the following matrix equation is derived.

Aq = b, A (a _ij ), b = (b _i ) (31-1) However, in the above, α _mki = ∂y _mk / ∂q _i (32-1) Here, an asterisk (*) represents a conjugate complex number, and m takes a value of 1 or 2.

The above matrix equation is calculated based on the solution of the present invention by dq / dt = −
Replace as λ (Aq-b). This equation is solved based on the flowchart of FIG. In step a3, measured data of the y parameter: η _11k , η _12k , η _21k , η
_22k ; read k = 1,2, ... N. Then the initial trial value q ₁ ⁰ unknown quantity at step b3, q ₂ ^0, ..., give q _M ^0. Step c3
, An initial time and a time interval are given, and time integration is started. In step d3, the theoretical value y _mk of the y parameter
_Mki , β using the measured value of the y parameter and η _mk
seeking _mk, seeking a _ij, b _i in step e3.

Next, in step f3, q _i ¹ is calculated for i = 1, 2,..., M. q _i ¹ is given by the following equation.

q _i ¹ = q _i ⁰ −δt · λ _i (Aq ₀ −b) (33) This means that the integration on the time axis has advanced by one step.

Next, in step g3, the absolute value | q _i ¹ of the variation of the variable
-Q _i ⁰ | is determined for i = 1, 2,..., M, and it is checked whether all of them are smaller than a predetermined error limit.

If any absolute value | q _i ¹ −q _i ⁰ | is larger than the error limit, in step h3, the correction value q _i ¹
As new trial values, let them be equal to q _i ⁰ ,
Furthermore the time t ¹ as t ^0, performing integration calculation on further time axis in the same manner as described above.

If the absolute value of the amount of change falls below the error limit at all i, the calculation is terminated assuming that M equivalent circuit constants q ₁ , q ₂ ,..., Q _M have been obtained. The resulting solution is such that at dq / dt = −λ (Aq−b), the time derivative term on the left side of the equation is zero, so the equation Aq =
It is nothing but the solution of b.

Next, regarding the problem of determining the equivalent circuit constant of a transistor,
The results of numerical calculations that have been specifically performed will be described. First, FIG. 4A
If the equivalent-emitter equivalent circuit of the bipolar transistor shown in (1) is regarded as a combination of three four-terminal networks as shown in FIG. 4B, the calculation of the y parameter can be easily performed. If the calculation results of the y-parameter y _i of the intrinsic part 1000 in FIGS. 4A and 4B are shown omitting the progress, the following is obtained.

_{y 11i = {g b [g} e g ee (1-α) -ω 2 c e c ci] + jωg b (c e g ee + c ci g e + c ci g ee)} / Δy i ... (34-1) _{y 12i = - [jωc ci g} b × g e + g ee + jωc e)] / Δy i ... (34-2) y 21i = [g b (αg e g ee + ω 2 c e c ci) -jωc ci g b _{(g e + g ee)]} / Δy i ... (34-3) y 22i = [- ω 2 c e c ci (g e + g ee) + jωc ci (g b g e + g e g ee + g ee g b)] / Δy _i (34-4) Δy _i = g _b g _e + g _b g _ee + g _e g _ee (1-α) −ω ² c _e c _ci + jω [c _e (g _ee + g _b ) + c _ci (g _e + g _ee )] (34-5) Subsequently, the y parameter y ^* of the combination of the intrinsic part 1000 and the external collector capacitance 1001 is given by the following equation according to the composition rule of the four-terminal network.

^{_{y 11 * = y 11i + jωc}} co ... (35-1) y 12 * = y 12i -jωc co ... (35-2) y 21 * = y 21i -jωc co ... (35-3) y 22 * = y 22i + Jωc _co (35-4) Further, the whole including the collector conductance g _c, that is, the y parameter of 1000 & 1001 & 1002 is calculated using the composition rule. If you omit the middle and show only the result, it will be as follows.

_{y 11 = (g c y 11} * + Δy *) / (g c + y 22 *) ... (36-1) y 12 = (g c y 12 *) / (g c + y 22 *) ... (36-2) _{y 21 = (g c y 21} *) / (g c + y 22 *) ... (36-3) y 22 = (g c y 22 *) / (g c + y 22 *) ... (36-4) Δy * ^{_{= y 11 * y 22 * -y}} 12 * y 21 * ... (36-5) FIG. 4A, the circuit constants to be determined in the Figure 4B is, g _e, g
_ee , g _b , g _c , α, c _e , c _ci , and c _co . The objective function that is an indicator of fitting is Given by This represents a relative error in fitting between the theoretical value y _mk of the y parameter and the measured value η _mk . Where m is
1 or 2, where k is the measurement frequency number k = 1, 2,
…, N

Next, for convenience in calculation execution, eight unknowns p _{1 to} p ₈ and respective normalization units r _{1 to} r ₈ are defined as follows.
Note the definition of _{standardized, q i = r i p i} , i = 1,2, note that made ... 8.

_{p 1 = g e (1-} α), r 1 = 1 p 2 = αg e, r 2 = 1 p 3 = g ee, r 3 = 1 p 4 = g b, r 4 = 1 p 5 = g c , r ₅ = 1 p ₆ = c _e , r ₆ = 10 ^-12 p ₇ = c _ci , r ₇ = 10 ^-12 p ₈ = c _co , r ₈ = 10 ^-12 The target measurement data is needed. Originally it should perform calculations using the raw data of the actual device, the demand using well known data properties, this appropriate value (hereinafter to circuit constants p ₁ ~p ₈ True The value is referred to as “value” and the y parameter is used as artificial “measurement data”. The true values of the circuit constants used in the following study are as follows.

g _e = 0.6 siemens p ₁ = g _e (1-α) = 0.012 siemens α = 0.98 = p ₂ = g _e α = 0.588 siemens g _ee = 0.08 siemens = p ₃ g _b = 0.04 siemens = p ₄ g _c = _{0.15siemens = p 5 c e = 200fF} = p 6 c ci = 3fF = p 7 c co = 50fF = p 8 "measurement" frequency gives 40 points for each 0.5GHz starting the 0.1 GHz. As a result, the final value is 19.6 GHz.

To start the calculation, the following numerical values are given as initial trial values of circuit constants.

g _e = 0.7 siemens p ₁ = 0.014 siemens α = 0.98 p ₂ = 0.686 siemens g _ee = 0.1 siemens = p ₃ g _b = 0.2 siemens = p ₄ g _c = 0.2 siemens = p ₅ c _e = 300 fF = p ₆ c _{ci = 5fF = p 7 c co} = 30fF = p true value given to ₈ more and when indicating the four y parameters corresponding to the trial value plotted in the complex plane, FIGS. 6A, second 6D
It looks like the figure. Since the measurement data is artificially generated as described above, it does not include the noisy variation peculiar to the raw data, but as can be seen at a glance, the difference between the true value and the trial value at the same frequency is quite large. Taking the square root of the objective function (37) to quantify the degree of the deviation, y
The magnitude of the relative error vector of the entire parameter is known.

Ie And 6A to 6D show that RE = 3.19 or 31
9%. Desired fitting conditions are roughly 1 error
%, The above initial error is very large.

The time step for integrating on the time axis is 10
^-9 s = 1 ns.

FIG. 7 shows five test runs actually performed by a computer. For the first job “S”, the relative error (RE) = 4.12% as a result of performing the calculation for 500 time steps using λ ₁ to λ ₈ shown in the leftmost column.
This is a significant improvement over the initial error of 319%, but it is still not a good fit. Then, referring to the progress of “S”, executing “W” in which the change in p ₆ (= c _e ) is large and the change in p ₇ is small, RE = 2.
47% improvement. In the next "I", if the final result of "W" is used as an initial condition and proceeds with 500 steps at the same λ, RE is further improved to 1.07%. In the same manner, for “Q”, λ ₇ = 5.E3 (old value = 5.E2) from the same initial condition as “W”.
Running with gives an excellent result of RE = 0.210%. As described above, as a general tendency, a good result with an error of about one digit% is obtained in the first run, and then, starting from the final result, using the same or larger λ as in the previous run, Performing a second run gives better results.

The last “U” starts from the last value of “Q”,
A minimum error of 0.1% is reached in 192 time steps. This is the best result obtained by a series of calculations, and when the obtained y parameters are plotted in FIGS. 6A to 6D, at least on the graph, the results obtained true values for p ₁ to p ₈ It seems to match perfectly. So this fitting is good enough.

Next, FIG. 8 shows a configuration example of the present apparatus for the problem of determining the transistor equivalent circuit constant. The outline of this will be described below. Logic cells corresponding to _eight circuit constants p _{1 to} p 8 are represented by large squares 1 to 8, and generally have a nonlinear input / output response V = g (v). The input terminals of each logic cell have amplifiers A1, A2, ... A8 and resistors, respectively.
R1, R2, ... R3 and capacitors C1, C2, ... C8 are connected, and amplifiers A1, A2,
... A8 has a variable gain corresponding to _{λ 1, λ 2, ... λ} 8. As mentioned above, adjusting and optimizing λ is an important point for the calculation of the present apparatus to be successful.

The capacitors C1, C2,... C8 are connected to a current C (dv
/ dt). The resistance gives a relaxation time constant τ = CR. Current sources S1, S2, ... S8 connected to the input terminals of amplifiers A1, A2, ... A8
Is, as described later, corresponds to the constant term b _i in equation (31-3).

Eight boxes 11, 21, ... connected to the output of each logic cell
81; ... 18, 28,... 88 represent a current source of T _ij V _j , and the quantity q _j obtained by normalizing the j-th parameter p _j with the normalization unit r _j is defined as the output voltage V _j of the j-th cell. Make it correspond. With the above circuit configuration,
The equation that determines the transient response of the network of FIG. 8 can be written as:

_{V i = g (v i)} ... (40) where i is the number of the logic cell, v _i, input of V _i are each i-th cell, representing the output voltage. R _i , C _i , and I _i represent a resistance, a capacitance, and a current source, respectively. T _ij represents the coupling strength between the input of the i-th cell and the output of the j-th cell, and the function g is the input / output response function of each cell. In equations (39) and (40), if the following conditions are given, it can be seen that dq / dt = −λ (Aq−b).

(1) the value obtained by dividing the inter-cell binding strength T _ij in capacitance C _i to correspond with a _ij.

(2) A value obtained by dividing the current I _i by the capacitance C _i is made to correspond to a constant b _i .

(3) The input / output response function g is simply replaced with the linear function V _i = v _i .

(4) The time constant C _i R _i is unnecessary in the equation dq / dt = −λ (Aq−b), so R _i = ∞.

That is, according to the above associations (1) to (4), the calculation of the expression dq / dt = −λ (Aq−b) is automatically performed using the circuit network of FIG.

However, in addition to the above, a current source that actually has voltage dependency
In order to realize T _ij V _j and I _i in an analog manner, a control device as shown at the right end of FIG. 8 is required. This device is
First, measurement data 104 provided by a y-parameter measuring device is required as reference data for fitting. Next, using this measurement data and the output voltage of each logic cell, an arithmetic unit that performs differential operation, multiplication, and addition according to the equations (31-1), (31-2), and (31-3) 103 is required. These calculations may be performed by either a digital or analog method, but a digital method is considered preferable in consideration of calculation accuracy.

Equations (31-1), (31-2),
When the data of (31-3) is a digital quantity, it is converted into an analog signal through the DA converter 101,
Further, using the address unit 100, a predetermined address (i, j)
This is transmitted to both i and j that vary from 1 to 8.

The above is an example of an apparatus for solving the transistor equivalent circuit constant determination problem based on the gist of the present invention.

[Effect of the Invention] As described above, according to the method of the present invention, the matrix equation
When finding the solution of Ax = b, by introducing a sensitivity coefficient λ and replacing it with the integration problem of the differential equation dx / dt = −λ (Ax−b), the inverse matrix calculation x = A ⁻¹ b Without doing
You can reach its purpose. This is a particularly effective solution when the matrix A is large and has many zero elements.

Most of the matrix equations derived from the differential equations that describe the behavior of the actual physical system have such features of large size and many zero elements, so this method is effectively applied to them. Needless to say,

According to the device of the present invention, the calculation of the solution of the matrix equation Ax = b is replaced by the integration problem of the differential equation dx / dt = −λ (Ax−b), and is executed by a device imitated by an electric circuit. As a result, when the circuit reaches a steady state, a solution can be obtained naturally and promptly. The present apparatus is particularly effective when the size of the matrix A is large and the number of zero elements is large (the majority of matrices derived from the differential equations that describe the actual behavior of the physical system are of this kind). give.

[Brief description of the drawings]

FIG. 1 is a flowchart showing a method for analyzing the operation of a semiconductor device according to a first embodiment of the present invention, and FIG.
FIG. 3 is a flow chart showing a method for analyzing the operation of a semiconductor device according to an embodiment, FIG. 3 is a circuit diagram showing a device for analyzing the operation of a semiconductor device according to a first embodiment of the present invention, and FIG. FIG. 4B is an equivalent circuit diagram of FIG. 4A, FIG. 5 is a flowchart showing a method for analyzing a problem of determining an equivalent circuit constant of a semiconductor device according to the third embodiment of the present invention, and FIGS. FIG. 6D is a graph in which four y parameters corresponding to a true value and a trial value are plotted on a complex plane, and FIG. 7 is a method for analyzing a problem of determining an equivalent circuit constant of a semiconductor device according to a third embodiment of the present invention. FIG. 8 is a circuit diagram showing an apparatus for analyzing a problem of determining an equivalent circuit constant of a semiconductor device according to a second embodiment of the present invention, and FIG. 9 is a diagram showing mesh division of a semiconductor device space. Fig. 10 is electronic It is a diagram showing a hole transport equation and matrix-vector equation of the simultaneous equations consisting of Poisson's equation. 1, 2, 12: Logic cell, 11, 12, 1212: Current detecting means, 100: Address unit, 101: D / A converter, 103, 104: Circuit.

──────────────────────────────────────────────────続 き Continued on the front page (58) Field surveyed (Int.Cl. ⁷ , DB name) H01L 21/66

Claims (6)

(57) [Claims] 1. A method for analyzing a system of simultaneous equations comprising electron and hole transport equations and Poisson equation in a semiconductor device modeling by a computer, wherein the system of simultaneous equations is converted into a time differential term dp / dt, dn / dt, dψ /
(a), including dt and sensitivity coefficients λ _p , λ _n , λ _{下 記}
(B), substituting the expression _{(c), dp / dt = -λ p f} p ... (a) dn / dt = λ n f n ... (b) dφ / dt = λ φ f ψ ... (c) semiconductor devices arrangement of the dividing points (M, n) and determines the sensitivity coefficient lambda _p, lambda _n, to lambda _[psi giving proper spatial position dependency which matches the physical characteristics of the semiconductor device, the (a) , (B), and (c) are replaced by the following (d), (e),
(F), dp (M, N) / dt = −λ _p (M, N) × f _p (M, N) (d) dn (M, N) / dt = λ _n (M , N) × f _n (M, N) (e) dψ (M, N) / dt = λ _ψ (M, N) × f _ψ (M, N) (f) The above (d), (e) A method for analyzing the operation of a semiconductor device, wherein a solution of the system of simultaneous equations is obtained by time-integrating the equations (f) and (f). 2. A method for analyzing a Poisson equation representing a relationship between a potential and a space charge in a semiconductor device modeling by a computer, wherein the Poisson equation includes a time differential term ∂ψ / ∂t and a sensitivity coefficient λ. substituting the _{equation, ∂ψ / ∂t = λ ψ f} ψ ... (a) dividing point arrangement of semiconductor devices (M, N) and determines the even sensitivity coefficient lambda proper which matches the physical characteristics of the semiconductor device giving a spatial position dependency, the converted formula (a) to (b) the following equation, ∂ψ (M, N) / ∂t = λ (M, N) × f ψ (M, N) (B) An operation analysis method of a semiconductor device, wherein a solution of the Poisson equation is obtained by time-integrating the above equation (b). 3. A method of analyzing the equivalent circuit constant determination problem of a semiconductor device based on the least square method by a computer, wherein the equivalent circuit constant determination problem includes a time differential term dq / dt and a sensitivity coefficient λ. Dq / dt = −λ (Aq−b) (a) The sensitivity coefficient λ is given a value suitable for the physical characteristics of each equivalent circuit constant, and the above equation (a) is integrated over time. A method for analyzing the operation of a semiconductor device, wherein the solution of the problem of determining an equivalent circuit constant is obtained. In 4. A semiconductor device modeling, electronic, following differential equation dp / dt = -λ and a hole transport equations and system of simultaneous equations consisting of the Poisson equation time differential terms and sensitivity coefficients _p f _p dn / dt _{= λ n f n dψ / dt} = λ ψ, the device for analyzing replaced with f _[psi, the device includes a single or plurality of m logic cells, the bound adjustable amplification gain for each logic cell lambda m amplifiers having _i or multipliers equivalent thereto, and detection means provided in a network for each of the logic cells and detecting that a transient phenomenon in the device has entered a steady state. An analyzer characterized by the above-mentioned. In 5. The semiconductor device modeling, an apparatus for differential equation ∂ψ / ∂t = replaced with lambda _[psi f _[psi analysis including a time differential terms and sensitivity coefficients a Poisson equation representing the relationship between the potential and space charge, wherein The apparatus comprises one or more m logic cells, m amplifiers or equivalent multipliers coupled to each logic cell and having an adjustable amplification gain λ _i , and An analyzing apparatus, comprising: detecting means provided in a network form for detecting that a transient phenomenon in the apparatus has entered a steady state. 6. A problem of determining an equivalent circuit constant of a semiconductor device based on a least square method is represented by a differential equation dq / dt = −λ × (Aq−b)
The apparatus comprises: m or more m logic cells and m amplifiers coupled to each of the logic cells and having an adjustable amplification gain λ _i or a multiplier equivalent thereto And a current detecting means provided in a network for each of the logic cells and detecting that a current state in the device has reached a steady state.