JP2009224677A - Ion implantation distribution generating method, and simulator - Google Patents

Abstract

<P>PROBLEM TO BE SOLVED: To provide an ion implantation distribution generating method and a simulator, wherein a computation time is shortened by effectively using approximation. <P>SOLUTION: In the ion implantation generating method, a projection (R<SB>p</SB>(E))<SP>(2)</SP>taking into consideration up to a secondary term of a projection R<SB>p</SB>of a range of an ion implanted in a semiconductor is obtained by using a secondary perturbation model using an approximation expression of (R<SB>p</SB>(E))<SP>(2)</SP>=(R<SB>p</SB>(E))<SP>(1)</SP>+Δ<SB>p</SB><SP>(2)</SP>(E) by using a perturbation term Δ<SB>p</SB>(2) when a known projection taking up to a primary term into consideration is denoted as (R<SB>p</SB>(E))<SP>(1)</SP>. <P>COPYRIGHT: (C)2010,JPO&INPIT

Description

  The present invention relates to an ion implantation distribution generation method and a simulator, and relates to a method for generating an ion implantation distribution in a shorter time with substantially the same accuracy as a Monte Carlo simulator.

In a silicon integrated circuit device, introduction of impurities into a silicon substrate is generally performed by ion implantation.
In the process construction of such a silicon integrated circuit device, it is necessary to determine ion implantation conditions for obtaining a necessary element structure, and such ion implantation conditions are determined by simulation.
There is Monte Carlo as a means for theoretically predicting the ion implantation distribution into the amorphous layer. In this method, the trajectory of incident ions is traced based on the physics of nuclear stopping power and electron stopping power for the interaction between the incident ions and the substrate.

This theory is also effective in a general case where ions are implanted into an arbitrary substrate, and the accuracy can be further improved by tuning the electron stopping power.
FIG. 15 shows a calculation model. An ion having a mass number M ₁ , an atomic number Z ₁ , and an energy T _1i (velocity v _1i ) interacts with an atom having a mass number M ₂ and an atomic number Z ₂ constituting the substrate. The energy T _2f to be transmitted is given by _assuming that the scattering angle is Φ, the velocity of ions after interaction is ν _1i , and the velocity of substrate atoms is ν _2i .
T _2f / T _1i = 2M ₂ ν _2i ² sin ² (Φ / 2) / [(1/2) M ₁ ν _1i ² ]
= (4M ₂ / M ₁ ) {[M ₁ v _1i / (M ₁ + M ₂ )] ² / ν _1i ² }
× sin ² (Φ / 2)
= [4M ₂ M ₁ / (M ₁ + M ₂ ) ² ] sin ² (Φ / 2) (1)
It is expressed as

つまり、注入されたイオンは、核との相互作用により、
ΔＥ_n ＝Ｔ_2f
のエネルギーを失う。 Here, assuming that the distance between the ions and the substrate atoms is r, the collision parameter is b, and the potential energy is V (r), the transfer energy T _2f is obtained as the following equation (2).

In other words, the implanted ions are interacted with the nucleus,
ΔE _n = T _2f
Lose energy.

ここで、半径ｂの円周上の位置、即ち角度θは、Ｒａｎｄ（ｎ）をｎが０から１の間の乱数とすると、
θ＝２πＲａｎｄ（ｎ）
の関係から求まる。 Further, the scattering angle Φ accompanying the interaction is expressed by the following formula (3).

Here, the position on the circumference of the radius b, that is, the angle θ is Rand (n), where n is a random number between 0 and 1,
θ = 2πRand (n)
It is obtained from the relationship.

This determines the energy and direction after the collision.
Next, the same calculation is repeated with the implantation energy newly to trace the ion trajectory.
However, each time an actual collision occurs, the integration cost represented by the above equation (3) is executed every time. Therefore, Ziegler proposes Magic formula and solves this integration with a simple protocol. It is.
In the simulator TSUPRM, the calculation results are tabulated and this value is evaluated by interpolation.

但し、ａ_B をボーア径とすると、
ρ＝ｒ／ａ_U ，
η＝ｂ／ａ_U ，
ａ_U ＝０．８８５４ａ_B ／（Ｚ₁ ^0.23 ＋Ｚ₂ ^0.23 ）
である。
また、エネルギーは、
ε＝Ｅ_C ／（Ｚ₁ Ｚ₂ ｅ² ／ａ_U ）
Ｅ_c ＝（１／２）Ｍ_c ｖ_1i ²
１／Ｍ_c ＝１／Ｍ₁ ＋１／Ｍ₂
である。
また、ポテンシャルエネルギーＶ（ｒ）として、下記のＺｉｅｇｌｅｒ−Ｌｉｔｍａｒｋ−Ｂｉｅｒｓａｋ（ＺＬＢ）のポテンシャルエネルギーを用いる（例えば、非特許文献１参照）。
Ｖ（ｒ）＝（ｅ² Ｚ₁ Ｚ₂ ／ｒ）ｆ（ρ）
但し、

In addition, when each parameter in said Formula (3) is represented by the universal variable which normalized, it represents by the following formula (4).

However, if a _B is Bohr diameter,
ρ = r / a _U ,
η = b / a _U ,
a _U = 0.8854a _B / (Z ₁ ^0.23 + Z ₂ ^0.23 )
It is.
Energy is also
ε = E _C / (Z ₁ Z ₂ e ² / a _U )
E _c = (1/2) M _c v _1i ²
_{1 / M c = 1 / M} 1 + 1 / M 2
It is.
Further, the potential energy of the following Ziegler-Litmark-Biersak (ZLB) is used as the potential energy V (r) (see, for example, Non-Patent Document 1).
V (r) = (e ² Z ₁ Z ₂ / r) f (ρ)
However,

  However, even if the above-mentioned calculation cost reduction method is introduced, the Monte Carlo simulation follows each particle trajectory, so it is necessary to perform tens of thousands of calculations to reduce statistical errors, which takes time. There is a problem.

Therefore, when the electronic stopping power S _e, nuclear stopping power S _n given, the integral equation to follow the ion implantation distribution, Lindhart, Scharf, proposed by Schiott (e.g., see Non-Patent Document 2), this model It is called LSS theory.

  In this theory, if the injection conditions are determined without tracking the particle trajectory, the energy dependence of the distribution moment of any order of the distribution can be calculated immediately by solving the integral equation.

In this case, it is necessary to expand the integral equation, reduce it to a differential equation, and solve it, and an approximation is entered at that step.
In other words, assuming that the transmitted energy is smaller than the initial energy, the moment is expanded in terms of energy, and the order of expansion is the rank of the differential equation to be solved.

The scattering angle also develops a Taylor expansion with respect to energy, which is related to the coefficient of the equation to be solved and is independent of the rank. Analytical models are generally obtained with arbitrary coefficients only for first-order linear differential equations.
In this LSS theory, Taylor expansion is performed up to the first order with respect to the moment and second order with respect to the scattering angle, and a linear first-order differential equation is solved to derive analytical model expressions for R _p and ΔR _p .

Here, with reference to FIG. 16, the integral equation regarding the range R of the implanted ions in the LSS theory will be described.
FIG. 16 is a schematic diagram of the ion range R. The ions implanted with energy E move as shown in the figure while losing energy while interacting with the substrate atoms, changing the direction, and losing all energy to the substrate. Quiesce inside.

