Stability Analysis of the First-Order Steady- State Solution in the Czochralski Crystal Growth Process Using Perturbation Techniques

The Czochralski crystal growth is the most commonly used technique for growing crystals. Fluctuations are often encountered in this process and these are believed to be the cause of in homogeneities in the crystals that are undesirable in many applications.

A number of theoretical and experimental studies were carried out in an attempt to reduce or control the fluctuations in the Czochralski crystal growth process. This paper extends the preceding efforts to develop methods to explore the source(s) of these fluctuations by demonstrating that the first-order steady-state solutions of the rotational velocity and temperature are stable and do not exhibit any unusual patterns for any values of the configuration parameters in the Czochralski crystal growth process.

The melted silicon is considered to be viscous and incompressible, with constant fluid properties except for the density change with temperature. It resides in a cylindrical crucible of radius $R_c$ and depth H, Fig. 1. The crystal rod used to form the silicon ingot has a radius $r_c$ and an angular velocity $\omega$.

The compressibility of the melted silicon is ignored except for the buoyancy effect due to the density change with temperature [1]. Based on these assumptions the governing equations simplify to:

Equation (1) represents the continuity equation, eqn (2)-(4) stand for the momentum equations in the $r, \theta$, and z directions respectively, and eqn (5) represents the energy equation.

At the bottom and wall of the crucible, the fluid velocity vector is zero and the temperature is the wall temperature ( $\mathrm{T}_{\mathrm{H}}$ ). At the melt surface and under the crystal, the temperature is the crystal temperature ( $\mathrm{T}_{\mathrm{c}}$ ), the fluid transverse velocity is $\omega$ times the radius of the crystal, and the remaining velocities are zero. At the free boundary, the first-order derivative of the temperature and velocity vector with respect to z are zero. At the surface of the melt, the pressure is atmospheric ( $\mathrm{p}_{\text {atm }}$ ). Below the center of the crystal, the fluid rotational velocity, the first-order derivative of the other velocities, and temperature with respect to $r$ are zero. The top annular fluid surface between the crucible and the crystal is a free surface.

Based on the above descriptions, the mathematical expressions for the boundary conditions are then as follows:

$\mathrm{z}=0, \quad 0<\mathrm{r}<\mathrm{R}_{\mathrm{c}}, \quad \mathrm{v}_{\mathrm{r}}=\mathrm{v}_\theta=\mathrm{v}_{\mathrm{z}}=0, \quad \mathrm{~T}=\mathrm{T}_{\mathrm{H}}$

$\mathrm{z}=\mathrm{H}, \quad 0<\mathrm{r}<\mathrm{r}_{\mathrm{c}}, \quad \mathrm{v}_{\mathrm{r}}=\mathrm{v}_{\mathrm{z}}=0, \mathrm{v}_\theta=\omega \mathrm{r}, \mathrm{p}=\mathrm{p}_{\mathrm{atm}} \quad \mathrm{T}=\mathrm{T}_{\mathrm{c}}$

$\mathrm{z}=\mathrm{H}, \quad \mathrm{r}_{\mathrm{c}}<\mathrm{r}<\mathrm{R}_{\mathrm{c}}, \quad \frac{\partial \mathrm{v}_{\mathrm{r}}}{\partial \mathrm{z}}=\frac{\partial \mathrm{v}_\theta}{\partial \mathrm{z}}=\frac{\partial \mathrm{v}_{\mathrm{z}}}{\partial \mathrm{z}}=0, \mathrm{p}=\mathrm{p}_{\mathrm{atm}} \quad \frac{\partial \mathrm{T}}{\partial \mathrm{z}}=0$

$r=0, \quad \frac{\partial v_r}{\partial r}=v_\theta=\frac{\partial v_z}{\partial r}=0, \quad \frac{\partial T}{\partial r}=0$

$\mathrm{r}=\mathrm{R}_{\mathrm{c}}, \quad \mathrm{v}_{\mathrm{r}}=\mathrm{v}_\theta=\mathrm{v}_{\mathrm{z}}=0, \quad \mathrm{~T}=\mathrm{T}_{\mathrm{H}}$

First-order perturbation with respect to $\omega$ and $\Delta$ are considered in the analysis, where $\Delta$ is the difference between $\mathrm{T}_{\mathrm{c}}$ and $\mathrm{T}_{\mathrm{H}}$.

In the steady-state condition, all partial derivatives with respect to time are zero, yielding the following equations:

$\Delta$-perturbation:

$\omega$-Perturbation:

