Evaluation of the Performance of Published Point Defect Parameter Sets in Cone and Body Phase of a 300 mm Czochralski Silicon Crystal

Abstract

Prediction and adjustment of point defect (vacancies and self-interstitials) distribution in silicon crystals is of utmost importance for microelectronic applications. The simulation of growth processes is widely applied for process development and quite a few different sets of point defect parameters have been proposed. In this paper the transient temperature, thermal stress and point defect distributions are simulated for 300 mm Czochralski growth of the whole crystal including cone and cylindrical growth phases. Simulations with 12 different published point defect parameter sets are compared to the experimentally measured interstitial–vacancy boundary. The results are evaluated for standard and adjusted parameter sets and generally the best agreement in the whole crystal is found for models considering the effect of thermal stress on the equilibrium point defect concentration.

1. Introduction

The Czochralski (Cz) process is the method of choice for the production of silicon (Si) crystals for microelectronics applications. As the feature size of electronic devices is steadily decreasing, reduction of the size and number of grown-in defects as voids and self-interstitial agglomerates in Si wafers is of utmost importance for the process development. Accordingly, the distribution of intrinsic point defects (PD: vacancies, V, and self-interstitials, I)—the source of possible larger defect agglomerates—must be predicted and controlled very precisely during the crystal growth process.

The PD concentration depends on the thermal gradient and pulling velocity. The prevalent point defect type can be simply predicted by the $v/G$ criterion as shown by Voronkov [1], where v is the crystal growth rate (equal to the pull rate V under steady-state conditions) and G is the temperature gradient in the crystal at the crystallization interface. The 2D axisymmetric numerical simulation of vacancy and interstitial concentrations considering incorporation, diffusion, convection and recombination validates the concept for different process parameters [2,3] and opens the way for process optimization regarding the defect distribution and reduction [4,5].

The global thermal modeling of the Cz process was started in the 1990s [6,7,8]. Quasi-steady and transient 2D models considered radiation heat transfer between heater, insulation, crucible and Si melt and crystal, as well as conduction in solid bodies and melt, allowing for the optimization of hot-zone and pulling velocity. Later simulation of the point defect distribution in the crystal domain was added [2,3,4,9,10]. Transient modeling of both the heat transfer and PD distribution is necessary since the Cz process is transient in reality, and the start and end cones are a considerable part of the whole growth process.

In the next step of PD model development [11], the effect of thermal stress on the point defect distribution was included. While the stress effect was first theoretically examined for a 150 mm crystal [12,13] and found out to be weak, later experimental evidence for a larger crystal with a diameter of 300 mm [14] demonstrated its importance for industrial-size crystals. This topic received greater attention only in the 2010s by the works of Vanhellemont [15,16,17]. Stress changes PD formation enthalpy, affecting the equilibrium concentration and the critical value of $v/G$, which has been confirmed experimentally [18] and explained theoretically [19]. Since the stress generally increases with the crystal size, this effect is more important for large crystals and is actively investigated nowadays [20,21,22,23].

While the PD parameter sets proposed by different groups are validated and fine-tuned to improve agreement with experiments in the body phase [3,4,20,23,24,25,26] (we are not aware of works comparing simulation results with experiment in the cone), they vary considerably. One of the possible reasons for such discrepancies is uncertainties in the global heat transfer (e.g., due to imprecise material properties) since the point defect distribution is sensitive to the temperature field in the crystal. Thus, two PD parameter sets cannot be simply compared if they were derived based on different heat transfer models. A detailed comparison of the performance of published PD parameter sets for the same experiment is not available in the literature.

In the present work, the software CZ-Trans is used for dynamic simulations of the temperature field, phase boundaries [27] and PD distribution in the crystal [28]. The novelty is that we (i) apply an improved and numerically more precise CZ-Trans model; (ii) simulate 12 different PD parameter sets using the same thermal conditions; (iii) compare simulation results with the experimental PD distribution in the whole crystal, including start and end cones; (iv) adjust PD parameters for 1 position in the crystal and find the parameter sets with the best agreement to experiment.

2. Growth Experiment

A 300 mm (12 inch) Cz silicon single crystal growth experiment carried out at Siltronic AG was selected, see Figure 1. The pull rate was slowly changed during the growth, which resulted in different radial positions of the I-V boundary (IVB) along the crystal length. The minority carrier lifetime in the cone, body (cylinder) and endcone was measured by the microwave photoconductivity decay (µ-PCD) method, which revealed the IVB position [3,29], while the lateral photovoltage scanning (LPS) measurements (not shown) were used to investigate the crystallization interface shape.

The cone of the grown crystal was completely interstitial-rich but starting from the shouldering stage (transition from cone to cylindrical growth), defined as the crystal length $L=0$, a vacancy-rich region developed near the crystal axis and was present throughout the whole body phase. A more complicated structure developed near the endcone and the vacancy region in the center reached the crystal surface around $L=1200$mm.