At this time, the line integral of the trajectory along which the ions have traveled is defined as a range R, the projection of the range R in the injection direction is R _p , and the spread in the lateral direction is R _T.
Further, the x component of the lateral spread is Δx, and the y component is Δy.
In the drawings and the subsequent integral equations or differential equations, the lateral spread is represented by a symbol with a suffix “vertical symbol” added to R. However, in the specification text, “ R _T ".

である。
また、Ｒに関してｍ次のモーメントを

と定義する。 Let P (E, R) be the probability that ions implanted with energy E perpendicular to the surface will stop between R and R + dR. in this case,

It is.
In addition, m-order moment with respect to R

It is defined as

The probability that an ion interacts with atoms and electrons as it travels dR from the surface is
NdR∫dσ _n + NdR∫dσ _e (8)
And the energy of the ions in their respective interactions is
E- _Tn , E- _Te
Assuming that
Here, dσ _n and dσ _e are differential cross sections related to nuclear stopping power and electron stopping power.

で与えられる。
また、ｄＲ進む間に衝突しない確率は、
１−（ＮｄＲ∫ｄσ_n ＋ＮｄＲ∫ｄσ_e ） ・・・（１０）
である。この間にエネルギーは失われないから、イオンがＲに止まる確率は、
［１−（ＮｄＲ∫ｄσ_n ＋ＮｄＲ∫ｄσ_e ）］Ｐ（Ｅ，Ｒ−ｄＲ） ・・・（１１）
となる。よって、

となる。 Therefore, the probability that this ion interacts with atoms and electrons, loses energy, and stops at the range R is

Given in.
Also, the probability of not colliding while traveling dR is
1- (NdR∫dσ _n + NdR∫dσ _e ) (10)
It is. Since no energy is lost during this time, the probability that the ion stops at R is
[1- (NdR∫dσ _n + NdR∫dσ _e )] P (E, R-dR) (11)
It becomes. Therefore,

It becomes.

となる。両辺にＲ^m をかけてＲについて０から∞まで積分すると、

となる。 Here, considering the limit of taking dR to infinity and small,

It becomes. If you integrate R ^m on both sides and integrate R from 0 to ∞,

It becomes.

式（１５）における左辺の第１項は０である。また

と表現するが、これはＲ^m の期待値に相当する。 When the left side of Equation (14) is partially integrated and the order of integration on the right side is changed,

The first term on the left side in equation (15) is zero. Also

This corresponds to the expected value of R ^m .

となり、これが飛程Ｒに関するｍ次の積分方程式である。
ここで、飛程Ｒに関しては１次に関してのみ解析する。式（１７）においてｍ＝１とおいて

を得る。これが、飛程Ｒの従うべき積分方程式である。 Here, using equation (16), equation (15) becomes

This is the m-th order integral equation for the range R.
Here, the range R is analyzed only for the primary. In equation (17), m = 1

Get. This is the integral equation that the range R should follow.

Next, Equation 17 is solved under the approximation of T _n , T _e << _E , leaving up to the first order term. This approximation assumes that the transmitted energy is small.
In this case, corresponds to the approximation ignoring angle scattering for interaction with nuclei, for interaction with the electrons, the more energy is high T _e increases. For this reason, the accuracy of approximation is improved when the energy is relatively low.
This solution is expressed as <R (E)> ^(1), considering that it is a first order solution.

となる。ここで、
Ｓ_n ＝∫Ｔ_n ｄσ_n ，Ｓ_e ＝∫Ｔ_e ｄσ_e ・・・（２０）
を利用している。式（１９）より、良く知られている、

が導出される。 Next, by dividing both sides of equation (18) by N and expanding the Taylor,

It becomes. here,
_{S n = ∫T n dσ n,} S e = ∫T e dσ e ··· (20)
Is used. Well known from Equation (19),

Is derived.

Next, with reference to FIG. 17, the geometric relationship between R _p (E, cos φ), R _T (E, cos φ), R _p (E), R _T (E) will be described.
FIG. 17 is an explanatory diagram of the geometric relationship of R _p (E, cos φ), R _T (E, cos φ), R _p (E), R _T (E), and is incident at energy E and angle φ. The situation where ions are stationary at point A in the substrate is schematically shown.
Assuming that R _T (E, cos φ) is a lateral spread on a plane perpendicular to the axis of incidence angle 0, the projection R _p (E) with respect to the axis perpendicular to the incidence direction is transverse to the plane perpendicular to the incidence direction. The direction spread can be considered as R _T (E). Therefore,

と定義する。 Next, based on the above geometrical relationship, an integral equation to be followed by R _p , R _T , R _p R _T is derived.
Let P _p (E, R _p , cos φ) be the probability that the projected range of ions scattered at an angle φ with energy E stops between R _p and R _p + dR _p .
Let P _p (E, R _p ) be the probability when φ = 0. In addition, the m-th moment with respect to R _p is

It is defined as

The probability that the incident ion interacts with atoms and electrons while traveling by dR _p is the same as in the range R,
NdR _p ∫dσ _n + NdR _p ∫dσ _e (25)
And the energy in this interaction is
E- _Tn , E- _Te
To decrease.

で与えられる。
また、ｄＲ進む間に衝突しない確率は、
１−（ＮｄＲ_p ∫ｄσ_n ＋ＮｄＲ_p ∫ｄσ_e ） ・・・（２７）
である。この間にエネルギーは失われないから、イオンがＲ_p に止まる確率は、
［１−（ＮｄＲ_p ∫ｄσ_n ＋ＮｄＲ_p ∫ｄσ_e ）］Ｐ（Ｅ，Ｒ_p −ｄＲ_p ）
・・・（２８）
となる。よって、

となる。 In this case, regarded as zero scattering angle by interaction with electrons, the ions interact with atoms and electrons, the probability of stopping in R _p loses energy T _n, a T _e is

Given in.
Also, the probability of not colliding while traveling dR is
1- (NdR _p ∫dσ _n + NdR _p ∫dσ _e ) (27)
It is. Since no energy is lost during this time, the probability that the ion stops at R _p is
_{[1- (NdR p ∫dσ n +} NdR p ∫dσ e)] P (E, R p -dR p)
... (28)
It becomes. Therefore,

It becomes.

となり、両辺にＲ^m _p をかけてＲ_p について０から∞まで積分すると、

とｍ次の〈Ｒ^m _p （Ｅ）〉に関する積分方程式が導出される。 Here, considering the limit of making dR _p infinitesimal,

Next, the sides in over R ^m _p integrating from 0 for R _p to ∞,

And an integral equation for m-th order <R ^m _p (E)> is derived.

である。よって、このイオンが相互作用しエネルギーＴ_n ，Ｔ_e を失ってＲ_T に止まる確率は、

で与えられる。 Next, let us consider R _T representing the spread of the distribution in the horizontal direction.
Let P _T (E, R _T , cos φ) be the probability that ions implanted at an energy E and an incident angle φ stop between R _T and R _T + dR _T.
In this case, the probability that an ion incident angle φ interact with atoms and electrons between proceeds dR _T is

It is. Therefore, the probability of stopping in the R _T The ions lose energy T _n, T _e interact,

Given in.

である。この間にエネルギーは失われず、角度も変化しないためＲ_T も変化しない。よって、イオンがＲ_T に止まる確率は、

となる。 Also in this case, it is assumed that there is no change in angle in the interaction with electrons, and in this case, _RT is not changed. In addition, the probability that does not conflict between advance dR _T,

