Algorithm for Calculating the Initial Defect Structure of Semiconductor Silicon

An algorithm for calculating the defect structure of semiconductor silicon crystals was proposed. The proposed approach makes it possible to calculate the sizes, distribution densities of grown-in microdefects at any point of the crystal. Calculations are performed by using and analyzing the thermal conditions of crystal growth in the temperature range from 1683 K to 300 K. The algorithm flowchart is given.


Introduction
Semiconductor silicon includes dislocation-free single crystals of silicon, which were grown by the Czochralski (CZ-Si) and floating zone melting (FZ-Si). At present, without such crystals, it is impossible to imagine the development of the electronics industry. The parameters and properties of silicon devices depend critically on the defective structure of CZ-Si and FZ-Si crystals. Various technological treatments of silicon for the creation of devices lead to the transformation of the initual defect structure.
Recently, physical foundations for the formation of the initial defect structure of semiconductor silicon [1,2]. Mathematical models of the formation of various structural defects have also been developed [3]. This means that defect formation during crystal growth can be a completely controlled process. To do this, it is necessary to create a software package that could analyze the processes of defect formation in silicon from a unified position. Thus, the task of developing virtual methods for investigating (controlling) the initial defect structure of silicon with the help of algorithms of formation and transformation, both individual types of defects, and the entire structure, depending on the thermal conditions of crystal growth is relevant. The solution of this problem is the purpose of this work.

Physical Model of Formation and Transformation of Structural Defects
In the process of growth of dislocation-free single crystals of silicon are formed structural imperfections. They are called grown-in microdefects. It is well known that the formation and transformation of grown-in microdefects depends critically on the thermal conditions of the crystal growth [1,2]. The thermal conditions for crystal growth include: (i) growing method, (ii) crystal growth rate ( ), axial and radial temperature gradients ( and ), crystal cooling rate ( ), crystal diameter. Experimental studies of growth microdefects made it possible to establish the existence of three types of structural imperfections: precipitates of impurities [1,3], dislocation loops [4], and microvoids [5]. It was found that microvoids and dislocation loops are formed in different regions of the crystal, depending on the growth parameter ⁄ = (where is a constant) [6]. It was established that irrespective of the growth method near the crystallization front in unalloyed silicon single crystals precipitates of oxygen (SiO 2 ) and carbon (SiC) are formed. The defect formation in accordance with the diffusion model occurs in three stages [3]: (i) the formation of impurity complexes near the crystallization front; (ii) the formation, growth, and coalescence of precipitates during cooling of the crystal from the crystallization (melting) temperature to room temperature; and (iii) the formation of microvoids or dislocation loops at temperatures ≤ 1423 depending on the growth parameter ⁄ . This physical model is called the diffusion model of the formation of grown-in microdefects. The mathematical apparatus of the diffusion model using modern information technologies can provide the basis for the development of a software complex for analysis and calculation of the formation of growth and postgrowth microdefects in dislocation-free silicon single crystals [2].

Calculation of the Initial Defect Structure of Semiconductor Silicon
Proceeding from the diffusion model of formation and transformation of grown-in microdefects, the unified algorithm for calculating the defect structure of silicon should be a complex of algorithms for the formation of individual grown-in microdefects, taking into account the thermal conditions for growing crystals. The experimental and theoretical investigations of ultrapure dislocation-free silicon crystals demonstrated that the formation of a defect structure during the crystal growth occurs in the direction from the high-temperature precipitation of impurities to the formation of secondary growth defects (microvoids or dislocation loops). The main role in the process of defect formation is played by the precipitation of impurities [2]. Therefore, in the center of the unified algorithm is the algorithm for the formation of precipitates.

Algorithm for the Formation of Precipitates
At present, two approaches are used to calculate the initial stage of formation of precipitates (formation of impurityintrinsic defect complexes): classical and probabilistic [3,7].

