A method is described to estimate the temperature dependent interaction between two uncharged point defects in Si based on DFT calculations. As an illustration, the formation of the uncharged di-vacancy *V _{2}* is discussed, based on the temperature dependent attractive field between both vacancies. For that purpose, all irreducible configurations of two uncharged vacancies are determined, each with their weight given by the number of equivalent configurations. Using a standard 216-atoms supercell, nineteen irreducible configurations of two vacancies are obtained. The binding energies of all these configurations are calculated. Each vacancy is surrounded by several attractive sites for another vacancy. The obtained temperature dependent of total volume of these attractive sites has a radius that is closely related with the capture radius for the formation of a di-vacancy that is used in continuum theory. The presented methodology can in principle also be applied to estimate the capture radius for pair formation of any type of point defects.

## I. INTRODUCTION

Silicon crystals grown by the Czochralski technique can contain voids leading to so-called crystal originated pits (COP’s) on polished wafer surfaces. These voids that are also observed in Czochralski-grown germanium crystals, are large vacancy clusters that are formed during cooling immediately after solidification.^{1} COP’s have a detrimental impact on gate oxide integrity and can cause also problems with epitaxial layer quality.^{2–4} Besides intrinsic point defect pairs and clusters, also a wide variety of impurity-intrinsic point defect pairs and clusters are known.^{5} In continuum theory, point defect pair formation kinetics are described by rate equations, assuming a temperature independent capture radius.^{6} The capture radius that is used in these rate equations, is mostly chosen arbitrarily in order to have a good agreement with experimental results. Furthermore, isotropic diffusion is mostly assumed in the theory of diffusion-limited reactions of reacting point defects in a solid,^{6} while in real crystals the point defect diffusivity depends on the crystallographic directions.

The reaction rate *k _{f}* for the formation of a point defect pair AB, is commonly written as:

where *r* is the capture radius of point defect A and B and *D _{A}* and

*D*are the diffusion coefficients of A and B, respectively. The capture radius is usually assumed to be of the order of the bond length of Si with values depending on the type of point defect pair.

_{B}^{7–10}0.5 nm is a value that is typically used to describe the formation kinetics of the di-vacancy.

^{7,10}

In the past, material related phenomena have mostly been described using continuum models with empirical model parameters, due to a lack of atomistic models and calculation methods. Nowadays, ab initio calculations can be used to reveal the atomistic details and processes behind the continuum models. In a previous paper,^{11} an atomistic picture of the diffusion of two uncharged vacancies forming a di-vacancy in Si was discussed, based on DFT calculations using 64-atoms supercells. Due to the 64 atoms supercell, it was not possible to investigate sufficiently accurate the long-range interaction between the two vacancies which is needed for the determination of the interaction volume between the two point defects.

In the present paper, the interaction of two vacancies is investigated by using a 216-atoms supercell for the DFT calculations, which allows calculating with high accuracy the temperature dependent radius of the attractive field surrounding the vacancy which is a good and non-empirical approximation of the capture radius used in continuum theory.

## II. CALCULATION DETAILS

### A. Irreducible point defect pair configurations

A program code was used allowing determining all basic irreducible atomic configurations of two vacancies in the 216-atoms supercell. The weight of each configuration is determined as the number of equivalent configurations and this is done by considering the symmetry of the Si crystal and the periodic condition of the supercell. This search is performed automatically by using a script program code based on a published algorithm,^{12} and is similar to a method used elsewhere.^{13} All the atomistic models are automatically produced, not including the local distortions although there is of course an influence of the initial local distortion of the position of atoms around the vacancies due to the Jan-Teller effect.^{14} The local distortions are neglected as they have only a small effect on the calculated interaction volume around the vacancy. The main goal of the present paper is indeed not a detailed study of the silicon di-vacancy but rather to illustrate an approach that can be used for all types of point defect pairs. When studying more in detail a specific point defect pair, such details can be taken into account in the ab initio calculations whenever considered useful.

