JP2964992B2 - Diffusion simulation method - Google Patents

Description

DETAILED DESCRIPTION OF THE INVENTION

[0001]

The present invention relates to a diffusion simulation method, and more particularly to a diffusion simulation method for performing a numerical simulation of impurity diffusion inside a device in a computer simulation of a semiconductor device manufacturing process.

[0002]

2. Description of the Related Art Computer simulation of thermal diffusion of impurities, which is one of the manufacturing processes of semiconductor devices, is analyzed as described in "Process Device Simulation Technology", edited by Danyo, pp. 26-28. A method is adopted in which a power region is divided into meshes, the diffusion equation is discretized at each mesh point, and further linearized by Newton's method or the like and converted into a system of linear equations to obtain a solution of the diffusion equation.

Recently, as a diffusion model that more faithfully reflects the physical phenomenon of impurity diffusion, for example, IEEE
Trans. Electron Devices, Vol. 40, No. 7, pp. 1215 to 1222; Han
e, H .; Matsumoto, "A Model f
or Boron Short Time Annea
ling After Ion Implantati
A point defect diffusion model has been proposed as described in Equations (6)-(11) on page 1216 of "on".

An equation describing such a point defect diffusion model is, for example, as described in equation (9) on page 1216 of the same document, the diffusion species concentration change contribution due to diffusion and drift:

[0005] The contribution of the diffusion species concentration change due to the interaction / reaction between each diffusion species:

[0006]

The configuration is as follows. From the viewpoint of the connection between the equations, it can be considered that f1 represents the connection between adjacent mesh points of the same diffusion type, and f2 represents the connection between each diffusion type at the same mesh point.

In general, in an impurity diffusion simulation in a manufacturing process of a semiconductor device, it is necessary to collectively solve the mutually connected diffusion equations.
The simultaneous linear equation obtained by linearizing the diffusion equation by the ewton method has a dimension of mesh number × (diffusion number + 1). +1 in parentheses is a contribution of an equation for obtaining an electrostatic potential な る serving as a driving force for drifting a diffusion species.

Usually, as shown in FIG. 3A, the coefficient matrix of this simultaneous linear equation first forms a mesh point block by grouping equations and arranging them in the order of mesh point numbers. Further, a method of arranging each equation in numerical order, or a method of arranging mesh point groups together and arranging them in equation number order to form an equation block as shown in FIG. It consists of either.

When this simultaneous linear equation is solved by using a Krylov subspace method such as a conjugate gradient method with preprocessing, FIG.
In the configuration method, the block incomplete LU decomposition method in which each mesh point block is virtually regarded as one matrix element can be used as an effective pre-processing method. Tend to be used. In the method of constructing the coefficient matrix as shown in FIG. 3A, the increase in the total amount of calculation is suppressed to approximately O (n) with respect to the number n of mesh points, but the amount of calculation for one mesh point block is reduced to the number m of equations. O for
(M ² ), the number of equations increases, the total amount of calculation becomes enormous, and it becomes difficult to find a solution in a practical time.

As one method for efficiently obtaining a solution by approximating a large-scale diffusion equation group as described above, for example, SIS
PAD'96, pages 93-94, M.P. Kawakam
i, M. Sugaya, S .; Kamohara, "A
New High-Speed Non-equil
Ibrium Point Defect Model
for Annealing Simulation ”
There is a method described in. In this conventional example, a coefficient matrix of a diffusion equation relating to diffusion type 1 is configured in the form of FIG.

As described above, the diffusion term f1 and the interaction / reaction term f
When the diffusion equation is discretized in the time domain and integrated, a method using an implicit method for f1 and an explicit method for f2 is proposed. More specifically, the previous time t _k-1, the current time and t _k, a backward Euler method implicit method, when the explicit method assuming the case of employing a forward Euler method, seeds 1 diffuser at the current time The concentration C ₁ (t _k ) is given by the equation:

Is obtained by solving At the current time, the concentration of each diffusion species at the previous time, C ₁ (t
_k-1 ), C ₂ (t _k-1 ),... C _m (t _k-1 ) are known, so the value of f2 is known, and the components resulting from f2 are included in the coefficient matrix of the diffusion equation. Not 1 / (t _k −t
_k-1 ) and f1 only contribute to the coefficient matrix. This means that no coupling component with other diffusion species is included in the diffusion equation coefficient matrix of the diffusion species 1 and the small matrix forming each mesh point block in FIG. 3A is a diagonal matrix. . Therefore, the calculation amount of one mesh point block is the number of equations m
However, only the increase in O (m) can be obtained, and a solution can be obtained in a practical time.

