Java Simulation of Au Diffusion in Si Affected by Vacancies and Self-Interstitials : Partial Differential Equations

A Java program in a GUI environment has been developed for the numerical solution of basic partial differential equations and applied to Au diffusion in Si affected by vacancies and self-interstitials. Text fields of selected parameters for the calculation are set on the display, and the calculation starts by checking the start button after putting values on the text fields. The calculated results are plotted immediately after the finish of the calculation as the concentration profiles of substitutional Au, interstitial Au, vacancies and self-interstitials, and their diffusion can be presented immediately, resulting in the identification of the diffusion mechanism. By changing the values of the text fields, new results can be represented immediately. The diffusion of Au in Si can be simulated correctly and easily by this program. Results from the program for one set of conditions are shown, including images produced on the display.


Introduction
Au atoms in Si occupy interstitial and substitutional sites, and the substitutional Au exists in three states depending on the heat treatment history [1]: high-temperature substitutional Au, low-temperature substitutional Au, and agglomerations of substitutional Au.During the heat treatment, high-temperature substitutional Au diffuses very slowly itself and the change in its concentration is dominated by an interchange mechanism with interstitial Au and substitutional Au associated with vacancies [2] and self-interstitials [3].In this case, the concentrations of substitutional Au, interstitial Au, vacancies and self-interstitials can be obtained from four partial differential equations [4,5].The concentration of substitutional Au exhibits a diffusion profile depending on the relative contributions of vacancies and interstitial Si atoms [6].The author has previously investigated Au diffusion using Java programming by obtaining a numerical solution to a single partial differential equation obtained from the above four differential equations under several approximations.This approach was adopted because of the difficulty of directly numerically solving the above four partial differential equations due to the limitation in the capacity of personal computers [7].However, recently the capacity of personal computers has been progressing rapidly, and the author has been able to directly solve the four partial differential equations involved in Au diffusion by Java programming.As a result, the diffusion of Au in Si can be simulated correctly and easily using Java in a GUI (Graphical User Interface) environment.

Basic Partial Differential Equations for Au Diffusion in Si
The basic diffusion equations for substitutional impurities, interstitial impurities, vacancies and self-interstitials in the interchange mechanism of interstitial and substitutional impurities associated with vacancies and self-interstitials are given as Here, the terms caused by the internal sink-sources such as dislocations in Equations ( 2)-( 5) based on the site conservation are added to the equations given in reference [4].In above equations, N, t, x, K and D are concentration, time, distance from surface, chemical reaction constant and diffusion constant.The subscripts s , i , V , I , Si , Vr , Ve , Ir , Ie , Ls , Fr , Fe , iSS , VSS, and 0 stand for substitutional impurity, interstitial impurity, vacancy, self-interstitial, Si atom at a lattice site, vacancy recombination, vacancy emission, self-interstitial recombination, self-interstitial emission, lattice site, recombination of Frenkel pair, emission of Frenkel pair, emission and annihilation of interstitial impurity at internal sink-sources, emission and annihilation of vacancy at internal sink-sources, emission and annihilation of self-interstitial at internal sink-sources, and thermal equilibrium value.Each concentration is normalized by its thermal equilibrium concentration and the distance, x, is normalized by the specimen thickness, L, for easier numerical calculation.


  Here, C is the normalized concentration and ξ is x/L.The coefficients are given as follows.

Crank-Nicolson's Implicit Method and Gauss-Seidel's Iteration Method
The partial differential terms in the equations are approximated to the difference equations by Crank-Nicolson's implicit method.Namely, an advanced differential-difference approximation is used for the partial differential by time, ∂C/∂t, and a neutral differential-difference approximation is used for the partial differential by distance, neutral value for the function of C, f(C), according to Crank-Nicolson. .
 , where e is . We obtain a criterion value for the convergence, and (39) used as an approximate value of .

Constants
Calculating the new value of a variable, such as   i, j 1 C  , involves using unknown values of variables such as in Equations ( 35)-(38).Here, we use Gauss-Seidel iteration method, in which the previously calculated value is used as the initial value in the first calculation of from i = 0 to i = i max , for example, .We note the value obtained from the The chemical reaction constants for the recombination reactions based on random diffusion to sinks are given by Damask and Dienes [8].
  first calculation as .For the second calculation of C  , , the first calculated values are used as

C
Here, r iV , r sI , and r VI are recombination distances be-tween interstitial impurity and vacancy, substitutional impurity and self-interstitial, and vacancy and self-interstitial.The chemical reaction constants for the creation reactions are given from thermal equilibrium conditions as follows, The chemical reaction constants of the internal sinksources, are given also by Damask and Dienes, for spherical internal sink-sources with a capture radius of r c and a concentration of N c , and as for cylindrical internal sink-sources with a density of n d .
We used a single atomic distance of 0.234 nm for the recombination distances and the capture radius, by assuming that the recombination and the capture occur for the nearest neighbours.The thermal equilibrium concentration of substitutional Au in Si is given by Collins, Carlson and Gallagher [9]   and the diffusion constant and thermal equilibrium concentration of interstitial Au in Si is given by Willcox and LaChapelle [10]   24 3 i0 2.52 5.95 10 exp cm , N kT Values of the thermal equilibrium concentrations and diffusion constants of vacancies and self-interstitials differ greatly between different authors, by as much as 6 orders of magnitude [11].Therefore, it is difficult to choose particular values.On the other hand, comparatively consistent products of the equilibrium concentration and the diffusion constant are available [12].The self-diffusion of host atoms is given by summing two vacancy and self-interstitial terms [11].In our case [6], take into account the correlation factors.Here, D self is the self-diffusion constant of Si atom, whose value does not differ much among the authors [13], and we use a D times value of the equation by Demond et al. [14], where 0.1 ≤ a D ≤ 10.We chose the contribution rate of self-interstitials, d Isd , as a constant related to the concentrations of vacancies and self-interstitials, where, f I is the correlation factor of the self-interstitials.
Then the contribution rate of vacancies, d Vsd , is obtained as where, f V is the correlation factor of the vacancies.In our case, f I = 0.7273 and f V = 0.5.A value of either one of the equilibrium concentrations or a diffusion constant is necessary in order to fix both of them from Equations ( 57) and (58).It is better to use a diffusion constant as it has a true physical nature, whereas fixing a equilibrium concentration at some arbitrary value can cause difficulties, such as unrealistically small diffusion constant if that contribution is in reality very small.We use a V times value of by Masters and Gorey [15], and a I times value of by Tan and Gösele [16], where 0.1 ≤ a V ≤ 10 and 0.1 ≤ a I ≤ 10.N Si0 is obtained by Equations ( 6) and (58), N I0 and N V0 are obtained by Equations ( 57) and (58), respectively.