It turns out that there are nineteen irreducible configurations of two vacancies in the 216-atoms cell while in a 64-atoms cell there are only nine.^{11} These configurations have been automatically extracted from the total numbers of 23220 (= _{216}C_{2}) possible configurations for a 216-atoms supercell. In some of these configurations, two vacancies are on a same zigzag bond along a <110 > direction in a {110} plane and have a stronger interaction similar as between the vacancy and an impurity atom.^{15}

### B. DFT calculation procedure

The DFT calculations are based on the standard approach, using the local density approximation^{16,17} with the ultrasoft pseudopotential method,^{18} and plane waves as basis set for efficient structure optimization. The generalized gradient approximation (GGA) is used for describing the exchange–correlation energy, in the Perdew–Burke–Ernzerhof (PBE) form.^{19} The *CASTEP* code is used to self-consistently solve the Kohn-Sham equation using a three-dimensional periodic boundary condition.^{20} The density mixing method^{21} and the BFGS geometry optimization method^{22} are used to optimize the electronic structure and the atomic configurations, respectively. Only the neutral charge state of the system, which corresponds to vacancies introduced in nearly intrinsic Si, is considered in the present study. Even in more heavily doped Si this is a useful approximation in the early stage of cooling after solidification during Si single crystal growth from a melt. For k-point sampling, the 2 × 2 × 2 special points of the Monkhorst-Pack grid^{23} are used. The cutoff energy of the plane waves is 350 eV. Under these conditions, the supercell of the perfect Si was optimized and is a cube with a 1.639 nm size. Each cell with these nineteen basic configurations of two vacancies is then created and relaxed fully. Furthermore, the binding energy of the two vacancies is calculated for each configuration by using the formation energy of an uncharged isolated *V* which is 3.58 eV for a 216-atoms supercell as calculated by the standard approach based on the Si atom self-energy in the perfect cell. Although the energy differences between this size cell and 64 atoms cell^{11} is about 0.3 eV for the second vacancy at positions #1 or #3, the sites at a distance in the range 0.6-0.8 nm, except site #6, #9 which are on the cell boundary in a 64 atom cell, have converged to values differing less than 0.1 eV from those obtained with the 216 atom supercell.

## III. RESULTS AND DISCUSSION

### A. Binding energy between two vacancies

Fig. 1(a) shows the relative configuration of two vacancies whereby one vacancy is approaching the second one which is assumed to be fixed at the position labeled #0 in Fig. 1(a). Figure 1(b) shows the binding energies of the nineteen basic configurations of two vacancies in a 216-atoms supercell as a function of the nominal distance between them. The different vacancy positions in Fig. 1(b) are labeled by the numbers corresponding to the numbers in Fig. 1(a). In Fig. 1(b), each *V* position is defined by the position of the Si atom before forming *V* and before structural relaxation. The most stable configuration is labeled #1. This configuration describes two vacancies located at the nearest neighbor site, i.e., the di-vacancy. The obtained binding energy *E _{b}* is 1.75 eV, significantly higher than the 1.43 eV obtained before for a 64-atoms supercell.

^{11}Both values are lower than the calculated 2 eV reported in literature

^{24,25}but are close to the reported experimental values of 1.5 eV

^{26}and ≳1.6 eV.

^{27}The second

*V*set at position #2 as an initial model of the atoms in 216-atoms cell moved to the #1 position (forming a di-vacancy) after the geometry optimization. There was a meta-stable position for the second

*V*at #2 and no trasition state was found between #1 and #2 positions for the second

*V*in 64-atoms cell.

^{11}The 216-atoms cell size is large enough to estimate the attractive volume size since the binding energy of the second

*V*at #9 become negative thus limiting the number of binding sites. Comparing to the 64-atoms supercell, the binding energies of the second

*V*at #3∼6 or #7-1 increased.

### B. Repulsion and attraction between two vacancies

Figure 1(a) shows the binding energies between one *V* fixed at the position #0 and a second *V* which is mobile. The fixed *V* is surrounded by attractive sites (blue or light blue circles), and by repulsive sites (red or light red circles). In some configurations like #0-#7-2, #0-#13 and #0-#18, as shown in Fig. 2, the two vacancies have a negative binding energy smaller than -0.45 eV. It is interesting to note that configurations #7-1 and #7-2 have the same distance between the two vacancies but have totally different binding energies. Since the sites that have a negative binding energy (= are repulsive) are limited in number, the effect of these sites is ignored in the following when estimating the attractive interaction volume around a vacancy. A second vacancy at a site with negative binding energy will indeed move to a neighboring site with a positive binding energy.