What is also important is the crystallization interface shape, which was extracted from the LPS measurements or, for crystal parts where they were unavailable, from µ-PCD data. The measured crystallization interface deflection $H_C$ as a function of crystal length is plotted in Figure 2, while the interface shapes are depicted in Figure 3. Since the melt flow was not directly modeled, this information was used to improve the agreement between simulations and experiment as the interface shape affects both the thermal stress and point defect distribution.

3. Numerical Model

3.1. Heat Transfer and Phase Boundaries

The program CZ-Trans  [27] was used to model the transient heat transfer in the Cz system and the changing phase boundaries. The model was axisymmetric, and the temperature field $T(r,z)$ was calculated in the crystal, melt, crucible and heat shield domains while a simple integral model was used for the heater to calculate its time-dependent temperature as a function of the heater power change. Radiative heat transfer between surfaces was considered. Since the crystal was pulled upwards, both the convection and diffusion terms were included in the temperature equation

$$\frac{\partial T}{\partial t}+\mathbf{V}\mathbf{·}\mathbf{\nabla }T=\mathbf{\nabla }\mathbf{·}\left({\mathit{\lambda }}_{\mathrm{c}}T\right),$$

(1)

where ${\mathit{\lambda }}_{\mathrm{c}}$ is the thermal conductivity and $\mathbf{V}$ is the crystal pull rate, which has only a vertical component ${\mathbf{V}}_z=V$.

While the melt flow was not simulated, to correctly capture the crystallization interface shapes its influence was approximated by (i) homogeneous but time-dependent melt thermal conductivity ${\mathit{\lambda }}_{\mathrm{m}}\left(t\right)$ and (ii) redistribution of the heat flux densities along the crystallization interface [30]. For the crystal ${\mathit{\lambda }}_{\mathrm{c}}\left(T\right)=98.89-9.43\times {10}^{-2}T+2.89\times {10}^{-5}{T}^2$ (in $\mathrm{W}/\mathrm{m}/\mathrm{K}$) was used; emissivities were 0.64 and 0.30 for the solid and liquid Si, respectively. The melting point of silicon $T_0=1685K$.

For a more precise numerical modeling of the heat transfer (especially at the crystallization interface, since it affects the crystal growth rate), the possibility of using high-order finite elements was implemented in CZ-Trans. The developed temperature solver uses the open-source finite element library deal.II [31]. While the solver supports arbitrary element orders, from a practical point of view a good balance between the accuracy and performance is achieved for second-order elements (9-node quadrilaterals).

Since a quad-only mesh is required for deal.II, it was introduced for the crystal, melt and crucible domains. Where possible, a structured grid was generated by CZ-Trans, otherwise an unstructured mesh was created using Gmsh [32].

Newton’s method with a step length of 1 was employed to deal with the non-linearities, such as the temperature-dependent thermal conductivity. The obtained temperature field on second-order finite elements in the crystal was used to calculate the thermal stress and point defect distribution on the same mesh.

3.2. Point Defects

3.2.1. Governing Equations

The conservation equations for point defects can be written in a general form as

$$\frac{\partial C}{\partial t}+\mathbf{V}\mathbf{·}\mathbf{\nabla }C=-\mathbf{\nabla }\mathbf{·}\mathbf{J}+k^{\mathrm{r}}\left(C_I^{\mathrm{eq}}{C}_V^{\mathrm{eq}}-C_I{C}_V\right),$$

(2)

with: $C=C_I$, $C_V$—point defect concentrations, $C^{\mathrm{eq}}$—equilibrium point defect concentration, D—diffusion coefficient, $k^{\mathrm{r}}$—reaction constant for recombination. The zero isoline of $\Delta =C_I-C_V$ is the I-V boundary (IVB) [1]. For convenience, we introduce two non-negative variables ${\Delta }_I=max(\Delta ,0)$ and ${\Delta }_V=max(-\Delta ,0)$ for interstitial-rich and vacancy-rich regions, respectively.

The diffusive flux of the point defects, $\mathbf{J}$, is given by

$$\mathbf{J}=-D\mathbf{\nabla }C+\frac{Q^{\ast }D}{k{T}^2}C\mathbf{\nabla }T,$$

(3)

where $Q^{\ast }$ is the activation enthalpy for thermal diffusion and $k=8.617\times {10}^{-5}\mathrm{e}\mathrm{V}/\mathrm{K}$ is the Boltzmann constant. $Q^{\ast }>0$ means an uphill PD drift—in the direction of ∇T, from a colder to a hotter crystal part since the drift velocity is given by ${\mathbf{v}}_{\mathrm{drift}}=-D\mathbf{\nabla }(Q^{\ast }/\left(kT\right))=D{Q}^{\ast }/\left(k{T}^2\right)\mathbf{\nabla }T$ [33]. Note that $Q^{\ast }<0$ is also a valid value. The definition of the thermal diffusion term (its sign) is not consistent in the literature, therefore some PD parameter sets [18,34,35,36] were converted to match the formulation used in this study.