It is. During this time, no energy is lost, and the angle does not change, so _RT does not change. Therefore, the probability that the ion stops at _RT is

It becomes.

となる。 From the above,

It becomes.

となり、両辺にＲ^m _T をかけてＲ_T について０から∞まで積分すると、

と〈Ｒ^m _T （Ｅ）〉に関する積分方程式が導出される。 Here, considering the limit of making dR _T infinitesimal,

Next, the sides in over R ^m _T be integrated from 0 for R _T to ∞,

And the integral equation for <R ^m _T (E)> is derived.

ここで、１次のテーラー展開した解を〈Ｒ_p （Ｅ）〉⁽¹⁾ とおくと、下記の式（４０）となる。

First, in order to obtain the projection R _p of the projected range, it becomes When m = 1 in formula (31) and equation (39) below.

Here, if the solution of the first-order Taylor expansion is set as <R _p (E)> ⁽¹⁾ , the following equation (40) is obtained.

但し、α（Ｅ）及びα（Ｅ′）は、下記の式（４２）で表される。

From this, the following equation (41) is obtained.

However, α (E) and α (E ′) are represented by the following formula (42).

Next, in order to obtain the second moments ΔR _p and ΔR _pt , if m = 2 in the equations (31) and (38), the following equations (43) and (44) are obtained.

Here, Lindhard introduced the variables <R ² _c (E)> and <R ² _r (E)> defined by the following formulas (45) and (46).

<ΔR ² _p (E)> and <R ² _T (E)> are the variables <R ² _c (E)> and <R ² _r (E)> and the following equations (47) and (48). Associated.

Further, the lateral deviation <ΔR ² _pt (E)> is related to <ΔX ² >, not <R ² _T (E)> itself. Therefore,
<R ² _T (E)> = <ΔX ² > + <ΔX ² > = 2 <ΔR ² _pt (E)> (49)
And from this,
<ΔR ² _pt (E)> = <R ² _T (E)> / 2 (50)
It becomes.

Lindhart succeeded in obtaining independent differential equations of the following equations (51) and (52) by introducing these variables.

Gibbons also substitutes the energy-independent nuclear stopping power into this model formula, evaluates R _p and ΔR _p of B, N, Al, P, Ga, As, In, and Sb in Si. Compared with the experimental data, the agreement is good.
Furthermore, numerical evaluation of moments up to the third order has been studied.
J. et al. F. Ziegler, J. et al. P. Biersack, and U.S. Littmark, The stopping and range of ions in Solid, Pergamon, 1885 J. et al. Lindhart, M .; Scharff, H.M. Schott, Mat. Fts. Med. Vid. Sclsk, vol. 33, pp. 1-39, 1963

However, in the LSS theory, approximation is necessary when solving the integral equation. In the models proposed so far, there is no problem with R _p, but ΔR _p has a problem that accuracy is significantly reduced or calculation is not possible.

  Further, as described above, it has been studied to numerically evaluate moments up to the third order, but these numerically analyzed data are not systematically compared. In comparison, there is also a problem that the currently used Pearson distribution cannot be used only with moments up to the third order.

  That is, the results derived from the LSS theory vary depending on the types of nuclear stopping power model, electron stopping power model, and how the integral equation is developed. The accuracy of the approximation handled in the LSS theory can be accurately evaluated by comparing with the Monte Carlo calculation using the same nuclear stopping power model and electron stopping power model.

However, no comparison of Monte Carlo and LSS theory has been made so far with high nuclear stopping power and electron stopping power. Furthermore, the actual distribution can be expressed with high accuracy only by using the fourth moment.
Therefore, in the conventional LSS theory, it cannot be associated with the currently used Pearson function.

Also, when obtaining the ion implantation distribution from the correlation between the collision parameter b and the scattering angle Φ, since Monte Carlo follows each trajectory of the particle, it is necessary to perform tens of thousands or more calculations in order to reduce the statistical error. take time.
The integration that evaluates the scattering angle takes the longest time to proceed with this calculation.

Also, how much the distance L between the collision and the next collision affects the calculation time.
Until now, the average interatomic distance L _min has been used. In this case, there is a problem that the calculation time becomes enormous because the scattering is calculated even when the interaction is negligible.

Therefore, the present invention is based on the integral equation of the LSS theory using the same nuclear stopping power and electron stopping power as the Monte Carlo simulator that reproduces the experimental data, and R _p , ΔR _p , ΔR _pt by second-order and third-order perturbation calculations. The purpose is to derive the analytical model of γ and significantly reduce the calculation time.

  Also, when using Monte Carlo, an object is to significantly reduce the calculation time by replacing the integral for evaluating the scattering angle with an approximate expression that expresses analytically.

According to one aspect of the present invention, the projection <R _p (E)> ⁽²⁾ that takes into account the quadratic term of the range R _p of the range of ions implanted into the semiconductor is known to the first term. When the projection is <R _p (E)> ⁽¹⁾ , using the perturbation term Δ _p ⁽²⁾ (E),
< _Rp (E)> ⁽²⁾ = < _Rp (E)> ⁽¹⁾ + _Δp ⁽²⁾ (E)
An ion implantation generation method obtained using a second-order perturbation model using the approximate expression is provided.

  Further, according to another aspect, when the scattering angle Φ of ions implanted into the substrate is determined by the Monte Carlo method, the scattering angle Φ is determined by the constant high-angle term, the low-angle term, and the intermediate-angle term. An ion implantation distribution generation method is provided that is expressed as an analytical expression by the reciprocal of the sum of three terms.

  Further, according to another aspect, a simulator is provided in which calculation formulas relating to the various ion implantation distribution generation methods described above are incorporated.

According to the present invention, a highly accurate analysis model can be derived, and higher-order moments can be predicted with high accuracy. This makes it possible to analytically determine the parameters of the Pearson distribution that is currently used as a standard.
Therefore, the energy dependence of the distribution of ion implantation of an arbitrary impurity into an arbitrary amorphous substrate can be predicted instantaneously.
The distribution generation software can also predict the depth dependency of the horizontal distribution.

  Furthermore, even when using Monte Carlo, the calculation time can be greatly shortened by replacing the integral for evaluating the scattering angle with an approximate expression expressing analytically.

また、ユニバーサルな核阻止能Ｓ_n （ε）は、下記の式（５４）で与えられる。

First, the primary to tertiary stopping powers are calculated.
The primary nuclear stopping power S _n (E) is given by the following equation (53).

The universal nuclear stopping power S _n (ε) is given by the following equation (54).

Here, when the above ZLB potential is used, the universal nuclear stopping power S _n (ε) is approximated by the following equation (55).

但し、ユニバーサルな拡張された核阻止能Ω_n ² （ε）は下記の式（５７）で表される。

これも、下記の式（５８）で近似的に表現される。

This stopping power is extended to higher order, and the expanded nuclear stopping power Ω _n ² corresponding to the second moment is expressed by the following equation (56).

However, the universal expanded nuclear stopping power Ω _n ² (ε) is expressed by the following equation (57).

This is also expressed approximately by the following equation (58).

但し、ユニバーサルな拡張された核阻止能Λ_n ³ （ε）は下記の式（６０）で表される。

これも、下記の式（６１）で近似的に表現される。

Further, the expanded nuclear stopping power Λ _n ³ corresponding to the third-order moment is expressed by the following equation (59).

However, the universal extended nuclear stopping power Λ _n ³ (ε) is expressed by the following equation (60).