### C. The impact of the attractive and repulsive fields on vacancy pairing

An isolated *V* has a migration energy of about 0.25 eV.^{28} When this moving *V* enters the interaction volume of another vacancy with several attractive sites (blue or light blue circles in Fig. 1(a)), the moving vacancy will slow down by being temporarily trapped in sites with low binding energy. When the trapped vacancies are released again, they will preferentially move to sites with larger binding energy and finally end up in the most stable configuration being the well-known di-vacancy. The range of the attractive field surrounding a vacancy, is therefore closely related to the capture cross-section or capture radius used in continuum modeling. The volume of this field (which is not spherical due to its dependence on crystallographic directions) can be converted into an attractive field radius, which is an approximation of the capture radius by averaging as discussed further below. In a previous paper,^{11} a first estimation of this volume was obtained using a 64-atom supercell which is however too small to contain to complete attractive volume for which an at least 216-atoms supercell (4.40 nm^{3}).

### D. The attractive field at 0 K

The extent of the attractive field around vacancy for a second vacancy is determined by the *E _{b}* > 0 sites listed in Table I. Table II shows that when assuming 0 = <

*E*<

_{b}*E*, all sites except #17 can contribute to the attractive field surrounding site #0. The ratio of the volume of the attractive field and the volume of the supercell is estimated by calculating the ratio of the summed weights of the attractive sites and the total number of sites in the supercell. Hereby the relaxations of atoms around the vacancies are neglected after the relaxation calculations. Using criterion (a) of Table II, yields a radius

_{b1}*r*of the attractive field at 0 K of 0.74 nm. Using the criteria (b) or (c) of Table II to define the radius of the attractive field, yields values of 0.70 and 0.68 nm, respectively. This gives an idea of the uncertainty on

*r*, depending the criterion one uses to define the boundary of the attractive field.

Site# . | Symbol E
. _{bi} | E (eV)
. _{b} |
---|---|---|

1 | E _{b6} | 1.75 |

2 | E _{b6} | (1.75)^{a} |

3 | E _{b5} | 0.25 |

5 | E _{b4} | 0.16 |

6 | E _{b3} | 0.12 |

4 | E _{b2} | 0.073 |

7-1 / 17 | E _{b1} | ≈0.055^{b} |

Site# . | Symbol E
. _{bi} | E (eV)
. _{b} |
---|---|---|

1 | E _{b6} | 1.75 |

2 | E _{b6} | (1.75)^{a} |

3 | E _{b5} | 0.25 |

5 | E _{b4} | 0.16 |

6 | E _{b3} | 0.12 |

4 | E _{b2} | 0.073 |

7-1 / 17 | E _{b1} | ≈0.055^{b} |

^{a}

The second V moved to the #1 position (forming a di-vacancy) after the geometry optimization.

^{b}

Values are very close to each other.

. | weight . | Binding . | (a)0 = < E < _{b}E_{b1}
. | (b) E_{b1} = < E < _{b}E_{b2}
. | (c) E_{b2} = < E < _{b}E_{b3}
. | (d) E_{b3} = < E < _{b}E_{b4}
. | (e) E_{b4} = < E < _{b}E_{b5}
. | (f) E_{b5} = < E < _{b}E_{b6}
. | (g) E_{b6} < E
. _{b} | |||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|

Second V site #
. | (23436 total) . | energy $ E b i $ (eV) . | Attractive ? . | Connected to site #0 ? . | Attractive ? . | Connected to site #0 ? . | Attractive ? . | Connected to site #0 ? . | Attractive ? . | Connected to site #0 ? . | Attractive ? . | Connected to site #0 ? . | Attractive ? . | Connected to site #0 ? . | Attractive ? . | Connected to site #0 ? . |

0 | 216 | - | (core) | Y | (core) | Y | (core) | Y | (core) | Y | (core) | Y | (core) | Y | (core) | Y |