The numerical finite difference solution of the first-order steady-state equations is now presented and discussed. As an example, Fig. 2a-d shows contour plots of the transverse velocity for fixed configurational parameters $\mathrm{H}, r_{\mathrm{c}}$ and $\mathrm{R}_{\mathrm{c}}$, with fixed difference order 8 and a variable spatial discretization step size which decreases by a factor of 2 sequentially from (a) to (d). Similarly, Fig. 3a-d shows an equivalent study for the temperature. It can be noted that for each case, the contours look very similar as the spatial step size decreases which suggest that convergence is achieved. This observation concerning convergence remains true for a range of configurational parameters ( $\mathrm{H}, r_{\mathrm{c}}$ and $\mathrm{R}_{\mathrm{c}}$ ) and may be numerically verified at all locations except for a small neighbourhood near the upper corner (i.e. at $r=\mathrm{R}_{\mathrm{c}}$ and $\mathrm{z}=\mathrm{H}$ ). This exception is discussed next.

To further explore convergence, and especially the high gradients near the upper right corner, the transverse velocity at the surface $\mathrm{z}=\mathrm{H}$ for a range of parameters is plotted in Fig. 4a-d. Similar results are also shown for the temperature in Fig. 5a-d. It is noted that in all cases, all the curves are practically the same almost everywhere as seen in the insets of Figs 4 and 5. To focus on the upper right corner, the main plots in a range that is about $0.2 R_c$ away from that corner is shown. It is noted that as the step size decreases, all the curves seem to tend towards a common curve (shown dashed) that seems to have a singularity at the upper right corner point. The behaviour of these curves is an interesting but typical behaviour of non-uniform convergence. To further characterize the behaviour near the upper right corner an asymptotic analysis of the first-order equations near that corner is considered next.

As a first step to studying the asymptotics near the upper right corner, we change the independent variables from ($r, \mathrm{z}$) to ($\rho, \phi$) where these are defined as:

$\rho=\sqrt{\left(r-R_c\right)^2+(z-H)^2}, \Phi=\tan ^{-1}\left(\frac{z-H}{r-R_c}\right)$

The geometry and boundary conditions associated with these new variables near the upper right corner are shown in Fig. 6a where $u$ can be either the transverse velocity or the temperature. In terms of these variables, the solutions are in the form

$u=c r^p \cos (s \Phi)$, where $\mathrm{c}, \mathrm{p}$, and s are constants to be determined.

Substituting the above into the first-order equations of either the transverse velocity or the temperature, considering the limit as $\rho \rightarrow 0$ and also applying the boundary conditions, it can be determined that c may be any constant, and $\mathrm{p}=\mathrm{s}= \pm 1, \pm 3, \pm 5$, etc. (negative or positive odd numbers).

In particular, the smallest negative exponent is $\mathrm{p}=-1$ and it is expected that only this singular term exists unless some additional source condition is applied at the corner. Therefore, both the transverse velocity and the temperature are expected to behave asymptotically as follows:

$v_\theta o r T \sim \frac{\cos (\Phi)}{\rho}$

Using the above functional form and matching the far-field behaviour in each case, the thick dashed curves shown in Figs 4 and 5 are obtained. It is noted that the match between the asymptotic form and the numerical results is very good in all cases with the numerical results further approaching this asymptotic behaviour as the step size decreases. The finite difference approach does not have a built-in singularity; therefore the corresponding results will never exactly match the singular asymptotic behaviour. The leading singularity which is the term with the smallest negative exponent in the asymptotic analysis for a sector with angle $\alpha$ (Fig. 6b) is given by:

$p=s=-\frac{\left(\frac{\pi}{2}\right)}{a}$

This indicates that the leading singularity decreases as $\alpha$ increases. For example, the leading singularity equals -1 and -0.75 for $\alpha=90^{\circ}$ and $\alpha=120^{\circ}$. This suggests that it is possible to mitigate or reduce the effect of the singularity by redesigning the crucible so that the incident angle of the vessel with the liquid at the free surface has an $\alpha$ that is larger than $90^{\circ}$. Verification of this mitigation potential requires either detailed experimental work or an analysis of the fully nonlinear Navier-Stokes equations.

Based on the above discussions, it can be concluded that the finite difference approach gives convergent results almost everywhere except near the upper right corner. At that point there is a singularity which has been asymptotically quantified and numerically verified. With these results, a discussion of the stability of these steady-state solutions is provided in the next section.

Several researchers have attempted to analyze the stability of incompressible flows [2-6].