The temperature dependences of the point defect properties are given by

$$\begin{array}{cc}\hfill D& =D^{0}exp(-\frac{H^{\mathrm{m}}}{kT}),\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill C^{\mathrm{eq}}& =C^{0}exp(-\frac{H^{\mathrm{f}}+\Delta H^{\mathrm{f}}}{kT}),\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill k^{\mathrm{r}}& =4\pi a^{\mathrm{r}}\left(D_I+D_V\right)exp(-\frac{\Delta G^{\mathrm{r}}}{kT}),\hfill \end{array}$$

(6)

where $D^0$ and $C^0$ are prefactors; $H^{\mathrm{m}}$ is the migration energy, $H^{\mathrm{f}}$ and $\Delta H^{\mathrm{f}}$ are formation enthalpy and its change due to the thermal stress, respectively; $a^{\mathrm{r}}$ is the capture radius and $\Delta G^{\mathrm{r}}$ is the energy barrier for recombination. The formation entropy is included in $C^0$.

The critical value of the Voronkov parameter $\xi =v/G$ corresponding to the present formulation, Equations (2) and (3), is [18,33,37]:

$${\xi }_{\mathrm{crit}}=\frac{C_I^{\mathrm{eq}}{D}_I(H_{\mathrm{av}}^{\mathrm{f}}-Q_I^{\ast })-C_V^{\mathrm{eq}}{D}_V(H_{\mathrm{av}}^{\mathrm{f}}-Q_V^{\ast })}{\left(C_V-C_I\right)k{T}_0^2},H_{\mathrm{av}}^{\mathrm{f}}=\frac{H_I^{\mathrm{f}}+H_V^{\mathrm{f}}}{2},$$

(7)

where all quantities are evaluated at the melting point (MP) of Si, $C_V=C_V^{\mathrm{eq}}$ and $C_I=C_I^{\mathrm{eq}}$. The crystal is interstitial-rich for $\xi <{\xi }_{\mathrm{crit}}$ and vacancy-rich for $\xi >{\xi }_{\mathrm{crit}}$. The radial position of the crystallization interface where $\xi ={\xi }_{\mathrm{crit}}$ corresponds to the IVB, which separates both point defect regions. Note that the Voronkov criterion is an approximation; for precise results Equation (2) should be solved numerically.

3.2.2. Numerical Aspects

The coupled point defect equations for $C_I$ and $C_V$ (2) are solved using second-order finite element method (FEM) on a quadrilateral mesh. Newton’s method is employed to deal with the non-linear recombination term [38]. The system of linear equations is solved using a direct method, obtaining the Newton updates $\delta C_I$ and $\delta C_V$. Two Newton iterations are typically sufficient for one time-step of a transient simulation.

Unless the simulation is steady-state, the pull rate V in Equation (2) is set to zero and the convection is taken into account by moving the crystal mesh upwards by $V\Delta t$. The fields are interpolated if the mesh is regenerated. The same procedure is also applied for the transient temperature field calculations, Equation (1) [27,39].

At temperatures below 1000 ${}^{\circ }\mathrm{C}$, the PD diffusion coefficients are so low that the defects are practically frozen-in. To avoid numerical diffusion (e.g., due to remeshing and interpolation) during unsteady modeling of defect transport in the crystal, the approach proposed by the STR Group and available in the commercial software package CGSim [40] has been used: a static (unchanged) grid is created in the low-temperature part of the crystal, which is automatically extended during the simulation as the crystal grows.

3.3. Considered Parameter Sets

Modeling was performed for twelve different PD parameter sets available in the literature, which are listed in

Table 1. Summary of considered point defect parameter sets defined by Equations (3)–(5). Suffix “-s” denotes parameter sets with the thermal stress influence, see Table 3 for the relevant parameters.

Parameter Set		Ref.		${\mathit{D}}^0({\mathbf{cm}}^2/\mathit{s})$	${\mathit{H}}^{\mathbf{m}}\left(\mathbf{eV}\right)$	${\mathit{C}}^0\left({\mathbf{cm}}^{-3}\right)$	${\mathit{H}}^{\mathbf{f}}\left(\mathbf{eV}\right)$	${\mathit{Q}}^{\ast }\left(\mathbf{eV}\right)$
A	2000-Nakamura	[41]	I	1.040 × 104	2.4	1.060 × 1022	2.4
			V	2.140	1.4	5.290 × 1022	2.6