Classic Approach
In [8], was calculated the formation of "impurity-intrinsic point defect" complexes during the growth of silicon crystals with the inclusion of the elastic interaction between the components of the complex. The calculation algorithm is based on approximation of the strong complexation for model of consistent diffusion, which allows a simple physical interpretation and can be the most adequate for the physical model. [9]. The solution to the problem of diffusion of a component A in a semi-infinite sample homogeneously doped with a component B, with the absence of evaporation of the component B from the sample and the presence of the free component A at the sample boundary, has the form [10]: ( ) where is the concentration of complexes; , , and are the diffusion coefficients of the free components , and complexes, respectively; is the coordinate (crystal length); and is time. The component is the background impurity (oxygen or carbon ) and the component is intrinsic point defects (vacancies or interstitials ). ! and ! are the concentrations of the free impurities and and " # and " # are the impurity concentrations at the interface [8]. In the calculations, the following quantities were used: 15 3 (0) The solution to equations has a physical meaning (" $ ≈ 10 $' − 10 $) cm -3 ) only at * = 0.01 . Note that, in the approximation of strong complex formation, * # is interpreted as the front boundary of the complex formation reaction. Since is the crystal length and # is the position of the crystallization front, we can conclude that complex formation occurs near the crystallization front [10].
This approximation is valid at the initial stages of the formation of nuclei, when their sizes are small and the use of Fokker-Planck continuity differential equations is impossible.
In paper [11] was considered the modern approach based on solving systems of coupled discrete differential equations of quasi-chemical reactions for the description of the initial stages of the formation of nuclei of new phases and a similar system of Fokker-Planck continuity differential equations. The nucleation and evolution of a complex system of grownin microdefects (which consists of oxygen and carbon precipitates) during cooling of the crystal are described by the systems of coupled differential equations for each type of defect: SiC In the system of equations (7), we took into account that the oxygen precipitates serves as sinks for oxygen atoms and vacancies and as sources of interstitial silicon atoms. Then, we can write the following equations: , SiO At the same time, the carbon precipitates, in turn, also serve as sinks for carbon atoms and interstitial silicon atoms and as sources for vacancies. Therefore, we can write In the general case, the proportionality factors , , , v i v i γ γ γ γ * * can depend on the quantities , o c n n and are determined by the conditions of thermodynamic equilibrium [12].
The corresponding system of coupled Fokker-Planck equations can be transformed into the dimensionless form: The normalized sizes of the precipitates are determined in the system of equations (10) as follows: SiC SiO SiC The fluxes of particles on the right hand sides of the system of equations (10) are described by the expressions in which the following notation is used for the normalized kinetic coefficients: The normalized rates of growth and dissolution of the precipitates in expressions (14) The critical size of the precipitates can be determined according to [12] from the expressions: are the supersaturations of the oxygen atoms, carbon atoms, intrinsic interstitial silicon atoms, and vacancies, respectively; σ is the density of the surface energy of the interface between the precipitate and the matrix; µ is the shear modulus of silicon; δ and ε are the relative linear and volume misfit strains of the precipitate and the matrix, respectively; i γ and v γ are the fractions of the intrinsic interstitial silicon atoms and vacancies per impurity atom attached to the precipitate, respectively; p V is the molecular volume of the precipitate; The number of impurity atoms in the compressed precipitates with the radii O r and C r is determined according to [13] from the formula: where is the fraction of impurity atoms per intrinsic defect, The calculations were performed using the following data: p V = 4,302·10 -2 nm 3 (SiO 2 ); p V = 2,04·10 -2 nm 3 (SiС); σ = 310 erg/cm 2 (SiO 2 ); σ = 1000 erg/cm 2 (SiС); µ = 6,41·10 10 Pa; δ = 0,3; ε = 0,15; SiO act The algorithm used for solving the problem of simulation of the simultaneous growth and dissolution of the oxygen and carbon precipitates due to the interaction of point defects during cooling of the crystal from the crystallization temperature is based on the monotonic explicit difference scheme of the first order accuracy as applied to the Fokker-Planck equations (10) [11].
The critical size distribution functions of the oxygen and carbon precipitates can be found only from the Fokker-Planck continuity differential equations. The model of dissociative diffusion makes it possible to analyze the processes occurring in the diffusion region near the crystallization front. These mathematical models, together with the experimental results obtained from the investigation of quenched crystals, have demonstrated that the nucleation processes occur very rapidly near the crystallization front.