[0014]

In the conventional diffusion simulation method described above, when the diffusion equation is discretized in the time domain and integrated, an explicit method is used for f2. There is an upper limit on the time step width that can be integrated into the equation. The third document, page 94, FIG.
FIG. 4 shows the oscillation of the total dose of the boron-interstitial silicon pair caused by exceeding the upper limit of the time step width. The authors of this example show that even if the dose of the boron-interstitial silicon pair oscillates over time, the final spatial distribution of boron they actually need is greater than if there is no oscillation. This vibration is accepted as there is no difference. However, such vibration is obviously non-physical, and this method is not applicable to those who intend to evaluate the temporal change of the boron-interstitial silicon pair concentration by simulation.

An object of the present invention is to perform a numerical simulation of a phenomenon in which a plurality of diffusion species diffuse while interacting / reacting with each other in a practical time, and an arbitrary length when time-integrating a diffusion equation. It is intended to provide a diffusion simulation method capable of obtaining a solution of a diffusion equation such that a stable solution without vibration / divergence is always obtained even when the time step width is used.

[0016]

According to a first diffusion simulation method of the present invention, a diffusion equation including an interaction term / reaction term between diffusion species is numerically solved by a computer, and the inside of the device is manufactured in a semiconductor device manufacturing process. In the first step, only the interaction / reaction term portion in the diffusion equation is solved by an implicit method to obtain the movable diffusion species concentration at the current time. Then, in the second step, the movable diffusion species concentration obtained in the first step is virtually regarded as the diffusion species concentration at the previous time, and only the diffusion term in the diffusion equation is solved by an implicit method. It is configured to obtain the final diffusion species concentration at the current time.

A second diffusion simulation method according to the present invention is a computer-implemented diffusion simulation method for solving the definition of the concentration of active diffusion species and the diffusion equation numerically by a computer to obtain a change in the impurity distribution inside the device in the semiconductor device manufacturing process. First, in the first step, the active diffusion species concentration definition equation is solved to determine the movable diffusion species concentration at the current time, and then, in the second step, the movable diffusion species determined in the first step is determined. The active diffusion species concentration at the current time is obtained by solving the diffusion equation implicitly assuming that the concentration is virtually the diffusion species concentration at the previous time, and the activity between the previous time and the current time is further obtained in the third step. The configuration is such that the final diffusion species concentration at the current time is obtained from the diffusion species concentration change.

[0018]

According to the present invention, the interaction / reaction term part and the diffusion term part in the diffusion equation are divided into two stages to determine the movable diffusion species concentration at the current time and the final diffusion species concentration at the current time. Seeking independent. Since there is no spatial differential operator in the diffusion term in the step of time-integrating the interaction / reaction term part, the coefficient matrix constructed in the form of FIG. It becomes an angular matrix. The total amount of computation Therefore this step interaction when the average mesh points adjacent to each mesh point and n _Adj /
Calculation amount O (n) when the reaction and diffusion terms are solved at once
It is about 1 / n _Adj of × O (m ² ). In the step of partial time integration of the diffusion term, there is no interaction / reaction term, so there is no coupling between diffusion species, and the small matrix of each mesh point block of the coefficient matrix configured in the form of FIG. It becomes a matrix. Therefore, the total operation amount of this step is the operation amount O when the interaction / reaction term and the diffusion term are collectively solved.
(N) × O (m ² ) is reduced to On × O (m). Therefore, in the present invention, the amount of operation O (n) × O (m ² ) when solving the diffusion equation collectively is O (n) ×
(O (m ² ) / n _Adj + O (m)), and a solution can be obtained in a practical time. On the other hand, first in the present invention,
Since the implicit solution method is used for the time integration of the differential equation in both the second step, a stable solution without vibration / divergence can be obtained for an arbitrary time step size.

[0019]

Embodiments of the present invention will now be described with reference to the drawings.

FIG. 1 is a flowchart of a processing procedure for solving a boron-point defect pair diffusion equation for explaining the first embodiment of the present invention.