B	2001-Wang	[34,35]	I	2.101 × 10−1	0.907	3.945 × 1026	3.943	−3.66
			V	1.000 × 10−4	0.489	2.675 × 1023	2.848	−2.66
C	2002-Nakamura	[36]	I	2.459 × 10−1	0.9	6.284 × 1026	4.05	−1.01
			V	3.513 × 10−4	0.3	3.951 × 1026	3.94	0.0
D	2007-Kulkarni-I	[25,42]	I	1.950 × 10−1	0.9	6.176 × 1026	4
			V	6.262 × 10−4	0.4	7.520 × 1026	4
E	2007-Kulkarni-II	[25,42]	I	4.000 × 10−3	0.3	4.725 × 1027	4.3492
			V	2.000 × 10−3	0.38	1.200 × 1027	4.12
F	2007-Sinno	[43]	I	2.370 × 10−1	0.937	6.365 × 1026	4.0
			V	7.870 × 10−4	0.457	9.931 × 1025	3.7
G	2009-Voronkov	[33]	I	3.667 × 10−3	0.2	1.884 × 109	4.95	4.5
			V	1.876 × 10−3	0.38	2.967 × 1026	3.95	29
H	2011-Nishimoto	[26,44]	I	2.590	1.18	5.992 × 1025	3.77
			V	9.918 × 10−4	0.4	1.400 × 1026	3.84
I	2013-Vanhellemont-s	[17]	I	3.800 × 10−2	0.88	6.400 × 1025	3.68
			V	1.200 × 10−3	0.45	2.580 × 1026	3.88
J	2014-Nakamura-s	[18]	I	2.459 × 10−1	0.9	6.284 × 1026	4.05	−1.01
			V	3.513 × 10−4	0.3	3.951 × 1026	3.94	0.0
K	2016-Kamiyama-s	[20,26]	I	2.590	1.18	5.992 × 1025	3.77
			V	9.918 × 10−4	0.4	1.400 × 1026	3.84
L	2019-Mukaiyama-s	[21,22]	I	2.459 × 10−1	0.9	6.284 × 1026	4.05
			V	3.442 × 10−4	0.3	4.030 × 1026	3.94

. The parameter values at the melting point of Si are given in

Table 2. Values of point defect parameters at the melting point.

Set	$\mathit{D}({10}^{-4}{\mathbf{cm}}^2/\mathit{s})$		${\mathit{C}}^{\mathbf{eq}}\left({10}^{14}{\mathbf{cm}}^{-3}\right)$		${\mathit{C}}_{\mathit{V}}^{\mathbf{eq}}-{\mathit{C}}_{\mathit{I}}^{\mathbf{eq}}$	${\mathbf{DC}}^{\mathbf{eq}}\left({10}^{10}{\mathbf{cm}}^{-1}{\mathit{s}}^{-1}\right)$
	I	V	I	V	$\left({10}^{14}{\mathbf{cm}}^{-3}\right)$	I	V
A	6.898	1.390	7.030	8.850	1.819	48.495	12.303
B	4.070	0.034	6.348	8.111	1.762	25.839	0.280
C	5.000	0.445	4.840	6.490	1.650	24.200	2.888
D	3.964	0.398	6.712	8.172	1.461	26.603	3.256
E	5.067	1.460	4.635	5.707	1.071	23.489	8.334
F	3.734	0.338	6.917	8.520	1.602	25.831	2.881
G	9.250	1.370	2.950	4.550	1.600	27.288	6.234
H	7.656	0.631	3.174	4.580	1.406	24.300	2.890
I	0.887	0.541	6.301	6.407	0.106	5.587	3.467
J	5.000	0.445	4.840	6.490	1.650	24.200	2.888
K	7.656	0.631	3.174	4.580	1.406	24.300	2.890
L	5.000	0.436	4.840	6.620	1.780	24.200	2.886

, where it can be seen that the variations in D and $C^{\mathrm{eq}}$ between different parameter sets are usually stronger than variations in the transport capacity $P=D{C}^{\mathrm{eq}}$, which can be measured experimentally and is typically used to extract the temperature-dependent D and $C^{\mathrm{eq}}$ [4]. Two outliers are set B, which has low $D_V$ and $P_V$, and set I, which has low $D_I$, $P_I$ and $C_V^{\mathrm{eq}}-C_I^{\mathrm{eq}}$.

The last four sets include the influence of the thermal stress as described in the next subsection and summarized in

Table 3. Thermoelastic parameters (Young’s modulus E, coefficient of thermal expansion $\alpha$, Poisson’s ratio $\nu$) at the MP and parameters for the thermal stress influence on the point defects $\mathrm{d}{H}^{\mathrm{f}}/\mathrm{d}\sigma$.