Probabilistic Approach
The classical concepts of a periodic structure of a crystal are based on physical limitations, such as the localization of each atom in the vicinity of a fixed site of the crystal lattice, the consistency of introducing the concept of probability and the mechanical description of the behavior of particles, and the assumption that the number of atoms in a crystal is an integer. Meanwhile, in [14,15] Vlasov developed a different model of a crystal. In this model, the periodic structure of a crystal not is a consequence of the restrictions on the freedom of movement of atoms in the crystal, but is determined by specific statistical laws of motion of particles, in accordance with which the periodic structure agrees with the freedom of movement of atoms, so that the probability of finding an atom in interstitials is always different from zero [14].
The Vlasov model is based on the following basic physical concepts [14]: (i) the rejection of the principle of the spatial and velocity localization of particles (in terms of the classical mechanics), which takes place regardless of the force interactions; (ii) the introduction of force interactions by analogy with the classical mechanics, but with the inclusion of a new principle of nonlocalization of particles; and (iii) the behavior of each particle of the system is described by means of an ,-function extended in phase space.
The stationary properties of a crystal are described using the distribution density of particles -./0 = 1 ,./, 3043. The molecular field is determined by only probable, rather than exact, positions of the atoms. This is expressed by the potential function containing the probability density of the particles with the inclusion of the temperature distribution of the particles [14]. The choice of the pair interaction potential depends on the problem under consideration. Then, the nonlocal model of a crystal is based on the following nonlinear equations, which allow us to calculate the molecular potential and the particle position density under thermal equilibrium conditions [14]: considered to mean such a value of the parameter * for which the equations of system (19) have solutions other than trivial.
If the position of one of the particles is taken as the origin of the coordinates, we can determine -.00 = *5 . The characteristic number * can be determined from the main criterion for the existence of the crystalline. In this case, the crystallization condition can be written as follows: where " is the number of particles; G is the melting (crystallization) temperature of the crystal; $,' * = − $,' .
For the evaluation of the parameters of the formation of silicon-carbon and silicon-oxygen complexes, the interatomic interaction can be represented in the form of the Mee-Lennard-Jones potential: where and / # are the depth and the coordinate of the minimum of the interatomic potential, respectively; N and L are the parameters, N > L. In the case of the formation of a stable bond between the atoms, we have / = √2 The results of the calculations showed that the classical theory of nucleation and growth of second-phase particles in crystals and the Vlasov model of crystal formation identically describe the processes of high-temperature precipitation, which provides the basis for the defect formation in crystals. It is important to note that the Vlasov model of a solid state does not require introduction of additional parameters, such as the growth rate and the axial temperature gradient of the crystal, into the calculations.