Step 101 is a step of generating a mesh used for spatially discretizing a diffusion equation inside a given semiconductor device structure. Step 102 is a step of setting an initial value of a diffusion species concentration on each mesh point. And 1
03 is a step of setting the initial value of the analysis time to 0. 104 is a step of setting a new time step width and updating the analysis time, and 105 is a step of taking out only the interaction / reaction terms from the original diffusion equation and performing time integration by an implicit method to obtain a movable diffusion species concentration. 106 is to extract only the inner diffusion term of the original diffusion equation,
In this step, the final diffused species concentration at the current time is obtained by integrating the movable diffused species concentration determined at step 05 with the diffused species concentration at the previous time virtually and performing time integration by the implicit method. 107
Is a step of checking whether or not the analysis time has reached a predetermined end time, and if the analysis time has arrived, terminates the processing; otherwise, returns the processing to step 104.

First, the diffusion equation of a boron-interstitial silicon pair is as follows.

Here, (6), (7), (8), (9)
Equations are the diffusion equations for vacancies, neutral interstitial silicon, + 1-valent interstitial silicon, and boron-interstitial silicon pairs, respectively. Equations (10) and (11) are the interactions of boron and boron clusters, respectively. It is a reaction equation. No. (1
Equation (2) is a charge neutral condition equation, and is a relational equation between electron / hole concentration and electrostatic potential:

This is used to determine a local electrostatic potential by cooperating with.

Next, this embodiment will be described with reference to the flowchart shown in FIG.

First, in step 101, a mesh is formed in an internal region of a semiconductor device for simulating impurity diffusion. Next, in step 102, an initial density value of each diffusion type is set on each mesh point generated in step 101. At step 103, the analysis time is initialized, and at step 104, the analysis time is advanced by one. In step 105, only the interaction / reaction term of the diffusion equation is taken out and time-integrated by an implicit method to obtain a movable diffusion species concentration C _i ^mov . When the receding Euler method is applied as an implicit solution to the boron / point defect pair diffusion model taken up in the present embodiment, the following equation is specifically solved.

Here, regarding the boron (C _B− ) and the boron cluster (C _BBcls ), since the equations (19) and (20) do not originally include a diffusion term, the concentration of the movable diffusion species in this step is And the final diffused species concentration is determined. Although not a differential equation, equation (21) is solved simultaneously with other differential equations in this step. Next, at step 106, the movable diffusion species concentration obtained at step 105 is virtually regarded as the diffusion species concentration at the previous time, and only the inner diffusion term of the diffusion equation is taken out and time-integrated by using the implicit method to obtain the final value. Request diffusion species concentration C _i. When the receding Euler method is applied as an implicit method to the above-described boron / point defect pair diffusion model, the following equation is solved.

Here, the boron (C _B− ), the boron cluster (C _BBcls ), and the charge neutrality condition have no diffusion term, so this step is skipped. Finally, in step 107, it is checked whether or not the analysis end time specified in advance has been reached. If the time has reached, the simulation ends, and if not, the flow returns to step 104.

In this embodiment, the concentration of the movable diffusion species is first determined, and then the concentration is diffused as a virtual initial state. Therefore, the finally determined diffusion species concentration does not strictly satisfy the original diffusion equation, and an approximation error occurs. However, this error has a time step size of O (δt _k ), and is almost the same as the truncation error of the backward Euler method used most often in the implicit time integration method. Therefore, the accuracy of the solution does not deteriorate significantly compared to the conventional example,
Rather, a more accurate solution can be obtained because the stability is improved. Further, as described in the section of the operation, in the present invention, the amount of calculation when handling the coefficient matrix can be reduced as compared with the case where the original diffusion equation is solved collectively, and the simulation can be performed in a practical time.

Next, an example in which the diffusion equations of boron, phosphorus, and arsenic described by the following Fermi model in the present embodiment are simultaneously solved will be described.

Here, equations (26), (28) and (30) are the diffusion equations for boron, phosphorus and arsenic, respectively. The diffusion term and drift term in each equation are the active impurity concentration C _Ba ,
Described by C _Pa and C _Asa . The second (2
Equations (7), (29), and (31) are the equations for defining the active impurity concentration of boron, phosphorus, and arsenic, respectively. The active impurity concentration of boron and phosphorus is limited by their respective capacities C _Bss and C _Pss . Is a model formula in which the rate of increase in active impurity concentration is suppressed by clustering.

Although these active impurity concentration definition formulas are originally to be described by an ordinary differential equation relating to time, when the response time constant of the movable impurity concentration is sufficiently shorter than the diffusion time, the system is brought into an equilibrium state. It can be regarded as a certain one, and thus can be replaced with a stationary active impurity concentration definition formula. In this case, there is no need to integrate the equations over time, and the active impurity concentration definition equation has no relation to the stability of the system. The last expression (32) is a charge neutral condition expression.

FIG. 2 is a flow chart for explaining a second embodiment of the present invention.

Reference numeral 201 denotes a step of generating a mesh used for spatially discretizing a diffusion equation inside a given semiconductor device structure; and 202, a step of setting an initial value of a diffusion species concentration on each mesh point. And 2
03 is a step of setting the initial value of the analysis time to “0”. 204 is a step of setting a new time interval and updating the analysis time. 205 is a step of solving the active diffusion species concentration definition formula to obtain the movable diffusion species concentration.
Step 6 is a step of obtaining the final active diffusion species concentration at the current time by integrating the diffusion equation with the implicit method by considering the movable diffusion species concentration obtained in step 205 virtually as the diffusion species concentration at the previous time. 207 is a step for obtaining a final diffusion species concentration from a change in the active diffusion species concentration between the previous time and the current time, and 208 is a step for checking whether or not the analysis time has reached a predetermined end time. If so, the process is terminated, and if not, the process returns to step 204.

First, in step 201 of FIG. 2, a mesh is formed in an internal region of a semiconductor device for simulating impurity diffusion. Next, in step 202, step 201
The initial density value of each diffusion species is set on each mesh point generated in the step. At step 203, the analysis time is initialized, and at step 204, the analysis time is advanced by one. In step 205, the active impurity concentration C _i ^mov is determined by solving the active diffusion species concentration definition equation. For the Fermi diffusion model taken up in the present embodiment, the following equation is specifically solved.

Here, although it is not an active concentration definition equation, equation (36) is also solved in this step together with the active concentration definition equation. Next, at step 206, the movable diffusion species concentration obtained at step 205 is virtually regarded as the diffusion species concentration at the previous time, and the diffusion equation is time-integrated using an implicit solution method to thereby obtain the final active diffusion species concentration C _ia. Ask for. When the backward Euler method is applied as an implicit solution to the Fermi model, the following equation is solved.

In step 207, the final diffusion species concentration at the current time is obtained from the change in the active diffusion species concentration between the previous time and the current time as follows.

Finally, in step 208, it is checked whether or not the analysis end time specified in advance has been reached. If it has reached, the simulation is ended, and if not, the flow returns to step 204.

[0039]

As described above, according to the present invention, the part of the interaction / reaction term and the part of the diffusion term in the diffusion equation are divided into two steps to determine the concentration of the movable diffusion species at the current time and the current time. Since the final diffusion species concentration is obtained independently, it is possible to obtain a solution in a practical time compared to the case where the diffusion equation is solved collectively. Since the implicit solution method is used, there is an effect that a stable solution without vibration / divergence can be obtained for an arbitrary time interval.

[Brief description of the drawings]

FIG. 1 is a flowchart of a processing procedure for explaining a first embodiment of the present invention.

FIG. 2 is a flowchart of a processing procedure for explaining a second embodiment of the present invention.

FIG. 3 is a diagram showing a method of forming a coefficient matrix of a simultaneous linear equation obtained by discretizing a diffusion equation in a diffusion simulation method.

[Explanation of symbols]

 101-107, 201-208 steps

Claims (2)

(57) [Claims] 1. A diffusion simulation method using a computer for numerically solving a diffusion equation including an interaction term / reaction term between diffusion species by a computer to obtain a change in impurity distribution inside the device in a semiconductor device manufacturing process, First, in the first step, only the interaction / reaction term portion in the diffusion equation is solved by the implicit method to obtain the concentration of the movable diffusion species at the current time, and then, in the second step, the first step The final diffused species concentration at the current time was obtained by solving the obtained movable diffused species concentration virtually as the diffused species concentration at the previous time and solving only the diffusion term in the diffusion equation by the implicit method. A diffusion simulation method, characterized in that: 2. A diffusion simulation method by a computer for numerically solving an active diffusion species concentration definition equation and a diffusion equation by a computer to obtain a change in impurity distribution inside a device in a semiconductor device manufacturing process. Solve the active diffusion species concentration definition formula to determine the movable diffusion species concentration at the current time, and then in the second step,
Obtain the active diffusion species concentration at the current time by solving the diffusion equation implicitly assuming the movable diffusion species concentration determined in the first step as the diffusion species concentration at the previous time virtually,
Further, in a third step, a final diffusion species concentration at the current time is obtained from a change in the active diffusion species concentration between the previous time and the current time.