Parameter Set		E	$\mathit{\alpha }$	$\mathit{\nu }$	$\mathbf{d}{\mathit{H}}^{\mathbf{f}}\mathbf{d}\mathit{\sigma }\left(\mathbf{eV}/\mathbf{GPa}\right)$
		$\left(\mathbf{GPa}\right)$	$\left({10}^{-6}{\mathit{K}}^{-1}\right)$	$(-)$	I	V
I	2013-Vanhellemont-s	100	4.66	0.25	−0.070	0.160
J	2014-Nakamura-s	139.7	4.66	0.225	−0.070	0.154
K	2016-Kamiyama-s	166	4.66	0.217	−0.0694	0.3069
L	2019-Mukaiyama-s	165	4.66	0.217	−0.0694	0.3069

. In particular, it is the only difference between sets C and J, as well as between H and K.

The so-called fast recombination [4,33], meaning that $k^{\mathrm{r}}$ is sufficiently large to ensure $C_I{C}_V\approx C_I^{\mathrm{eq}}{C}_V^{\mathrm{eq}}$ near the crystallization interface, was applied for all cases by setting $a^{\mathrm{r}}=1$nm and $\Delta G^{\mathrm{r}}=1.5$eV. The first-type (Dirichlet) boundary condition $C=C^{\mathrm{eq}}$ was applied at the crystallization interface and crystal side surface. On the symmetry axis $r=0$, the no-flux boundary condition (BC) $\partial C/\partial n=0$ was used. The influence of $\Delta G^{\mathrm{r}}$ and BC at the crystal surface was checked for some parameter sets and is discussed in Section 4.3.

3.4. Thermal Stresses

At each time step the axisymmetric thermal stresses ${\sigma }_{rr}$, ${\sigma }_{zz}$, ${\sigma }_{\varphi \varphi }$ and ${\sigma }_{rz}$ in the crystal are calculated using second-order finite elements on a quadrilateral mesh. Since the largest temperature gradients and highest stresses are near the crystallization interface, the thermoelastic parameters at the melting point are used, neglecting their temperature-dependence. A summary of the relevant parameters is given in Table 3. Note that the last two thermoelastic parameter sets are almost identical but the PD parameters differ as given in Table 1.

The mean thermal stress ${\sigma }_{\mathrm{ave}}$ affects the equilibrium PD concentration by changing the formation enthalpy that is used by the point defect solver:

$$\Delta H^{\mathrm{f}}=\frac{\mathrm{d}{H}^{\mathrm{f}}}{\mathrm{d}\sigma }{\sigma }_{\mathrm{ave}},{\sigma }_{\mathrm{ave}}=\frac{{\sigma }_{rr}+{\sigma }_{zz}+{\sigma }_{\varphi \varphi }}{3}.$$

(8)

Typically, ${\sigma }_{\mathrm{ave}}<0$ (compressive stress) near the axis, which, according to the coefficients $\mathrm{d}{H}^{\mathrm{f}}/\mathrm{d}\sigma$ in Table 3, increases $C_V$ and shifts the IVB to the outer part of the crystal, compared to the case without the stress influence. The reason is that under compressive stress the distance between atoms in a crystal lattice is smaller and the stress is relaxed by the increase of $C_V$ or decrease of $C_I$.

Depending on whether the stress is isotropic or planar, the theoretically predicted $\mathrm{d}{H}^{\mathrm{f}}/\mathrm{d}\sigma$ can vary significantly [19,22]. The latter parameter set with larger $\mathrm{d}{H}_V^{\mathrm{f}}/\mathrm{d}\sigma$ is more realistic, as confirmed experimentally by analyzing the relationship ${\xi }_{\mathrm{crit}}\left({\sigma }_{\mathrm{ave}}\right)$  [18].

4. Results and Discussion

4.1. Heat Transfer and Phase Boundaries

The transient CZ-Trans simulations were started from the beginning of the cone phase. PID control of the crystal pull rate was applied according to the difference between the actual crystal radius R and its setpoint, which was maintained within a few millimetres. Figure 1 demonstrates a good agreement between the modeled pull rate curve and experimental data. One reason for discrepancies can be due to long-scale temperature fluctuations caused by turbulent melt flow. While the pull rate oscillations are much stronger in the experiment during the cone phase, it practically does not influence the heat transfer, as verified by performing a simulation with non-smoothed experimental data as the pull rate setpoint; this also holds for the point defect distribution.