This is also expressed approximately by the following formula (61).

Here, based on the above formulas (22) and (23), the parameters R _p ² (E, cos φ), R _p ³ (E, cos φ), R _p R _T ² (E , Cosφ) is calculated in advance as follows.

Here, as shown in FIG. 17 above, a point on a circle having a radius R _T (E) centered on C on the plane where point C exists can be considered equivalent to A because of the symmetry of space. it can.
Therefore, it can be assumed that β takes a value from 0 to 2π with the same probability.
Therefore, the value obtained by integrating the equations (62) to (64) from 0 to 2π with respect to β is divided by 2π and the average value is evaluated.

を参考にしてこれらの平均値は、

となる。ここで、

である。

These average values with reference to

It becomes. here,

It is.

が得られる。 Next, an integral equation related to the cross term necessary for obtaining a higher-order moment such as γ was derived as the following equation (73).

Is obtained.

となる。
但し、Ω_n ² ＝∫Ｔ_n ² ｄσ_n ，Ω_e ² ＝∫Ｔ_e ² ｄσ_e
を利用している。
なお、一般に、Ω_e ² ≪Ω_n ² であるので、実際の計算ではΩ_e ² は無視する。 Next, for the range projection R _p , the solution obtained by expanding the above equation (39) to the second order is expressed as <R _p (E)> ⁽²⁾ .

It becomes.
However, Ω _n ² = ∫T _n ² dσ _n , Ω _e ² = ∫T _e ² dσ _e
Is used.
In general, since Ω _e ² << Ω _n ² , Ω _e ² is ignored in the actual calculation.

The above equation (74) is a second-order differential equation and cannot be solved analytically. Therefore, a second-order approximation formula <R _p (E)> ⁽²⁾ is obtained by perturbation approximation, using the known <R _p (E)> ⁽¹⁾ and perturbation term Δ _p ⁽²⁾ that can be calculated.
<R _p (E)> ⁽²⁾ = <R _p (E)> ⁽¹⁾ + Δ _p ⁽²⁾ (E) (75)
It expresses.

が得られ、これより、摂動項Δ_p ⁽²⁾ は、下記の式（７７）で表される。

但し、式（７７）におけるζ_p ⁽²⁾ （Ｅ）は下記の式（７８）で表される。

Substituting this equation (75) into the above equation (74) and dropping the product of Ω ² and Δ _p ^{(2) as} a third-order minute term, the equation (40) relating to the first-order expansion is replaced with By using

Thus, the perturbation term Δ _p ⁽²⁾ is expressed by the following equation (77).

However, ζ _p ⁽²⁾ (E) in the equation (77) is expressed by the following equation (78).

Thus, by approximating the unknown <R _p (E)> ⁽²⁾ with a known computable <R _p (E)> ⁽¹⁾ and the perturbation term Δ _p ⁽²⁾ , A value <R _p (E)> ⁽²⁾ based on a second-order perturbation model of the projection R _p can be numerically calculated.

となる。
但し、Λ_n ³ ＝∫Ｔ_n ³ ｄσ_n ，Λ_e ³ ＝∫Ｔ_e ³ ｄσ_e
を利用している。 Similarly, for the range projection R _p , the solution obtained by expanding the above equation (39) to the third order is denoted as <R _p (E)> ⁽³⁾ .

It becomes.
However, Λ _n ³ = ∫T _n ³ dσ _n , Λ _e ³ = ∫T _e ³ dσ _e
Is used.

Next, since <R _p (E)> ⁽³⁾ cannot be calculated as it is, <R _p (E)> ⁽²⁾ and the perturbation term Δ _p ^(3), which are already known from the above equation (77 ^), are changed. make use of,
<R _p (E)> ⁽³⁾ = <R _p (E)> ⁽²⁾ + Δ _p ⁽³⁾ (E) (80)
It expresses.

但し、式（８１）におけるζ_p ⁽³⁾ （Ｅ）は下記の式（８２）で表される。

When calculated in the same manner as in the second-order case using this approximate expression, the perturbation term Δ _p ⁽³⁾ is expressed by the following expression (81).

However, ζ _p ⁽³⁾ (E) in the equation (81) is expressed by the following equation (82).

In this case as well, the range of unknown <R _p (E)> ⁽³⁾ is approximated by <R _p (E)> ⁽²⁾ and perturbation term Δ _p ⁽³⁾ , which are known and can be calculated. It is possible to numerically calculate the value <R _p (E)> ^{(3) according} to the third-order perturbation model of the projection R _p of.

ここで、未知の〈Ｒ_c ² （Ｅ）〉⁽²⁾ を、計算可能な既知の〈Ｒ_c ² （Ｅ）〉⁽¹⁾ と摂動項Δ_c ²⁽²⁾を用いて、
〈Ｒ_c ² （Ｅ）〉⁽²⁾ ＝〈Ｒ_c ² （Ｅ）〉⁽¹⁾ ＋Δ_c ²⁽²⁾ ・・・（８４）
とおいて、上記の式（８３）を整理すると、左辺は、
２〈Ｒ_p （Ｅ）〉⁽²⁾ ／Ｎ＝２｛〈Ｒ_p （Ｅ）〉⁽¹⁾ ＋Δ_p ⁽²⁾ （Ｅ）｝／Ｎ であるので、下記の式（８５）となる。

Next, ΔR _p and ΔR _pt are examined.
First, when a solution obtained by expanding the integral equation of Equation (51) to the second order is expressed as <R _c ² (E)> ⁽²⁾ , the following Equation (83) is obtained.

Here, an unknown <R _c ² (E)> ⁽²⁾ is converted into a computable known <R _c ² (E)> ⁽¹⁾ and a perturbation term Δ _c ^{2 (2)} .
<R _c ² (E)> ⁽²⁾ = <R _c ² (E)> ⁽¹⁾ + Δ _c ^{2 (2)} (84)
When the above equation (83) is rearranged, the left side is
Since 2 <R _p (E)> ⁽²⁾ / N = 2 {<R _p (E)> ⁽¹⁾ + _Δp ⁽²⁾ (E)} / N, the following equation (85) is obtained.

但し、式（８６）におけるζ_c ⁽²⁾ （Ｅ）は下記の式（８７）で表される。

Therefore, the perturbation term Δ _c ^{2 (2)} is expressed by the following equation (86).

However, ζ _c ⁽²⁾ (E) in the equation (86) is represented by the following equation (87).

ここで、未知の〈Ｒ_c ² （Ｅ）〉⁽³⁾ を、上記の結果により既知となっった〈Ｒ_c ² （Ｅ）〉⁽²⁾ と摂動項Δ_c ²⁽³⁾を用いて、
〈Ｒ_c ² （Ｅ）〉⁽³⁾ ＝〈Ｒ_c ² （Ｅ）〉⁽²⁾ ＋Δ_c ²⁽³⁾ ・・・（８９）
とおいて、上記の式（８８）を整理すると、左辺は、
２〈Ｒ_p （Ｅ）〉⁽³⁾ ／Ｎ＝２｛〈Ｒ_p （Ｅ）〉⁽²⁾ ＋Δ_p ⁽³⁾ （Ｅ）｝／Ｎ であるので、下記の式（９０）となる。

Similarly, when a solution obtained by expanding the integral equation of equation (51) to the third order is represented as <R _c ² (E)> ⁽³⁾ , the following equation (88) is obtained.