Algorithm for Calculating the Growth and
Coalescence of Precipitates In the classical theory of nucleation and growth of newphase particles, the process of precipitation in a crystal is treated as a firstorder phase transition and the kinetics of this process is divided into three stages: the formation of newphase nuclei, the growth of clusters, and the coalescence stage [18]. At the second stage of the precipitation process, clusters grow without a change in their number. This growth is accompanied by a considerable decrease in the degree of supersaturation of the solid solution.
The nucleation centers attach and detach monomers, which is described by the system of equations [9]: where i N is the volume average concentration of nucleation centers that attach i particles; N is the monomer concentration; i k N is the rate of attachment of a monomer for a nucleation center; i g is the rate of detachment of a monomer for a nucleation center. At the initial instant of time, the system contains only monomers and nucleation centers. The growth of precipitates is limited by the monomer diffusion. The kinetic coefficients are given by the formula 4 , where i R is the radius of attachment of a free particle by a cluster consisting of i particles; D is the diffusion coefficient of a free particle.
The system of equations (23) In the case of the diffusion controlled precipitation, the mathematical expectation ( ) i t can be described by the macroscopic equation [19]: where 0 4 i k R D π = , m is the initial size of precipitates; α is the parameter dependent on the cluster geometry. Expressions (23) and (25) allow us to write the differential equation that describes the variation in the monomer concentration during the decomposition of the solid solution [9]: By numerically solving Eq. (26) simultaneously with Eq. (23) we can calculated the average radius of the precipitate at the growth stage: At the third stage of the precipitation process, when the particles of the new phase are sufficiently large, the supersaturation is relatively low, new particles are not formed and the decisive role is played by the coalescence, which is accompanied by the dissolution of small-sized particles and the growth of large-sized particles. The condition providing for changeover to the coalescence stage is the ratio  [20].
In the general case, the size distribution of precipitates has the following form: The average size of precipitates at the stage of the coalescence is proportional to the cube root of time [20]: where D is the diffusion coefficient of impurity atoms; ( ) 0 cr R t is the initial critical radius; σ is the surface tension at the precipitate-solid solution interface; Ω -atomic volume.

Algorithm for Formation of Microvoids and Dislocation Loops
The nucleation of microvoids and dislocation loops begins at specific temperatures (in the temperature range 1423…1173 К), while precipitates are already formed and participate in competing processes of the growth of particular precipitates and dissolution of the other precipitates [1].

Algorithm for the Formation of Microvoids
The model describing the defect dynamics in a crystal involves the kinetics of Frenkel reactions, nucleation of point defects, growth of clusters, and point defect balance [21].
, , Here, concentrations are due to their diffusion (first term), convection (second term), Frenkel reactions (third term), their consumption for the existing clusters (fourth term), and the formation of new clusters (fifth term). The rate of consumption of point defects for the formation of new clusters is negligible and, subsequently, can be ignored [21]. The diffusionlimited growth rates of clusters (for any b and ) formed at the corresponding quantities d and c are described by the equations: where ψ is the density of monomers in the cluster.
The equations for the nucleation rates are represented in the form: where e * is the maximum change in the free energy. The initial length or height of the crystal is taken to be 0. For the crystal with a finite length ℎ , the equilibrium conditions are assumed to be dominant over the entire surface, including the crystal-melt interface: The initial size of the steady state critical cluster is negligible as compared to the microdefect size. Consequently, the initial size of the stable cluster insignificantly affects its final size. The initial size of the stable cluster (nucleus) is calculated from the number of monomers in the critical cluster: where avg stands for the average volume. The model of point defect dynamics under consideration is onedimensional in nature, and, hence, the influence of the radial diffusion that is dominant in the vicinity of the surface of the crystal is disregarded. Therefore, the model can be applied to an axial distribution of defects in the crystal for fixed radial positions far from the surface.
The system of equations is solved using the exact spacetime discretization. The algorithm involves the solution to the An analysis of the experimental and calculated data within microvoid formation in accordance with the diffusion model of the formation of grown-in microdefects revealed the following reasons for the occurrence of microvoids in dislocation-free Si single crystals [22]: (i) a sharp decrease in the concentration of background impurity that was not associated into impurity agglomerates (formed in the cooling range of 1685…1423K); (ii) a large (over 80 mm) crystal diameter (in this case vacancies fail to drain from the central part of the crystal to the lateral surface); (iii) crystals of large diameter generally contain a ring of D-microdefects which forms due to the emergence of the (111) face on the crystallization front and which depletes the region inside with impurity atoms.