An important result is the crystallization interface shape, which evolves in time according to the growth rate v which, in its turn, depends on the heat flux densities $q_{\mathrm{c}}$ and $q_{\mathrm{m}}$ in the crystal and melt, respectively. While the heat transfer in the crystal and $q_{\mathrm{c}}$ is straightforward to calculate, $q_{\mathrm{m}}$ can be strongly influenced by the melt flow, which is not solved explicitly. Its influence is included in the model by introducing the heat flux density correction $q_{\mathrm{corr}}$ [30], which describes redistribution of heat sources along the interface, calculated based on the experimental interface shapes. The values of $q_{\mathrm{corr}}\left(r\right)$ are obtained at certain positions in the crystal and are linearly interpolated according to the time-dependent crystal length. In order not to tamper with the global heat transfer, the integral correction should be small; in the present results it does not exceed 0.2 kW while the typical heat flux from the melt is 5 kW.

To verify the accuracy of the phase boundary calculations, the crystallization interface deflection was compared with experiment, obtaining a rather good agreement within a few millimeters, see Figure 2. The non-smooth behavior at $L\approx 1100mm$ is likely to be caused by the change in pull rate (Figure 1).

Also shown in Figure 2 is the $v/G$ ratio at $r=0$ and $r=R$. Assuming the typical value of ${\xi }_{\mathrm{crit}}=1.3\times {10}^{-3}\mathrm{c}{\mathrm{m}}^2{\min }^{-1}{\mathrm{K}}^{-1}$, the crystal outer surface is expected to always be interstitial-rich. On the axis, a transition from the interstitial-rich beginning of the cone to the fully vacancy-rich inner crystal region is predicted at $L\approx -95$mm. Note that these results were obtained without the thermal stress influence and the Voronkov criterion does not consider the interface shape and might be inaccurate for strongly transient growth (cone). The full transient PD simulations are presented in the following subsection.

Besides the interface deflection, the whole interface shape $z\left(r\right)$ should be correct as it affects the temperature gradient, thermal stress and PD distribution. The comparison to the experimental data plotted in Figure 3 confirms the good agreement. Since the calculated interfaces are plotted at regular intervals, one can visually estimate v and, e.g., verify that v at the beginning of the body phase is higher than at the end. Although the difference in v between the crystal center and TP is harder to see with the naked eye, it is noticeable during shouldering ($L=0$).

4.2. Thermal Stresses

The thermal stresses depend strongly on the interface shape; in case of a concave interface typical for the body phase (Figure 4) the stress at the interface is compressive (${\sigma }_{\mathrm{ave}}<0$) near the axis and tensile close to the crystal surface. Relevant for the point defect simulations is the change of PD formation enthalpy $\Delta H_I^{\mathrm{f}}>0$ and $\Delta H_V^{\mathrm{f}}<0$ according to the parameters in Table 3. The stress effect leads to the increase of $C_V^{\mathrm{eq}}$ and decrease of $C_I^{\mathrm{eq}}$ near the center. Further away from the interface the magnitude of ${\sigma }_{\mathrm{ave}}$ becomes lower and the compressive/tensile regions switch places.

As seen in Figure 4, the described tendencies are opposite during the cone growth, thus the IVB will be shifted towards the axis. Towards the end of the body phase the interface stress becomes entirely compressive but in the endcone reverts to the standard body phase behavior.

The difference in thermoelastic parameters between parameter sets is primarily due to the Young’s modulus, E, therefore by estimating ${\sigma }_{\mathrm{ave}}\sim E$ the stresses for sets K and L are 65–66% higher than for set I, see Figure 4. Regarding the PD, $\mathrm{d}{H}_I^{\mathrm{f}}\mathrm{d}\sigma$ is roughly the same for all sets but $\mathrm{d}{H}_V^{\mathrm{f}}\mathrm{d}\sigma$ for sets K and L is almost twice as large than for I and J, therefore the formation enthalpy change $\Delta H_V^{\mathrm{f}}$ for sets K and L is about 3.2 times as high as for set I.

4.3. Point Defects

4.3.1. Overview

As mentioned previously, the temperature, stress and point defect fields in the crystal were solved on a quadrilateral mesh using second-order FEM. Before all parameter sets were investigated, the influence of numerical parameters was checked for parameter set F. Two additional calculations were carried out—one with the time step reduced about 2 times and another with a crystal mesh that was around twice as fine. In both cases no noticeable influence on the IVB position was observed.

The energy barrier for recombination $\Delta G^{\mathrm{r}}$ was changed from its default value of 1.5 eV to 2 eV and 1 eV but the difference was small, especially for the latter case. This corresponds to the generally accepted idea that in the fast recombination regime $k^{\mathrm{r}}$ should be sufficiently high but its exact value is unimportant.

The influence of the boundary condition at the crystal side surface was checked as well. When using the no-flux BC $\partial C/\partial n=0$ instead of the default $C=C^{\mathrm{eq}}$, the crystal outer surface became more interstitial-rich. The maximum of $\Delta$, which was previously located a few centimeters away from the crystal surface, shifted to the surface. However, the BC affected the results only at distances up to 5 cm from the crystal surface while the experimental IVB was almost everywhere further inside (Figure 1). In the cone $C_I$ increased but the axial position at which $\Delta =0$ remained unchanged.