From this paper, one can notice that the stability of the full nonlinear equations in the vicinity of the first-order solutions may be unstable even if the first-order solutions themselves turn out to be stable. This observation is important because the results show that the first-order solutions are stable but it is suspected that the singularity at the upper right corner to be a source of turbulence and therefore, at least near that corner, the steady-state solutions may be unstable.

The stability of the first-order solutions depends on the boundary conditions as well as on the governing equations. This is the main reason why a numerical approach to study the stability is adopted. Although the governing partial differential equations are relatively simple, the boundary conditions are of the mixed type at the top boundary and is the main source of difficulty in this problem.

The approach adopted in studying the stability of the steady-state solutions is to spatially discretize but to leave the time dependence as an analytic expression. The spatial discretization is done as before using finite differences and includes the application of the boundary conditions. The result of this approach is a dynamical system of coupled ordinary differential equations where the variables are the transverse velocity or temperature at the grid points. This leads to linear dynamical systems; therefore the steady-state solutions are stable if all real parts of all the eigenvalues are negative if the maximum real part of all eigenvalues is negative. Similar to the steady-state solution, the convergence of the maximum real part of all eigenvalues is studied. Examples of such a convergence study are presented in Table 1 for the transverse velocity and temperature. These examples indicate that, in all cases shown, the steady-state solution is convergent and stable.

Having considered convergence, the relative stability of the solutions as the parameters change is now considered. It is noted that the magnitudes of the eigenvalues, though neither their signs nor their ordering, depend on the scaling used for the time and spatial variables. To study the relative stability, the maximum (or least negative) real part of all eigenvalues over a grid of values of $r_c$ and $H$ were calculated but for fixed $R_c=1$, which sets the relative spatial scale. The results for either transverse velocity or temperature are similar and are shown in Fig. 7. The results indicate that the steady-state becomes relatively less stable, in the sense that it takes relatively longer to be achieved, as the height H of the cylinder decreases and as the rod cylinder $r_{\mathrm{c}}$ increases. Speculatively, the hypothesis is that under such conditions, the solution to the full nonlinear Navier-Stokes equations would show more turbulence near the upper right corner which is the location of the singularity in the first-order equations. In other words, it is suggested that the results may indicate that instability, probably in the form of turbulence, becomes more prevalent as the cylinder height decreases or as the rod radius increases. The effect of increasing turbulence as the rod radius increases is consistent with previous experimental results, although not directly observed. This effect of rod radius on turbulence is also consistent with previous numerical studies of the fully nonlinear Navier-Stokes equations of this problem. It would be interesting to investigate whether the speculation concerning the cylinder height, or rather the liquid height H , can also be corroborated

by either experiments or the numerical solution of the fully nonlinear

Navier-Stokes equations.

Having discussed the convergence of the numerically-based solutions and their stability, sample solutions for each of the transverse velocity and the temperature are presented. Examples of the steady-state solutions are illustrated using contour plots in Fig. 8a-i for the transverse velocity and in Fig. 9a-i for the temperature. The solutions are all topologically similar with the most interesting feature being a very high gradient near the upper right corner (i.e. near $r=\mathrm{R}_{\mathrm{c}}$ and $\mathrm{z}=\mathrm{H}$) as discussed before. One can also notice that the gradients in the z direction are larger as the height of the liquid decreases. This last observation may be related to the decrease in relative stability as the height of the liquid decreases.

In this paper, a numerical study of the first-order equations associated with the Czochralski crystal growth process is presented. Steady-state solutions are presented and it was shown that the numerical results are convergent and that the steady states are stable. Two observations that may be useful in elucidating the behaviour of this process:

There is a singularity in the first-order equations near the upper right corner. This singularity does not exist in the solution of the fully nonlinear equations but rather indicate a source of turbulence in the process at that corner. If this is the case, a way to mitigate the singularity was suggested, and hence the level of turbulence, by redesigning the lip of the cylinder at the upper right corner.

The steady-state solutions become relatively less stable as either the height of the liquid decreases or as the size of the rod increases. Most likely this change in relative stability is related to the process problem encountered with the Czochralski crystal growth process when the rod radius to cylinder radius ratio increases beyond some threshold.

While the results presented in this paper are suggestive, the solution to the first-order equations does not lead to definite conclusions. Therefore, the authors recommend further studies to extend the current one into the nonlinear range either through a second order analysis or through a numerical analysis of the nonlinear Navier-Stokes equations for this process. In parallel, the authors recommend an experimental study to validate and provide insight into the nature of the process with special emphasis to identify what is happening near the upper right corner of the cylinder.