1 | 432 | 1.75 | Y | Y | Y | Y | Y | Y | Y | Y | Y | Y | Y | Y | Y | Y |

2 | 1296 | (1.75) | Y | Y | Y | Y | Y | Y | Y | Y | Y | Y | Y | Y | N | N |

3 | 1296 | 0.25 | Y | Y | Y | Y | Y | Y | Y | Y | Y | Y | N | N | N | N |

4 | 648 | 0.07 | Y | Y | Y | Y | N | N | N | N | N | N | N | N | N | N |

5 | 1296 | 0.16 | Y | Y | Y | Y | Y | Y | Y | Y | N | N | N | N | N | N |

6 | 2592 | 0.12 | Y | Y | Y | Y | Y | Y | N | N | N | N | N | N | N | N |

7-1 | 1296 | 0.06 | Y | Y | N | N | N | N | N | N | N | N | N | N | N | N |

7-2 | 432 | -0.46 | N | N | N | N | N | N | N | N | N | N | N | N | N | N |

9 | 1296 | -0.10 | N | N | N | N | N | N | N | N | N | N | N | N | N | N |

10 | 2592 | -0.08 | N | N | N | N | N | N | N | N | N | N | N | N | N | N |

11 | 1296 | -0.13 | N | N | N | N | N | N | N | N | N | N | N | N | N | N |

12 | 1296 | -0.27 | N | N | N | N | N | N | N | N | N | N | N | N | N | N |

13 | 864 | -0.47 | N | N | N | N | N | N | N | N | N | N | N | N | N | N |

14 | 1296 | -0.03 | N | N | N | N | N | N | N | N | N | N | N | N | N | N |

15 | 2592 | -0.06 | N | N | N | N | N | N | N | N | N | N | N | N | N | N |

16 | 1296 | -0.16 | N | N | N | N | N | N | N | N | N | N | N | N | N | N |

17 | 324 | 0.06 | Y | N | N | N | N | N | N | N | N | N | N | N | N | N |

18 | 432 | -0.50 | N | N | N | N | N | N | N | N | N | N | N | N | N | N |

19 | 648 | -0.18 | N | N | N | N | N | N | N | N | N | N | N | N | N | N |

Attractive & connected to #0 | 9072 | 7776 | 7128 | 4536 | 3240 | 1944 | 648 | |||||||||

Volume of attractive field (nm^{3}) | 1.71 | 1.46 | 1.34 | 0.85 | 0.61 | 0.37 | 0.12 | |||||||||

Maximum separation r (nm) _{i} | 0.74 | 0.70 | 0.68 | 0.59 | 0.53 | 0.44 | 0.31 |

. | weight . | Binding . | (a)0 = < E < _{b}E_{b1}
. | (b) E_{b1} = < E < _{b}E_{b2}
. | (c) E_{b2} = < E < _{b}E_{b3}
. | (d) E_{b3} = < E < _{b}E_{b4}
. | (e) E_{b4} = < E < _{b}E_{b5}
. | (f) E_{b5} = < E < _{b}E_{b6}
. | (g) E_{b6} < E
. _{b} | |||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|

Second V site #
. | (23436 total) . | energy $ E b i $ (eV) . | Attractive ? . | Connected to site #0 ? . | Attractive ? . | Connected to site #0 ? . | Attractive ? . | Connected to site #0 ? . | Attractive ? . | Connected to site #0 ? . | Attractive ? . | Connected to site #0 ? . | Attractive ? . | Connected to site #0 ? . | Attractive ? . | Connected to site #0 ? . |

0 | 216 | - | (core) | Y | (core) | Y | (core) | Y | (core) | Y | (core) | Y | (core) | Y | (core) | Y |

1 | 432 | 1.75 | Y | Y | Y | Y | Y | Y | Y | Y | Y | Y | Y | Y | Y | Y |

2 | 1296 | (1.75) | Y | Y | Y | Y | Y | Y | Y | Y | Y | Y | Y | Y | N | N |

3 | 1296 | 0.25 | Y | Y | Y | Y | Y | Y | Y | Y | Y | Y | N | N | N | N |