$C(t=0)=C^{\mathrm{eq}}$ was used as the initial condition. Due to the small size of the seed crystal it did not influence the final PD distribution, which was confirmed by applying the calculated steady-state concentration field at the initial time.

Two series of the PD simulations are carried out in the following subsections. First, the original parameter sets given in Table 1 are considered in the calculations and the obtained IVB compared to the experimental data. Since the values of ${\xi }_{\mathrm{crit}}$ differ between each parameter set, see

Table 4. Adjustment of $C_I^{\mathrm{eq}}$ and the corresponding values of ${\xi }_{\mathrm{crit}}$ (in ${10}^{-3}\mathrm{c}{\mathrm{m}}^2{\min }^{-1}{\mathrm{K}}^{-1}$) at ${\sigma }_{\mathrm{ave}}=-20$MPa for original and adjusted PD parameters.

Parameter Set		$\Delta {\mathit{S}}_{\mathit{I}}^{\mathbf{f}}/\mathit{k}$	${\mathit{C}}_{\mathit{I}}^{\mathbf{eq}}$ Scale	${\mathit{\xi }}_{\mathbf{crit}}$ orig.	${\mathit{\xi }}_{\mathbf{crit}}$ adj.
A	2000-Nakamura	0.037	1.038	1.22	1.50
B	2001-Wang	−0.1802	0.835	2.51	1.31
C	2002-Nakamura	−0.056	0.946	1.63	1.32
D	2007-Kulkarni-I	−0.0322	0.968	1.57	1.32
E	2007-Kulkarni-II	−0.031	0.969	1.47	1.24
F	2007-Sinno	−0.0041	0.996	1.35	1.32
G	2009-Voronkov	−0.6635	0.515	2.32	1.23
H	2011-Nishimoto	−0.014	0.986	1.42	1.36
I	2013-Vanhellemont-s	0.0134	1.013	0.60	0.85
J	2014-Nakamura-s	−0.0306	0.970	1.45	1.29
K	2016-Kamiyama-s	0.0356	1.036	1.20	1.35
L	2019-Mukaiyama-s	0.0795	1.083	0.97	1.31

, they were manually adjusted to compensate for the different thermal models and ensure the correct IVB position in the middle of the body phase ($L=600$mm). Finally, the adjusted PD parameter sets are considered in the full transient simulations and the obtained results are once again compared to the experimental data.

4.3.2. Original PD Parameters

The results with the original PD parameters summarized in Table 1 and Table 3 are depicted in Figure 5. The temperature of 1000 ${}^{\circ }\mathrm{C}$ below which the defects are frozen-in is located in the endcone, meaning that almost the whole crystal can be compared with the experiment.

Only a few PD parameter sets—F, H, I and K—produce reasonable results but none of them is perfect. For sets A and L the crystal center is too vacancy-rich and for the remaining sets too interstitial-rich. The prediction of the IVB during shouldering is especially difficult—its modeled axial position at $r=0$ is typically too high.

These results qualitatively correspond to the theoretical ${\xi }_{\mathrm{crit}}$ in Table 4 and modeled $v/G$ curve in Figure 2. For example, neglecting the stress influence, a fully vacancy-rich center in the body phase is expected for parameter sets with ${\xi }_{\mathrm{crit}}<1.8\times {10}^{-3}\mathrm{c}{\mathrm{m}}^2{\min }^{-1}{\mathrm{K}}^{-1}$. The actual PD results suggest this value is around $1.45\times {10}^{-3}\mathrm{c}{\mathrm{m}}^2{\min }^{-1}{\mathrm{K}}^{-1}$ (between sets E and H). Assuming a typical stress level of $-20\mathrm{M}\mathrm{P}\mathrm{a}$ (sets K and L), the same tendencies are also observed for sets with the stress influence.

The thermal stress influence can be seen by comparing parameter set pairs C & J, and H & K, which have otherwise identical PD parameters: the IVB is shifted away from the axis and becomes straighter (i.e., having weaker radial variations), improving agreement with experiment.

4.3.3. Adjusted PD Parameters

To allow for a more direct comparison of all parameter sets and to improve agreement with experiment, all PD parameter sets were adjusted using quasi-stationary PD and thermal stress simulations on the temperature field obtained in the dynamic simulations. The equilibrium concentration of self-interstitials was adjusted by changing the formation entropy by $\Delta S_I^{\mathrm{f}}$, to match the experimentally observed I-V boundary at $L=600$mm. The corresponding scale factor for $C_I^{\mathrm{eq}}$ in Equation (5) is

$$C_I^{\mathrm{eq}}\mathrm{scale}=exp(\Delta S_I^{\mathrm{f}}/k).$$

(9)