Algorithm for the Formation of Dislocation Loops
We calculated the formation of microvoids and dislocation loops according to rigorous approximation for the point defect dynamics model subject to no recombination of intrinsic point defects at high temperatures [23]. It was proved that the process of microvoid formation has a homogeneous nature. However, the formation of dislocation loops is determined, mainly, by the deformation mechanism. This conclusion was made on the basis of a three order of magnitude discrepancy between the experimentally observed concentration of interstitial dislocation loops and their estimated value. [23]. Precipitates originate from elastic interaction between point defects. They are, initially, present in coherent, elastic and deformable state, when lattice distortions close to the precipitate-matrix boundary are not large, and one atom of the precipitate corresponds to one atom of the matrix [16]. Elastic deformations and any mechanical stress connected with them cause a transfer of excessive (deficient) substance from the precipitate or vice versa. Storage of elastic strain energy during the precipitate growth results in a loss of coherence by matrix. [24]. This leads to structural relaxation of precipitates. Structural relaxation of precipitates causes the formation of dislocation loops.
In the volume of silicon the precipitate produces a stress field caused by mismatch between the lattice parameters of precipitate ( ) 1 a and surrounding matrix ( ) 2 a [25]. Then, the intrinsic deformation of the precipitate is defined as described bellow where J is the shear modulus; υ is the Poisson's ratio.
Starting from a certain critical value of the crit R radius, the elastic strain energy mechanism begins to work. This mechanism results in the formation of a prismatic interstitial dislocation loop. The energy criterion for such mechanism is the initial final constitute elastic energy of the system with precipitate before and after relaxation.
In respect of a spherical precipitate with equiaxial intrinsic deformation, the calculation of elastic fields of the precipitate is substantially simplified. Let us assume that the intrinsic elastic strain energy of the precipitate before and after the formation of a dislocation loop of mismatch remains constant where d is loop diameter; f is the radius of the core loop; b is the magnitude of the Burgers vector. A critical value of the precipitate radius corresponds to a value at which the loop is formed on the precipitate [27] ( ) where α is a constant contribution of the dislocation core.
Formula (52) is approximate and can be used only to determine the value of crit R . critical radius. The dislocation loops with a radius of crit R R > become bigger in size at the coalescence stage, while small dislocation loops with a radius of crit R R < will dissolve [28,29]. The growth of dislocation loop in course of consequent as grown silicon crystal' cooling occurs both due to the dissolution of small loops with the sizes less than critical, and the oversaturation of silicon self-interstitials. In this case, the crystal growth ratio is g where ( ) D t is the diffusion coefficient of intrinsic interstitial silicon atoms; t is the time cooling the crystal; j is the proportionality factor. The crystal cooling time value is defined from the dependence: is the cooling rate of the crystal.
Let us assume that the formation of dislocation loops is defined only by deformation mechanism. Then, a concentration of dislocation loops during the crystal cooling shall be considered as the concentration function of precipitates. We have shown earlier that at the stage of growth and coalescence of precipitates their concentration is the function of the crystal cooling time [18]. Then, the loop concentration depends on the crystal cooling time [29]: relaxation by precipitate results in occurrence of one or more dislocation loops [25]. The critical radiuses of precipitates which are energy favorable to the nucleation of dislocation loop shall be crit R = 3,028 µm (SiO 2 ) and crit R = 3,402 µm (SiC).
The process of interaction of intrinsic point defects with impurities begins near the crystallization front. It is decisive in the formation of the defect structure of highly perfect dislocation-free single crystals of silicon during their growth. During the cooling of the crystal, depending on the thermal growth conditions, the impurity atoms die to form precipitates. Conditions are created for exceeding the equilibrium values of the concentrations of intrinsic point defects. This process leads to the formation of microvoids or interstitial dislocation loops in different regions of the crystal.