Here, the unknown <R _c ² (E)> ⁽³⁾ is converted into the known <R _c ² (E)> ⁽²⁾ and the perturbation term Δ _c ^{2 (3).} ,
<R _c ² (E)> ⁽³⁾ = <R _c ² (E)> ⁽²⁾ + Δ _c ^{2 (3)} (89)
When the above equation (88) is rearranged, the left side is
Since 2 <R _p (E)> ⁽³⁾ / N = 2 {<R _p (E)> ⁽²⁾ + Δ _p ⁽³⁾ (E)} / N, the following equation (90) is obtained.

但し、式（９１）におけるζ_c ⁽³⁾ （Ｅ）は下記の式（９２）で表される。

Therefore, the perturbation term Δ _c ^{2 (3)} is expressed by the following equation (91).

However, ζ _c ⁽³⁾ (E) in the equation (91) is represented by the following equation (92).

ここでも、未知の〈Ｒ_r ² （Ｅ）〉⁽²⁾ を、計算可能な既知の〈Ｒ_r ² （Ｅ）〉⁽¹⁾ と摂動項Δ_r ²⁽²⁾を用いて、
〈Ｒ_r ² （Ｅ）〉⁽²⁾ ＝〈Ｒ_r ² （Ｅ）〉⁽¹⁾ ＋Δ_r ²⁽²⁾ ・・・（９４）
とおいて、上記の式（９３）を整理すると、左辺は、
２〈Ｒ_p （Ｅ）〉⁽²⁾ ／Ｎ＝２｛〈Ｒ_p （Ｅ）〉⁽¹⁾ ＋Δ_p ⁽²⁾ （Ｅ）｝／Ｎ であるので、下記の式（９５）となる。

Next, when a solution obtained by expanding the integral equation of the equation (52) to the second order is _expressed as <R _r ² (E)> ⁽²⁾ , the following equation (93) is obtained.

Here again, the unknown <R _r ² (E)> ⁽²⁾ is transformed into the computable known <R _r ² (E)> ⁽¹⁾ and the perturbation term Δ _r ^{2 (2)} .
^{_{<R r 2 (E)>}} (2) = <R r 2 (E)> (1) + Δ r 2 (2) ··· (94)
When the above equation (93) is rearranged, the left side is
Since 2 <R _p (E)> ⁽²⁾ / N = 2 {<R _p (E)> ⁽¹⁾ + Δ _p ⁽²⁾ (E)} / N, the following equation (95) is obtained.

但し、式（９６）におけるζ_r ⁽²⁾ （Ｅ）は下記の式（９７）で表される。

Therefore, the perturbation term Δ _r ^{2 (2)} is expressed by the following equation (96).

However, ζ _r ⁽²⁾ (E) in the equation (96) is expressed by the following equation (97).

ここで、未知の〈Ｒ_r ² （Ｅ）〉⁽³⁾ を、上記の結果により既知となっった〈Ｒ_r ² （Ｅ）〉⁽²⁾ と摂動項Δ_r ²⁽³⁾を用いて、
〈Ｒ_r ² （Ｅ）〉⁽³⁾ ＝〈Ｒ_r ² （Ｅ）〉⁽²⁾ ＋Δ_r ²⁽³⁾ ・・・（９９）
とおいて、上記の式（９８）を整理すると、下記の式（１００）となる。

Similarly, when a solution obtained by expanding the integral equation of Equation (52) to the third order is represented as <R _r ² (E)> ⁽³⁾ , the following Equation (98) is obtained.

Here, the unknown <R _r ² (E)> ⁽³⁾ is converted into the known <R _r ² (E)> ⁽²⁾ and the perturbation term Δ _r ^{2 (3).} ,
^{_{<R r 2 (E)>}} (3) = <R r 2 (E)> (2) + Δ r 2 (3) ··· (99)
If the above equation (98) is rearranged, the following equation (100) is obtained.

但し、式（１０１）におけるζ_r ⁽³⁾ （Ｅ）は下記の式（１０２）で表される。

Therefore, the perturbation term Δ _r ^{2 (3)} is expressed by the following equation (101).

However, ζ _r ⁽³⁾ (E) in the equation (101) is expressed by the following equation (102).

Using the above parameters, ΔRp ⁽²⁾ and ΔR _pt ^{(2) based} on the ^second- order perturbation model are used by Lindhardt in the first-order perturbation model.
<ΔR _p ² (E)> = <R _p ² (E)> + <R _p (E)> ²
= {<R _c ² (E)> + <R _r ² (E)>} / 3− <R _p (E)> ² (103)
<R _T ² (E)> = 2 {<R _c ² (E)> − <R _r ² (E)>} / 3
... (104)
<ΔR _pt ² (E)> = <R _T ² (E)> / 2 (105)
Evaluate with.
However, in the second-order case, different from the first-order, the first equation of the above-described equation (103) is expressed by considering that it is a second-order perturbation model.
[<ΔR _p ² (E)> ⁽¹⁾ + Δ _p ^{2 (2)} ]
= [<R _p ² (E)> ⁽¹⁾ + Δ _p ^{2 (2)} (E)]
− [<R _p (E)> ⁽¹⁾ + Δ _p ⁽²⁾ (E)] ² ... (106)
It becomes.
The same applies to equations (104) and (105).

Next, as for ΔRp ⁽³⁾ and ΔR _pt ^{(3) based} on the ^third- order perturbation model, as Lindhard used in the first-order perturbation model,
<ΔR _p ² (E)> = <R _p ² (E)> + <R _p (E)> ²
= {<R _c ² (E)> + <R _r ² (E)>} / 3− <R _p (E)> ² (107)
<R _T ² (E)> = 2 {<R _c ² (E)> − <R _r ² (E)>} / 3
... (108)
<ΔR _pt ² (E)> = <R _T ² (E)> / 2 (109)
Evaluate with.
However, in the case of the third order, the first formula of the above formula (107) is expressed by considering that it is a third order perturbation model unlike the first order.
[<ΔR _p ² (E)> ⁽²⁾ + Δ _p ^{2 (3)} ]
= [<R _p ² (E)> ⁽²⁾ + Δ _p ^{2 (3)} (E)]
− [<R _p (E)> ⁽²⁾ + Δ _p ⁽³⁾ (E)] ² ... (110)
It becomes.
The same applies to equations (108) and (109).

Next, in order to obtain γ which is a third-order moment, if m = 3 in equation (31), the following equation (111) is obtained.

Further, when m = 1 and k = 2 are substituted in the above-mentioned cross term formula (73), the following formula (112) is obtained.

この変数の導入により、下記の式（１１５）及び式（１１６）のように、それぞれ変数が〈Ｒ_c ³ （Ｅ）〉及び〈Ｒ_r ³ （Ｅ）〉のみの独立な微分方程式を得ることができる。

Here, it tries to arrange | equalize Formula (111) and Formula (112) by introducing the variable represented by the following formula | equation (113) and Formula (114).

By introducing this variable, independent differential equations having only variables <R _c ³ (E)> and <R _r ³ (E)> are obtained as in the following formulas (115) and (116). Can do.

The variables <R _c ³ (E)> and <R _r ³ (E)> in this case are related to the target moment by the following equations (117) and (118).