The reason why this approach works is that even small changes in $C_I^{\mathrm{eq}}$ considerably affect ${\xi }_{\mathrm{crit}}$ (7) since the difference of two close values $C_V^{\mathrm{eq}}$ and $C_I^{\mathrm{eq}}$ is in the denominator. Table 4 provides a summary of the adjustments. Except for sets B and G, which originally had substantially higher ${\xi }_{\mathrm{crit}}$ than the rest, the largest change of $C_I^{\mathrm{eq}}$ was only about 8%, therefore adjustment of the diffusivity to keep the same transport capacity was not performed.

Also shown in Table 4 is ${\xi }_{\mathrm{crit}}$ for both the original and adjusted parameter sets. While ${\xi }_{\mathrm{crit}}$ varies, a value of $1.32\times {10}^{-3}\mathrm{c}{\mathrm{m}}^2{\min }^{-1}{\mathrm{K}}^{-1}$ is typical for adjusted cases both without and with a stress influence, which is more consistent, compared to the original PD parameters. As illustrated in Figure 6, after adjusting parameter sets with a stress influence, they are in a good agreement regarding ${\xi }_{\mathrm{crit}}$ for compressive stresses around ${\sigma }_{\mathrm{ave}}=-20\mathrm{M}\mathrm{P}\mathrm{a}$. The different slopes for sets I and J are due to low $C_V^{\mathrm{eq}}-C_I^{\mathrm{eq}}$ and low $\mathrm{d}{H}_V^{\mathrm{f}}/\mathrm{d}\sigma$, respectively.

Figure 7 shows the obtained PD distributions, which are now much more consistent and generally are in a better agreement with the experiment. There are two exceptions: set G could not reproduce the IVB in the cone and body phases due to too high $C_V$ (even for a higher ${\xi }_{\mathrm{crit}}$ the IVB in the beginning of body phase is expected to be missing), and set I which, despite $C_I^{\mathrm{eq}}$ changes by only 1.3%, has an unrealistically strong stress influence on ${\xi }_{\mathrm{crit}}$ (due to an order of magnitude lower $C_V^{\mathrm{eq}}-C_I^{\mathrm{eq}}$, compared to other sets) (Figure 6) and a rather “flat” interstitial distribution with almost no vacancies. These sets will be omitted from further discussion.

As intended, the adjusted parameter sets predict the correct radial position of the IVB in the body phase at $L=600$mm. However, for all cases without the thermal stress influence, the IVB (i) was shifted away from the axis at other body phase positions and (ii) was not present in the shoulder region.

For set J with a thermal stress influence the agreement in the body phase is improved but the IVB is still absent at the shoulder. Finally, the best results are achieved by the two last sets with a stress influence—K and especially L—although some differences are still present in the shoulder and initial part of the body phase.

The thermal stress influence is crucial for the correct prediction of the IVB, as can be seen by comparing set pairs C & J, and H & K. As was already observed for the original PD parameters, the IVB becomes straighter, which is in a better agreement with experiment. A worse performance of set J, compared to K and L, is most likely due to weaker stress parameters, see Table 3 and Figure 4. We believe that adjustment of $\mathrm{d}{H}^{\mathrm{f}}/\mathrm{d}\sigma$ can further improve agreement with the experiment, but such investigation is beyond the scope of this paper.

Since various authors use different thermophysical properties and modeling software in their heat transfer and PD simulations, the PD parameters were adjusted in the present study. The same global heat transfer was considered for all cases, allowing us to focus solely on the point defects. As a result, the differences between the adjusted PD parameter sets were rather small compared to the original ones, which could indicate their “universality”, as confirmed by a more consistent value of ${\xi }_{\mathrm{crit}}$ in Table 4. The main takeaway is that the PD parameters alone are insufficient—they have to be considered in conjunction with the thermophysical parameters.

Although our simulations were based on a single experiment, a wide variety of thermal field and pull rate combinations occurred during the strongly transient full crystal growth, therefore we believe the adjusted model can be applied to other hot-zones and growth conditions. In future, it would be beneficial to validate the PD parameter sets against several crystals, especially of larger diameter, to fine-tune the parameters for thermal stress influence.

5. Conclusions

A detailed comparison of 12 published point defect parameter sets has been carried out for Cz growth of an industrial-size Si crystal. The PD defect distributions obtained by different parameter sets were very diverse and none was able to reproduce the experiment, therefore the parameters were adjusted. The IVB, especially at the shoulder and early body phase, could not be satisfactorily predicted by adjusted PD sets without the thermal stress influence. A good agreement with experiment was achieved only for parameter sets with the stress influence and sufficiently high formation enthalpy change. The contribution of an inaccurate heat transfer simulation to the observed discrepancy between the calculated and experimental IVB at the shouldering cannot be ruled out.