Selection of Parameters
The calculated results can be presented immediately as figures by changing the value of the text fields.We select T, L, d Isd , c r and j max as the values of the text fields, where c r indicates the ratio of k to h, and 2 r . 2 Here, c r = 1 is the condition of convergence in the explicit method.Such a limitation for the convergence is not necessary in the implicit method, but the calculation does not converge in the Gauss-Seidel iteration if we use too large a value of c r .It is better to use a larger value of c r to provide a shorter computing time until a required diffusion time is achieved and then it is better to choose a larger convergence value after a tentative calculation with a small j max using several c r values.

Boundary and Initial Conditions
The surface concentration of impurities during heat treatment for impurity diffusion differs depending on the surface conditions [17].Therefore, we adopt a surface condition with sufficient impurities to put the boundary conditions of the surface into their equilibrium concentrations, and the surface concentrations of vacancies and self-interstitials should be at their equilibrium concentrations due to the presence of sufficient sink-sources of the surface.The boundary conditions at the other back surface are the same as the top surface mentioned above, and we put i = i max at x = L/2 and the concentrations at max is equal to those at max .We use the boundary conditions such as the thermal equilibrium concentration at i = 0 and the same concentration at and i i .
The initial concentrations of impurities contributing to impurity indiffusion depend on the quality of the asgrown crystal, and we used 5 × 10 10 cm -3 as the initial concentration of substitutional impurities assuming the quality as 99.9999999999%.We used 5.0 × 10 10 × N i0 /N s0 cm -3 as the initial concentration of interstitial impurities.The initial concentration of substitutional impurities contributing to out-diffusion can be set to the measured value after the pre-indiffusion of the impurities and that of interstitial impurities is set to the equilibrium concentration at the pre-indiffusion temperature.The initial concentrations of vacancies and self-interstitials in impurity indiffusion depend on the growth rate and on the longitudinal temperature gradient near the crystallization front [18], and we used N V ≤ 7.3 × 10 12 cm -3 , N I ≤ 2.6 × 10 13 cm -3 [11].We used N V = a × 7.3 × 10 12 cm -3 and N I = b × 2.6 × 10 13 cm -3 as the initial concentrations for impurity indiffusion, here, a ≤ 1 and b ≤ 1.Those for impurity out-diffusion can be set as the thermal equilibrium concentrations at the temperature for the pre-indiffusion heattreatment. 1  

Java Simulation
One example of the program used here is shown below.
The function of the main part of program is written as a comment line.Five text fields for temperature T, thickness L, fraction of interstitial mechanism d Isd , the ratio of k to h, c r , and number of calculations (to determine the time) j max are set on the display, and the calculation starts by checking the start button after putting the value in the text fields.The results are shown immediately as figures on the display after the end of the calculation, and the calculated results can be identified immediately.By changing the values of the text fields, the new results can be again be identified immediately.The distributions of the logarithmic concentrations of C s , C i , C V , C I and C Si are plotted on the display automatically in the range from (xx mini ,yy mini ) to (xx max ,yy max ).Six distributions at j = 0, j = j max /5, j = 2 × j max /5, j = 3 × j max /5, j = 4 × j max /5 and j = j max for each concentration are plotted in our case.One of the computed results for the indiffusion process is shown in Figure 1 as an image on the display, and that for the annealing is shown in Figure 2.

Discussion and Results
The author has numerically solved Equations ( 6)-( 10) by Java programming and the distribution of each concentration can be simulated correctly and easily using a usual personal computer employing the GUI method of Java.Namely, simultaneous partial different equations can be calculated numerically under known boundary and initial conditions by this program, and the distribution of each unknown can be easily simulated.The time for the computation until the required diffusion time depends on c r , the ratio of k to h, and its maximum value for the convergence depends on the constants in the equations.In our case (d Isd = 0.1 and x L = 0.1 cm), the maximum value is about 300 for the indiffusion at 900˚C and about 40 at 1160˚C.It is better to choose the largest c r value by tentative calculations with small j max values.

Figure 1 .
Figure 1.An example of the computed result for the indiffusion of Au in Si of 1 mm thickness at 900˚C.The figure is shown as an image on the display.

Figure 2 .
Figure 2.An example of the computed result for the annealing of supersaturated Au in Si of 1 mm thickness at 900˚C.The figure is shown as an image on the display.