On the other hand, the third-order moment μ ₃ is defined as a concentration distribution obtained by standardizing N (x).
μ ₃ = ∫ (x− <R _p (E)>) ³ N (x) dx (119)
It is represented by
Therefore, the third moment μ ₃ is
μ ₃ = ∫ (x− <R _p (E)>) ³ N (x) dx
= ∫ (x ³ −3x ² <R _p (E)> + 9x <R _p (E)> ²
− <R _p (E)> ³ ) N (x) dx
= <R _p ³ (E)>-3 <R _p ² (E)><R _p (E)> + 2 <R _p (E)> ³ (120)
It becomes.

Therefore,
γ = {∫ (x−R _p (E)) ³ N (x) dx} / ΔR _p ³ (121)
The third moment γ representing the asymmetry of the distribution defined by is evaluated by the following equation (122).

摂動項Δ_c ³⁽ⁱ⁾についての一般解が、式（１２５）として得られる。

但し、式（１２５）におけるζ_c ³⁽ⁱ⁾は、それぞれ、下記の式（１２６）乃至式（１２８）で表される。

When the solution of the above differential equation is put into the following formula (123) and formula (124),

A general solution for the perturbation term Δ _c ^{3 (i} ) is obtained as equation (125).

However, ζ _c ^{3 (i} ) in the equation (125) is represented by the following equations (126) to (128), respectively.

但し、式（１２９）におけるζ_r ³⁽ⁱ⁾は、それぞれ、下記の式（１３０）乃至式（１３２）で表される。

For the other variable <R _r ³ (E)>, a general solution for the perturbation term Δ _r ^{3 (i} ) is obtained as equation (129).

However, ζ _r ^{3 (i} ) in the equation (129) is represented by the following equations (130) to (132), respectively.

Next, the depth dependence of the lateral spread ΔR _pt will be examined.
Spread [Delta] R _pt lateral are said to have a depth-dependent, the dependency, upon the [Delta] R _pt0 and [Delta] R _pt in x = R _p,
ΔR _pt (x) = ΔR _pt0 + m (x−R _p ) (133)
It is expressed by

Therefore, an attempt is made to evaluate m from <R _pT2 ³ (E)>.
First, R _T ² and ΔR _pt are as described above,
R _T ² = Δx ² + Δy ²
= 2ΔR _pt ²
= 2 [ΔR _pt0 + m (x−R _p )] ^{2
_{≒ 2ΔR pt0 2 + 4mΔR pt0 (}} x-R p) ··· (134)
Associated with.
However, in the last approximation, since x is the concentration at a point far from R _p decreases suddenly, because m the vicinity R _p is important, terms related to m is assumed to be very small.

したがって、ｍは、

で表現することができる。 On the other hand, <R _pT2 ³ (E)> has a relationship of m and the following equation (135).

Therefore, m is

Can be expressed as

  On the other hand, when the ion implantation distribution is generated by the Monte Carlo simulator, the integration of the above equation (3) is not performed after each collision, but the relational expression between the collision parameter b and the scattering angle Φ is numerically solved. Then, an analytical approximate expression is derived by fitting the result using a graph.

In this approximation, it may be expressed by a constant high angle term (π) and a low angle term, or by a constant high angle term (π), a low angle term, and a medium angle term. good.
In this case, the low angle term and the medium angle term are represented by the product of the term including the normalized collision parameter η and the normalized energy ε.

  Further, after adding the constant high angle term (π), the low angle term, and the medium angle term to powers of numbers greater than one, the scattering angle Φ is approximated by taking the power of the reciprocal thereof.

Based on the above, the ion implantation distribution generation method according to the first embodiment of the present invention will now be described with reference to FIGS.
FIG. 1 is a flowchart of an ion implantation distribution generation method according to Embodiment 1 of the present invention. An implantation condition including a substrate type, an implanted impurity species, implantation energy, and a dose is input. Next, as explained in Equation 75 to Equation 106 above,
b. Using the second-order perturbation model, the projection R _{p of the} range, which is the first-order moment, the deviation ΔR _p of the projection R _p , which is the second-order moment, and the deviation ΔR _{pt in the} lateral direction are obtained. Then
c. Together to generate a Gaussian distribution of 1-dimensional distribution of R _p and [Delta] R _p obtained, R _p, generates a Gaussian distribution of a two-dimensional from [Delta] R _p and [Delta] R _pt.

FIG. 2 compares R _p and ΔR _p obtained by the ion implantation distribution generation method of Example 1 of the present invention with R _p and ΔR _p obtained by other simulations.
Here, along with the result by the Monte Carlo method, the result by the first-order perturbation model by the LSS theory and the result by the above-mentioned Gibbons are shown as the 1.5th-order perturbation model.

As is apparent from the figure, all the simulation results are in good agreement with respect to the range projection R _p which is the first moment.
However, with respect to the deviation ΔR _p of the projection R _p that is a second-order moment, the result of the second-order perturbation model of Example 1 of the present invention is the Monte Carlo method compared to the result of the first-order perturbation model according to the LSS theory. There is a very good agreement with the results of.
Note that in the 1.5th order perturbation model, the deviation ΔR _p of the projection R _p cannot be calculated because it includes a divergence term in the denominator of the integral equation.

FIG. 3 shows a one-dimensional distribution obtained by the ion implantation distribution generation method according to the first embodiment of the present invention. In this example, P (phosphorus) ions are implanted into a Si substrate at 10 keV at 1 × 10 ¹⁵ cm ^−2. Results are shown.
As is clear from the figure, a beautiful Gaussian distribution is shown.

FIG. 4 is a two-dimensional distribution obtained by the ion implantation distribution generation method according to the first embodiment of the present invention. Here, again, P (phosphorus) ions are implanted into a Si substrate at 10 keV at 1 × 10 ¹⁵ cm ^−2. Results are shown.
In this case, the simulation is performed assuming that the ion implantation mask does not transmit ions at all.

As described above, in the first embodiment of the present invention, the projection R _p , the deviation ΔR _p , and the lateral deviation are obtained in a very short time by using a second-order perturbation model as compared with the Monte Carlo method. ΔR _pt can be obtained with an accuracy close to the Monte Carlo method.

In the first embodiment, since the second-order perturbation model is used, only the second-order moment can be obtained with high accuracy, and the ion implantation distribution that can be generated is a Gaussian distribution.
However, the Pearson distribution that is actually used in the LSS process simulation as a standard requires higher-order moment information such as γ and β.

Here, a method for generating an ion implantation distribution according to the second embodiment of the present invention using a third-order perturbation model in order to obtain third-order and fourth-order moments and generate a Pearson distribution closer to the actual measurement result will be described. .
FIG. 5 is a flowchart of an ion implantation distribution generation method according to the second embodiment of the present invention. An implantation condition including a substrate type, an implanted impurity species, implantation energy, and a dose is input. Next, as explained in Equation 79 to Equation 135 above,
b. Using a third-order perturbation model, the range projection R _p , which is the first moment, the deviation ΔR _p of the projection R _p , which is the second order moment, and the lateral deviation ΔR _pt have depth dependence. Asking. The third-order moment γ is obtained, and the fourth-order moment β is obtained empirically. Then
c. Obtained R _p, ΔR _p, together with generating the gamma and β Karakara 1-dimensional _{distribution, R p, Δ R p,} γ, and generates a two-dimensional distribution of β and depth dependent [Delta] R _pt.

Here, a method for obtaining β, which is a fourth-order moment, will be described.
β can be obtained in principle by solving the LSS integral equation in the same manner as other moments, but becomes more complicated than the case of the third-order moment γ and is numerically unstable.
Therefore, since the ion implantation distribution does not change so peculiarly, it is indicated by Monte Carlo that universal is related between β and γ, and an empirical expression is used.