4 | 648 | 0.07 | Y | Y | Y | Y | N | N | N | N | N | N | N | N | N | N |

5 | 1296 | 0.16 | Y | Y | Y | Y | Y | Y | Y | Y | N | N | N | N | N | N |

6 | 2592 | 0.12 | Y | Y | Y | Y | Y | Y | N | N | N | N | N | N | N | N |

7-1 | 1296 | 0.06 | Y | Y | N | N | N | N | N | N | N | N | N | N | N | N |

7-2 | 432 | -0.46 | N | N | N | N | N | N | N | N | N | N | N | N | N | N |

9 | 1296 | -0.10 | N | N | N | N | N | N | N | N | N | N | N | N | N | N |

10 | 2592 | -0.08 | N | N | N | N | N | N | N | N | N | N | N | N | N | N |

11 | 1296 | -0.13 | N | N | N | N | N | N | N | N | N | N | N | N | N | N |

12 | 1296 | -0.27 | N | N | N | N | N | N | N | N | N | N | N | N | N | N |

13 | 864 | -0.47 | N | N | N | N | N | N | N | N | N | N | N | N | N | N |

14 | 1296 | -0.03 | N | N | N | N | N | N | N | N | N | N | N | N | N | N |

15 | 2592 | -0.06 | N | N | N | N | N | N | N | N | N | N | N | N | N | N |

16 | 1296 | -0.16 | N | N | N | N | N | N | N | N | N | N | N | N | N | N |

17 | 324 | 0.06 | Y | N | N | N | N | N | N | N | N | N | N | N | N | N |

18 | 432 | -0.50 | N | N | N | N | N | N | N | N | N | N | N | N | N | N |

19 | 648 | -0.18 | N | N | N | N | N | N | N | N | N | N | N | N | N | N |

Attractive & connected to #0 | 9072 | 7776 | 7128 | 4536 | 3240 | 1944 | 648 | |||||||||

Volume of attractive field (nm^{3}) | 1.71 | 1.46 | 1.34 | 0.85 | 0.61 | 0.37 | 0.12 | |||||||||

Maximum separation r (nm) _{i} | 0.74 | 0.70 | 0.68 | 0.59 | 0.53 | 0.44 | 0.31 |

### E. The temperature dependent attractive field

With increasing temperature, dissociation of the pair of point defects into two separate point defects will become more and more important. In continuum theory, the reaction rate *k _{r}* of this reverse pairing or dissociation reaction is expressed as:

k_{f} is the reaction rate of the pairing process, while *C _{A}*,

*C*,

_{B}*C*are the concentrations of the point defects

_{AB}*A*,

*B*and of the point defect pair

*AB*, respectively. In the atomistic picture of the migration of the second

*V*, the pairing and dissociation reactions always occur successively with the probability of both reactions determined by the activation energy of the point defect for the migration. Both reactions can be described by probabilities, which have to be introduced into the concept of the interaction field in the atomistic picture.

Therefore, in this paper, the interaction field is determined using the capture probability *p* of the second mobile *V* by the fixed one at #0, having a binding energy *E _{b}*. This binding energy is defined with respect to the energy of the two vacancies at infinite distance from each other. Figure 3 shows schematically the relation between

*p*and

_{i}*E*assuming an activation energy

_{bi}*E*for diffusion of the second

_{ai}*V*at the site with binding energy

*E*. Here, this activation energy

_{bi}*E*takes into account all barriers, including the repulsive sites around the attractive sites in Fig. 1(a), for the second

_{ai}*V*to diffuse from infinity to the site with a binding energy

*E*. Assuming that the vacancy concentration, which can be applied for eq. (2), around the barriers is uniform, the prefactors describing the inbound and out bound diffusion can be assumed to be the same and the capture probability

_{bi}*p*can be written as 1-exp[-

_{i}*E*/

_{bi}*k*]. When the temperature

_{B}T*T*approaches to 0 K,

*p*becomes 1 (100%) in case of

_{i}*E*> 0 as discussed above. The capture probability

_{bi}*p*on the site with binding energy