Unified Algorithm for the Formation of the Original Defect Structure
A unified algorithm for the formation of defect structure is based on the diffusion model for the formation and transformation of grown-in microdefects. The diffusion model fully describes the kinetics of the diffusion decay of supersaturated solid solutions of point defects during the cooling of the crystal after growing. All parameters of precipitates, microvoids and dislocation loops are determined through thermal parameters of crystal growth. Determination of the type of defects and calculation of the formation of grown-in microdefects is performed depending on the value of crystal growth rate, temperature gradients and cooling rate of the crystal [3].
Using a parabolic distribution of the axial temperature gradient, it is possible to calculate the defect structure at any point of the crystal. The parabolic radial distribution of the axial temperature gradient has the form where h is the axial temperature gradient at the crystal edge; j k is the crystal radius and / is the current coordinate in the range of 0 … j k . The dependence of the critical growth rate on the crystal radius is plotted proceeding from the relation to the so-called oxidation-induced stacking fault (OSF) ring, which is observed in crystals after thermal treatments [1].
The algorithm of its application is based on the hightemperature precipitation process, which determines the subsequent course of the formation of a defect structure of crystals ( Figure 1). Experimental database (EDB) should contain all available experimental data on the crystal, for example, (i) data on the crystal structure, growth method, and thermal parameters of the crystal growth; (ii) data on the types of defects, their size and concentration, and temperatures of formation of defects; (iii) data on the physical model of defect formation; etc. The necessary data from EDB are supplied to input module of the algorithm for the calculation of the defect structure in high-temperature precipitation model.
Computation of precipitates is carried out within the framework of probabilistic and classical theories of nucleation, growth and coalescence of precipitates using analytical and approximate calculations. Parameters characterizing the processes of precipitation of oxygen and carbon are determined. Among them, for example, are the critical radii of precipitates, the size distribution of precipitates, the change in the average size of the precipitates during crystal cooling, and others. Accounting for other impurities is carried out in the form of carbon and oxygen components of the calculation. For example, the sum of the concentrations of impurities of carbon, boron, aluminum and hydrogen form the carbon component of the calculation. Similarly, the sum of the concentrations of impurities of oxygen, gallium, phosphorus, antimony, arsenic, germanium, nitrogen and iron form the oxygen component of the calculation.
Next are the calculation and analysis of the formation of the secondary defect structure of crystals (for example, dislocation loops, microvoids, etc.). The results of calculations are compared with the EDB data. In the case of good agreement with the EDB data, the results of calculations can serve as the characteristics of the initial defect structure of crystal and used for the calculation and analysis of the defect structure of the crystal after different heat treatments. If a satisfactory agreement is not achieved, the calculations are repeated again until the results are obtained with a desired accuracy. In the process, the mathematical models and (or) experimental data from EDB can be refined.
In the calculation of microvoids, a check is first made to determine the conditions for their formation. так как This is due to the fact that microvoids are not formed at a cooling rate of the crystal ≥ 40 K/min [31] and in crystals with a diameter less than 70 mm [22]. Calculation of microvoids and dislocation loops makes it possible to determine for each of these types of defects such parameters as the critical radius, concentrations and sizes.
The algorithm makes it possible to develop software for the analysis and calculation of the defect structure of silicon. The output of information can be done in tabular, graphical and animation forms.

Conclusion
Mathematical models of the formation of grown-in microdefects gave the equivalent of the object (the process of defect formation during crystal growth), which in mathematical form reflects its most important properties. Mathematical models are presented in a form convenient for the application of numerical methods. The computational algorithm does not distort the basic properties of the model and, therefore, of the original object. The adequacy of the developed package "mathematical model -algorithm" of the original physical model is verified by carrying out computational experiments and their comparison with the results of physical studies.
The advantage of calculating the formation and growth of precipitates, microvoids and dislocation loops is that the input and control parameters of the computational and graphical part are the growth parameters: the rate of crystal growth and the axial temperature gradient. The axial temperature gradient is chosen in a certain range of values depending on the diameter of the grown crystal. Of the three growth parameters, the two parameters are strictly defined (the diameter of the crystal and the rate of its growth), and the axial temperature gradient is specified in a certain range of values.
In this paper, we proposed an algorithm for the calculation and analysis of a defect structure of crystals and presented the schematic diagram of this algorithm. Based on the schematic diagram, it is possible to develop a software product for analysis and calculation of a defect crystal structure.