FIG. 6 shows the simulation result by Monte Carlo as the relationship between β and γ.
Here, the results are shown when B, As, and P are ion-implanted into the Si substrate at 10 keV, 20 keV, 40 keV, 80 keV, 160 keV, and 320 keV, respectively.

This correlation is
β = 2.661 + 1.852γ ² (137)
It is expressed as
Using this approximate expression, β may be obtained from γ obtained by a third-order perturbation model to generate the above-described one-dimensional distribution and two-dimensional distribution.

However, this relational expression (137) does not exist in the commonly used Pearson IV region as shown in the figure, but crosses the Pearson I, II, III, and VI regions.
Therefore, as shown in the figure as β _D2 , β is increased by 5% and forcibly moved to the critical value β _{D2 represented} by the following formula (137) of Pearson IV.
β _D2 = [3 (13γ ² +16) +6 (γ ² +4) ^3/2 ] / (32−γ ² )
... (138)

Next, a method for generating an ion implantation distribution according to the third embodiment of the present invention that analytically expresses the relational expression between the collision parameter b and the scattering angle Φ to reduce the calculation cost will be described. When T _2f becomes small to some extent, for example, when it becomes about 10 eV, it is not necessary to consider a larger b.
In other words, it can be considered that there is no interaction with the nucleus.