_{i}*E*depends on the binding energy

_{bi}*E*in Table I as

_{bi}*p*= 1-exp[-

_{i}*E*/

_{bi}*k*]. Figure 4(a) shows the

_{B}T*p*dependences on

_{i}*T*for sites having each binding energy

*E*, which classified the nominal capture radius

_{bi}*r*in Table II. The reaction rate

_{i}*k*in eq.(1) would be over-estimated when using

*r*. This is because

_{i}*r*is deduced under the condition

_{i}*p*< 1. To avoid this problem, we evaluated

_{i}*r*in eq. (1) as

*p*. Figure 4(b) shows the calculated

_{i}r_{i}*p*for each binding energy

_{i}r_{i}*E*(i.e., each second

_{bi}*V*site) in Fig. 4(a). As we consider all

*p*,

_{i}*p*has a different value for each

_{i}r_{i}*p*at each temperature. With increasing temperature, all

_{i}*p*decrease and lead to lower

_{i}*p*. When

_{i}r_{i}*p*starts to decrease,

_{i}r_{i}*p*

_{i+1}

*r*

_{i+1}is larger than

*p*. This means that with increasing temperature, the second vacancy will be captured at sites closer to the first vacancy and will form a vacancy pair there with high probability due to the higher binding energy as shown in Fig. 4(b). Therefore, the black curve, which is the maximum of these values for each temperature in Fig. 4(b), is the most probable radius of the effective attractive field and a good approximation for

_{i}r_{i}*r*in eq. (1). It may, however, be an overestimate considering all other possible geometries of the second vacancy after the

*p*process mentioned above. The real value of the radius is therefore expected to lie between the curves corresponding with the maximum and the averaged r values, also shown in Fig. 4(b).

_{i}The radius of the attractive field is nearly constant at temperatures below 100 K. At these low temperatures, the thermal energy becomes so small that it drops below the 0.1 eV that we assumed as a kind of cut off energy due to the inaccuracy of the ab initio calculation. This leads to an apparent constant radius.

A fourth order polynomial fit (R^{2} = 0.998) of the initio calculated values of the attractive field radius *r _{att}* corresponding with the maximum

*p*values between 0 K and melting temperature is shown in Fig. 4(c). The thus obtained polynomial fit

_{i}r_{i} can in first order approximation replace *r* in Eq. (1) between 0 K and melting temperature. Eq. (3) yields a value of 0.69 nm at room temperature, which is close to the capture radius for a di-vacancy used in the continuum models. The ab initio approach has the advantage that is based on an atomistic model, without making empirical choices. Furthermore it describes the temperature dependence of *r* which should allow more accurate modeling of rate equations.

As illustrated in Fig. 1(b), the attractive field for the second *V* contains also repulsive sites, in which some sites even have a negative binding energy smaller than -0.45 eV depending on the crystallographic directions around the fixed *V*. The presented methodology can in principle also be applied for calculating the interaction field leading to pair formation of any type of point defects, also those that have a longer interacting range than the 0.6 nm for *V*_{2} formation^{29,30} in which case also a larger supercell should be used for the calculations. Longer interaction ranges will in most cases also lead to larger capture radii.

## IV. CONCLUSIONS

A method is proposed to estimate the radius of the interaction field for the formation of a point defect pair in Si based on an atomistic picture. This radius is closely related to the capture radius used in continuum theory. The application of the method is illustrated by the simplest case which is the formation of a di-vacancy by pairing of two uncharged vacancies. Describing the point defect association and dissociation reactions by probabilities, the temperature dependence of the radius *r* of the interaction field can be obtained showing that it gradually decreases as the temperature increases as should therefore also be the case for the capture radius in the continuum model.

## ACKNOWLEDGMENT

This work is partially supported by the Advanced Low Carbon Technology Research and Development Program (ALCA) of the Japan Science and Technology Agency (JST). We are also indebted to Dr. Abhijit Chatterjee of Accelrys Software Inc., for stimulating discussions and for developing the script program code that is used in the present study to calculate the irreducible vacancy configurations.

## REFERENCES

*CASTEP code*is available from Accelrys Software Inc.