With such a transmission energy T _2fc, so can be considered to be the scattering angle [Phi _c in this case small, sin ² (Φ / ²⁾ can be approximated by Φ _c ^2/4.
Therefore, the above equation (1) is
T _2fc / T _1i = _{[4M 2 M 1 / (M 1} + M 2) 2 ] × Φ _c ^2/4. . (139) 8 From this, the scattering angle Φ _c is determined as Φ _c = {[(M ₁ + M ₂ ) ² / M ₂ M ₁ ] × (T _2fc / T _1i )} ^1/2 (140).
Here, when the collision parameters corresponding to the scattering angle [Phi _c and b _max, scattering angle [Phi _c is expressed by the following equation (141).

但し、ここでは、 η＝ｂ_max ／ａ_U ・・・（１４３）
であり、ｂ_max に対応するρがρ_min である。 Actually, since it is computationally expensive to obtain b _max from this integration for each collision, the relational expression between b and Φ is numerically solved in advance from equation (3), and the result is expressed analytically by fitting. Use the formula.
In the fitting, the following equation (142) represented by universal variables similar to the above equation (3) is used.

Here, η = b _max / a _U (143)
Ρ corresponding to b _max is ρ _min .

FIG. 7 is an explanatory diagram of the energy ε dependence of the scattering angle Φ, and the result of numerically solving the above equation (142) is shown by a solid line.
Here, for each η, the result of fitting,
Φ _c = [(1 / π) + 1 / (a (η) ε ^−0.98 )] ⁻¹ (144)
The result of combining with is shown by a broken line.

As shown in the figure, the fitting is not successful in the region where the angle increases and approaches saturation. However, there is no problem because the angle actually used is an angle of several degrees, that is, an area of 10 ⁻² radian. In addition, this area is simple,
Φ _c = a (η) ε ^-0.98 (145)
It is expressed as However, a ((eta)) is represented by the following formula | equation (146).

これから、
ａ＝Φ_c ／ε^-0.98 ・・・（１４８）
の場合のηが、上記の式（１４６）から求まる。
さらに、ηが求まると、上記の式（１４２）からｂ_max が求まる。 However, since η (a) in which η is represented by a is convenient for practical use, fitting is performed from the relational expression of a−η.
FIG. 8 is an explanatory diagram of the η dependence of a in the above formula (146), and FIG. 9 shows this as the a dependence of η in reverse.
Η (a) obtained by fitting from FIG. 9 is expressed by the following equation (147).

from now on,
a = Φ _c / ε ^−0.98 (148)
In this case, η is obtained from the above equation (146).
Furthermore, when η is obtained, b _max is obtained from the above equation (142).

The collision cross section σ _n of the interaction with the nucleus is σ _n = πb _max ² (149)
Can be considered.
Where the concentration of substrate atoms is N:
LNσ _n = 1 (150)
It can be considered that a collision occurs once the distance L is reached.

Therefore, a phenomenon is described in which an ion of energy E travels without interaction with the nucleus during the distance L and then interacts with the nucleus. It is considered that it only interacts with electrons while traveling L.
That is, while traveling L, interaction with electrons causes ΔE _c = S _e (E) L (151)
Lose energy.
The interaction with the nucleus after proceeding L selects the collision parameter b at a point on the circle from 0 to b _max at random. That is, b = b _max [Rand (n)] ^1/2 (152)
And
However, Rand (n) is a random number between 0 and 1.

When L is smaller than the average distance L _min (= 1 / N ^1/3 ) between atoms,
L = L _min (153)
And shall be. In this case, the area related to the collision is a rectangle, and σ _n = L _min ² (154)
It becomes.
When this is re-baked into the collision parameters that have been treated so far, b _max in this case is obtained from the above equation (149).
b _max = L _min / π ^1/2 (155)
It becomes.
Therefore, it is only necessary to generate a random number according to the above equation (152) using this b _max to determine b, and the injection energy becomes a predetermined minimum critical energy at which the above equation (141) is established with high accuracy. Calculate until it reaches.

FIG. 10 shows each L of the trajectory of 10 particles. Here, L is normalized by L _min , and L gradually decreases from about 10 times L _min .
Therefore, the number of calculations is remarkably reduced compared with the case where L _min is always used, and the calculation time is shortened.
In the figure, L of 10 particles are indicated by L1 to L10, respectively.

FIG. 11 is an explanatory diagram of an ion implantation distribution, and shows an ion implantation distribution when B ions are implanted at 1 × 10 ¹⁵ cm ⁻² at an acceleration energy of 40 keV.
As is clear from the figure, the ion implantation distribution obtained with the variable L in Example 3 of the present invention is almost the same as the calculation result when L _min is used, and the measurement result by SIMS is well expressed. ing.

Here, handling in the case where the substrate is composed of a plurality of atoms will be considered.
Suppose there are n _i atoms i. The critical angle Φ _ci of each atom is obtained as Φ _ci = {[(M ₁ + M _2i ) ² / M ₁ M _2i ] × (T _2fc / T _1i )} ^1/2 ... (156). The area σ _ni can be calculated, and the total collision cross-sectional area σ _n is expressed by the following equation (156).

Here, when the concentration of the substrate atoms is N,
LNσ _n = 1 (158)
It can be considered that a collision occurs once the distance L is reached.

そこで、その確率に対応させて乱数を発生させて原子ｉを選択する。選択した後は構成原子が一種類の場合と同じである。
また、Ｌが、下記の式（１６０）より小さくなった場合には、Ｌ＝Ｌ_min としなければならない。

The probability of colliding with atom i by the collision is expressed by the following equation (159).

Therefore, an atom i is selected by generating a random number corresponding to the probability. After selection, it is the same as in the case of one kind of constituent atom.
When L becomes smaller than the following formula (160), L = L _min must be set.

In this case, since the collision cross section is determined by the atomic density, if the probability of collision with the atom i is the same as in the case of two types, b _maxi is
b _maxi = (rσ _ni / π) ^1/2 (161)
It becomes. However, r is represented by the following formula (162).

  As described above, in Example 3 of the present invention, when the scattering angle is small, the approximate expression obtained by fitting is used instead of integration, and the ion implantation distribution is analytically calculated using the variable L. Therefore, the number of calculations can be remarkably reduced without reducing the calculation accuracy.

Next, in order to improve the accuracy in the transition region so as to be effective even when the scattering angle is large, an ion implantation distribution generation method according to the fourth embodiment of the present invention in which the equation (144) is expanded will be described.
FIG. 12 is an explanatory diagram of the energy ε dependency of the scattering angle Φ, and the result of numerically solving the above equation (3) is shown by a solid line.
Moreover, FIG. 13 is the figure which expanded the range of 1-3 radians in FIG.
Here, for each η, the result of fitting,
Φ = 1 / {(1 / π) ^3/2 + [1 / (a (η) ε ^−0.98 )] ^3/2
+ [1 / g (η) ε ^−0.2 ] ^3/2 ] ^2/3 (163)
The result of combining with is shown by a broken line.

Here, the term g (η) corresponding to the intermediate scattering angle is
g (η) = [1 / (4.0η ^−0.46 ) + 1 / (9.0η ^−1.3 )] ^{−1
= 4.0η -0.46 /(1+0.44η 0.84) ··· (} 164)
Η dependence shown in FIG. 14 is seen.

As is clear from FIG. 12 and FIG. 13, the result of numerical calculation by integration is reproduced with high accuracy from a low angle to a high angle.
Note that (1 / π) ^3/2 in the above equation (163) corresponds to a high angle term.

It is a flowchart of the ion implantation distribution generation method of Example 1 of this invention. Is an illustration of R _p and [Delta] R _p determined by ion implantation distribution generating method of Example 1 of the present invention. It is the one-dimensional distribution calculated | required by the ion implantation distribution generation method of Example 1 of this invention. It is the two-dimensional distribution calculated | required by the ion implantation distribution generation method of Example 1 of this invention. It is a flowchart of the ion implantation distribution generation method of Example 2 of this invention. It is the figure which showed the simulation result by Monte Carlo as the relationship between (beta) and (gamma). It is explanatory drawing of energy (epsilon) dependence of scattering angle (PHI). It is explanatory drawing about (eta) dependence of a in Formula (146). It is explanatory drawing about a dependence of (eta). It is explanatory drawing of L of each locus | trajectory. It is explanatory drawing of ion implantation distribution. It is explanatory drawing of energy (epsilon) dependence of scattering angle (PHI). It is an enlarged view of the range of 1-3 radians in FIG. It is explanatory drawing of (eta) dependence of g. It is a calculation model. It is a schematic diagram of the ion range R. _{R p (E, cosφ),} R T (E, cosφ), an illustration of the geometric relationships _{R p (E), R T} (E).

Claims (6)

Projection of considering quadratic term projection R _p of the projected range of ions implanted into the semiconductor <R _p (E)> ^(2), a known projection of considering 1 order term <R _p (E) > When using ⁽¹⁾ , the perturbation term Δ _p ⁽²⁾ (E)
< _Rp (E)> ⁽²⁾ = < _Rp (E)> ⁽¹⁾ + _Δp ⁽²⁾ (E)
A method for generating ion implantation, characterized in that it is obtained using a second-order perturbation model using the approximate expression given above. The deviation ΔR _p (E) ⁽²⁾ considering the second order term of the range R _p of ions to be implanted into the semiconductor is _{defined as} the lateral extent of the known range R _p taking into account the second order term <R _{When T} ² (E)> ⁽²⁾ ,
<R _c ² (E)> ≡ <R _p ² (E)> + <R _T ² (E)>, and
<R _r ² (E)> ≡ <R _p ² (E)>-<R _T ² (E)> / 2
<R _c ² (E)> and <R _r ² (E)> defined in the above, and known <R _c ² (E)> ⁽¹⁾ , <R _r ² (E)> ⁽¹⁾ , perturbation term Δ _c ^{2 (2)} (E) and perturbation term Δ _r ^{2 (2)} (E),
<R _c ² (E)> ⁽²⁾ = <R _c ² (E)> ⁽¹⁾ + Δ _c ^{2 (2)} (E) and
<R _r ² (E)> ⁽²⁾ = <R _r ² (E)> ⁽¹⁾ + Δ _r ^{2 (2)} (E)
Approximate to the ^{_{<R c 2 (E)>}} (2), <R r 2 (E)> (2) and the ^{_{<R p 2 (E)>}} (2) using <ΔR _p ² (E)> ⁽²⁾ = {<R _c ² (E)> ⁽²⁾ +2 <R _r ² (E)> ⁽²⁾ } / 3
− <R _p ² (E)> ⁽²⁾
The ion implantation generation method according to claim 1, wherein the ion implantation is obtained using a second-order perturbation model using the approximate expression. Projection <R _p (E)> ⁽³⁾ taking into account the third-order term of the range R _p of the range of ions implanted into the semiconductor, the above <R _p (E)> ⁽²⁾ and the perturbation term as Δ _p ^{(3) Using} (E),
< _Rp (E)> ⁽³⁾ = < _Rp (E)> ⁽²⁾ + _Δp ⁽³⁾ (E)
The ion implantation generation method according to claim 1, wherein the ion implantation is obtained using a third-order perturbation model using the approximate expression. The deviation ΔR _p (E) ⁽³⁾ taking into account the third order term of the range R _p of ions to be implanted into the semiconductor is _{defined as} the lateral extent of the known range R _p taking into account the third order term <R _{When T} ² (E)> ⁽³⁾ ,
<R _c ² (E)> ≡ <R _p ² (E)> + <R _T ² (E)>, and
<R _r ² (E)> ≡ <R _p ² (E)>-<R _T ² (E)> / 2
<R _c ² (E)> and <R _r ² (E)> defined in the above, and known <R _c ² (E)> ⁽²⁾ , <R _r ² (E)> ⁽²⁾ , perturbation term Δ _c ^{2 (3)} (E) and perturbation term Δ _r ^{2 (3)} (E),
<R _c ² (E)> ⁽³⁾ = <R _c ² (E)> ⁽²⁾ + Δ _c ^{2 (3)} (E) and
^{_{<R r 2 (E)>}} (3) = <R r 2 (E)> (2) + Δ r 2 (3) (E)
Approximate to the ^{_{<R c 2 (E)>}} (3), <R r 2 (E)> (3) wherein the ^{_{<R p 2 (E)>}} (3) using a <ΔR _p ² (E)> ⁽³⁾ = {<R _c ² (E)> ⁽³⁾ +2 <R _r ² (E)> ⁽³⁾ } / 3
− <R _p ² (E)> ⁽³⁾
The ion implantation generation method according to claim 3, wherein the ion implantation is obtained using a third-order perturbation model using the approximate expression described above. The lateral spread parameter <ΔR _pt ² > ⁽³⁾ taking into account the third order term of the range R _p of ions to be implanted into the semiconductor, using the above < _RT ² (E)> ⁽³⁾ ,
<ΔR _pt ² > ⁽²⁾ = < _RT ² (E)> ⁽²⁾ / 2
The ion implantation distribution generation method according to claim 3, wherein the ion implantation distribution is obtained using a second-order perturbation model. A simulator incorporating a calculation formula relating to the ion implantation distribution generation method according to any one of claims 1 to 5.

Images