Optimal Cooling System Design for Increasing the Crystal Growth Rate of Single-Crystal Silicon Ingots in the Czochralski Process Using the Crystal Growth Simulation

Abstract

Here, we report a successfully modified Czochralski process system by introducing the cooling system and subsequent examination of the results using crystal growth simulation analysis. Two types of cooling system models have been designed, i.e., long type and double type cooling design (LTCD and DTCD) and their production quality of monocrystalline silicon ingot was compared with that of the basic type cooling design (BTCD) system. The designed cooling system improved the uniformity of the temperature gradient in the crystal and resulted in the significant decrease of the thermal stress. Moreover, the silicon monocrystalline ingot growth rate has been enhanced to 18% by using BTCD system. The detailed simulation results have been discussed in the manuscript. The present research demonstrates that the proposed cooling system would stand as a promising technique to be applied in CZ-Si crystal growth with a large size/high pulling rate.

1. Introduction

Most single crystal silicon ingots are grown using the Czochralski process. In this process, monocrystalline silicon is manufactured by contacting a silicon seed crystal from molten silicon and pulling it out slowly using a pulling system with crucible rotation [1]. The uptime of the Czochralski process is quite long and expensive. Therefore, the current industry aims to increase productivity and economic efficiency by increasing the speed of crystal growth. To improve productivity, it is necessary to reach a high crystal growth rate, but it is difficult to maintain the quality of the silicon crystal at a high crystal growth rate from a predetermined process. Furthermore, the twisting phenomenon of the crystal has limited attempts to achieve a high crystal growth rate in a given hot zone. A higher crystal growth rate than the given value can result in a distortion of the lattice inside the crystal, which causes dislocations and deteriorates the crystal quality [2,3,4,5]. Therefore, our research proposes an implementation plan to improve the growth rate of the crystal as well as internally improve the chamber. The entire process of the Czochralski method proceeds in the following order: dipping, necking, shouldering, body growth, and tailing. To grow high-quality silicon ingots, certain variables and von Mises stress must be considered. Moreover, the crystal defect area V/G (crystal growth rate (V)/temperature gradient (G) of the silicon melt–crystal interface region) must be considered [6,7]. The crystal growth rate must be able to pull faster while maintaining quality from the crystal’s internal temperature gradient and the free defect area of the melt–crystal interface. Improving the pulling speed of a single crystal silicon ingot will need to consider structural improvements. CG-Sim software is a simulation software for growing single and polycrystals from a melt of crystals. The software is used to achieve the optimal design, process conditions, and crystal quality improvement through heat transfer, temperature analysis, as well as gas and melt convection analysis of the crystal defects and stress calculations. The mathematical modeling represents the equations for the convection of gas flow and fluid in Czochralski furnace and the heat transfer phenomena through hot zone [8].

2. Mathematical Modeling and Experiment

The Czochralski process structure is shown in Figure 1. We obtained cooling system using the Crystal Growth Simulation. The geometry of the Czochralski growth system is symmetrical in the axial direction, but the nature of the growing crystal is isotropic. Therefore, it is not easy to achieve full axial symmetry since it is too complicated to match the central axis of the triple point in the growth process. Therefore, this study assumes axial symmetry by reducing the asymmetrical pattern using rotation. The Czochralski growth simulation utilizes two-dimensional cylindrical coordinates [9]. In the Czochralski process, studies on computational models, such as the Reynolds number (Re), Grashof number (Gr), and Prandtl number (Pr), have been conducted [10]. To analyze the temperature distribution and fluid flow of the Czochralski process, there are (1) continuous equations of differential equations, (2) momentum equations for each velocity component, and (3) energy equations. By assuming incompressible fluids and expressing them in tensor form in Cartesian coordinates, the explanation can be expressed as follows:

$$\frac{\partial \mathsf{\rho }}{\partial \mathrm{t}}+\nabla ·\left(\mathsf{\rho }\mathrm{u}\right)=0$$

(1)

$$\mathsf{\rho }\frac{Du}{Dt}=-\nabla p+\nabla ·\tau +\rho g$$

(2)

$$\mathsf{\rho }\frac{\partial e}{\partial t}+p\left(\nabla ·u\right)=\nabla \left(\mathrm{k}\nabla \mathrm{T}\right)+\mathsf{\phi }$$

(3)

where $\frac{D}{Dt}$, ρ, u, $\nabla$, p, and t are the convective derivative, the density, velocity, divergence, pressure, and time, respectively. $\mathsf{\tau }$ represents the deviatoric stress tensor, which has an order of two, while g is the gravitational acceleration acting on the continuum, φ is the viscous-dissipation function, and k is the thermal conductivity. For the dimensionless number, in this article, the rate of crystallization of silicon corresponds to heat flux, transition in crystals and melts from the following equation [11]:

$${\mathsf{\rho }}_{\mathrm{crystal}}∆{\mathsf{H}\mathrm{u}}_{\mathrm{n}}={\left(\mathsf{\lambda }\frac{\mathrm{dT}}{\mathrm{dn}}\right)}_{\mathrm{melt}}-{\left(\mathsf{\lambda }\frac{\mathrm{dt}}{\mathrm{dn}}\right)}_{\mathrm{crystal}}$$

(4)

where u_n is the average rate of partial crystallization at the melt–crystal interface, ∆H is the heat of crystallization of 1.8 × 10⁶ (J/kg), and ρ_crystal is the solid silicon thermal conductivity at the melting point temperature 66.5 (W/km).

The average crystallization rate υ was calculated by integrating the space S with u_v representing the perpendicular crystallization rate to the partial crystallization rate u_n at the melt–crystal interface from following the equation:

$$\mathsf{\upsilon }=\frac{1}{\mathrm{s}}{\displaystyle \int }{u}_{\mathrm{v}}\mathrm{ds}$$

(5)

Flow due to surface tension occurs in the direction of the low-temperature part at the high-temperature part where the surface tension is small, i.e., flow occurs in the direction of the interface from the crucible wall. Therefore, flow caused by a difference in density and the flow direction due to a difference in surface tension are mutually synergistic. Hence, the overall flow occurs from the crucible to the interface at the surface portion and from the interface to the bottom at the center. Rotation of the crucible and crystal primarily produces circumferential flow. The dimensionless number that can indicate the effect of this rotation is the Reynolds number $ℛ$e for the rotation. If the speed of crystal and crucible and the characteristic length are Ω_x, Ω_c, and h respectively, the Reynolds number for rotation is defined as follows:

$$\mathcal{R}{\mathrm{e}}_x=\frac{{\mathsf{\Omega }}_x{h}^2}{v}$$

(6)

$$\mathcal{R}{\mathrm{e}}_c=\frac{{\mathsf{\Omega }}_c{h}^2}{v}$$

(7)

While the shape of the flow in forced convection depends on Re, the shape of the flow in natural convection depends on Gr, which represents the ratio of the buoyant to viscous forces acting on the fluid.

$$\mathrm{Gr}=\frac{g\beta \Delta T{R}^3}{v^2}=\frac{\mathrm{Buoyant}\mathrm{force}}{\mathrm{Viscous}\mathrm{force}}$$

(8)

where R refers to the diameter.

Pr indicates the contribution of the momentum diffusivity relative to thermal diffusivity in layer flow. Thermal diffusivity dominates when Pr < 1, and momentum diffusivity dominates when Pr > 1.

$$\mathrm{Pr}=\frac{v}{a}=\frac{\mathrm{viscous}\mathrm{diffusion}\mathrm{rate}}{\mathrm{thermal}\mathrm{diffusion}\mathrm{rate}}=\frac{C_p\mu }{k}$$

(9)

In the Czochralski process, single-crystal silicon crystal growth directly interacts with the flow of dissolved silicon in quartz crucible. The impurities and argon gases introduced into melt–crystal interfaces, heat transfer, molten silicon convection, and temperature gradient are important parameters, and the crystal growth rates influences the temperature gradient inside the growth determination [12]. At the equal temperature gradient, the high rate of crystalline growth resulted in high in vacancy defect with the melt–crystal interface in concave form, whereby the low rate of crystalline growth appeared in high interstitial defect in convex form. Therefore, it provides information to predict the number of defects in the crystal, as it helps to determine the shape of the melt–crystal interface depending on the growth rate of the crystal. The improving in the growth rate should be carried out by considering the quality of the appropriate crystal [13,14,15,16,17].

The numerical simulation used for this experiment was CG-Sim^® (Crystal Growth Simulation) of the STR Group in St. Petersburg, Russia. CG-Sim is a crystal growth data that is a package that collects data from real processes and combines them into a single software. The efficiency of the process can be improved by comparison with actual data. The CG-Sim (melt) software can simulate the process of growing single crystals and multicrystals from the melt of a crystal and can produce different results depending on the internal structure, size of the crystal, and shape of the quartz crucible. The simulation was carried out considering the quality determination of the crystal grown Si-ingot according to the design of a new cooling system and the method to realize a high pulling rate.

Table 1. Properties of the material.

Materials	Heat Capacity (J/kg)	Heat Conductivity (W/m)	Emissivity
Si(crystal)	1000	75.22	0.7
Si(melt)	915	66.50	0.3
Graphite	720	105.26	0.7
Ar	521	3.06	–
Quartz crucible	1232	4	0.85
Steel	438	15	0.45
Carbon felt	–	0.02	0.9
Water	4200	0.6	0.5

lists the properties of the main substances, including the thermal properties and emissivity. The unique properties of silicon, graphite as a cooling system pipes, inert gases, such as argon, quartz crucibles (silica) and iron (steel), carbon felt as a thermal insulator, and water as a cooling system, which all influence the growth of the crystals, are used [18,19,20].

3. Results and Discussions

Initially, the crystal growth simulations (CG-Sim) were performed for basic model crystal design (BTCD) of the Czochralski process (CZ-process) for clear understanding of the parameters that affect single-crystal production yield. Figure 2 describes the thermal gradient in the CZ-process with different crystal position (CP) heights. The growth rate and defects of the single-crystal silicon ingot were studied using CG-Sim analysis at various crystal position heights with a pulling speed of 1.0 mm/min. Figure 3 illustrates the V/G ratio, von Mises stress, and power consumption with different CP heights at various crystal lengths at 1.0 mm/min of fixed pulling speed in CZ-process. The corresponding results are listed in

Table 2. Analysis result of the crystal positions.

Crystal Position	V/G (cm2/min/K) (Average)	von Mises Stress (MPa) (Average)	Power (kW)	Crystal Mass (kg)
CP400	0.00232	15.09	92.12	29.25
CP800	0.00239	15.33	89.08	62.41
CP1200	0.00236	15.07	87.56	95.77
CP2000	0.00235	15.08	86.59	162.52
CP2400	0.00236	14.99	86.31	196.60
CP3000	0.00238	14.71	87.75	246.07

. von Mises stress has been considered as an effective or equivalent stress. When an object is subjected to a force or moment from outside, it would tolerate that force to some extent, but if it is beyond certain size, it would no longer support the external force, hence gets destroyed. The condition that predicts such destruction is termed as “yield criterion.” The yield criterion is a value that indicates the maximum distortion energy at each point of an object under stress. It is based on various failure theories, but the maximum distortion energy theory, which has been considered the most accurate, uses von Mises equivalent stress. Plasticity occurs when the calculated equivalent stress value is greater than the yield strength of the crystal. The simulation package of the Czochralski process would judge the stress based on the number of Re, Gr, and Pr. The amount of change in stress increased slightly from CP400 to CP800 and decreased rapidly with increasing crystal height. When the stress was high, the uniformity of the temperature gradient difference based on the diameter of the crystal was exceeded. The results suggest an increase in the crystal length and the stress in the crystal has decreased with increased cooling effect [21,22,23]. It indicates that the increase of cooling time would improve the crystal quality. Generally, the stress at 30 MPa or less has not lead to stress defects in the crystal growth [24]. However, the silicon crystal growth conditions are defined based on the critical region in which the interstitial region is transformed in the vacancy region. The V/G ratio identified the defect domination area of vacancy or interstitial as shown in Equation (10). The critical V/G can be determined by concentration of vacancy (Cv) − concentration of interstitial (Ci) = 0. Vacancy is agglomerated to form void, Oxidation Induced Stacking Fault (OISF), Crystal Originated Pit (COP), and interstitial agglomerate to form a large dislocation loop [25]. The obtained defect-free area indicated that the vacancy and interstitial are very low concentrations without COP, OISF, and dislocation loops [26,27,28,29,30]. Therefore, when the V/G value is in range between 0.00213 and 0.00219 (cm²/min/k) along with radial and axial direction, the crystal growth is free of defect area and shown in Figure 4 [31,32,33,34].

The V/G ratio has initially increased from CP400 to CP 800 after which it decreased with the increase in CP height. Moreover, the power consumption was decreased from 92.12 to 87.75 kW. The CG-Sim analysis results for various CP heights of single crystal ingot in CZ-process suggest that the increase in the crystal length would decrease the stress in the crystal due to increased cooling time. Eventually, it indicates that the increase of cooling time would improve the crystal quality.

$$\begin{array}{cc}& \frac{V \to Velocity}{G\to Temperature gradient}\hfill \\ \hfill \Rightarrow & {\left(\frac{V}{G}\right)}_{Critical}=0.00216m{m}^{2}mi{n}^{-1}{K}^{-1}\hfill \\ \hfill \Rightarrow & \left(\frac{V}{G}\right)>{\left(\frac{V}{G}\right)}_{Critical}={\left(\frac{V}{G}\right)}_{Interstitial}\hfill \\ \hfill \Rightarrow & \left(\frac{V}{G}\right)<{\left(\frac{V}{G}\right)}_{Critical}={\left(\frac{V}{G}\right)}_{Vacancy}\hfill \\ \hfill \Rightarrow & {\left(\frac{V}{G}\right)}_{large dislocation loop}<0.00213{\left(\frac{V}{G}\right)}_{Interstitial}<0.00213~0.00219{\left(\frac{V}{G}\right)}_{I/V}<0.00219{\left(\frac{V}{G}\right)}_{Vacancy}<{\left(\frac{V}{G}\right)}_{OISF}\hfill \end{array}$$

(10)

Cooling System Length Adjustment and Dual-System Design

In order to understand the cooling effect on the single-crystal production yield and quality, we modified the cooling system in CZ-process with two types of designs and studied their effects using CG-Sim analysis. Figure 5 describes the cooling system BTCD, the long type cooling design (LTCD), and the other one is double type cooling design (DTCD) installed in the Czochralski process, where LTCD was just a 10% extended length of cooling system compared to that of BTCD. Figure 6 shows the CG-Sim 2D view of thermal gradients nearby crystal and illustrated point of advantage by cooling design in CZ-process. The decreased thermal gradients were observed in vertical and horizontal DTCD than LTCD, due to faster cooling rate in DTCP design, which decreased the internal heat of the crystal, thereby maintaining the uniformity.

To clearly understand the crystal quality, the V/G ratio, von Mises stress, power consumption, and crystal mass were calculated by CG-Sim analysis. The corresponding results are displayed in Figure 7 by comparing with the BTCD model.

Table 3. Analysis result of the simulation.

	Type	CP400	CP800	CP1200	CP2000	CP2400	CP3000
V/G (cm2/min/k)	BTCD	0.00232	0.00239	0.00236	0.00235	0.00236	0.00238
	LTCD	0.00230	0.00236	0.00234	0.00234	0.00234	0.00237
	DTCD	0.00222	0.00219	0.00217	0.00216	0.00217	0.00219
von Mises stress (MPa)	BTCD	15.09	15.33	15.07	15.08	14.99	14.71
	LTCD	14.75	14.80	14.98	14.94	15.00	14.49
	DTCD	14.20	13.84	13.82	14.13	14.01	13.59
Power (kW)	BTCD	92.12	89.08	87.56	86.59	86.31	87.75
	LTCD	92.01	89.64	88.20	86.22	86.83	88.15
	DTCD	94.04	91.51	90.58	89.45	89.09	90.86

illustrates the comparative data of BTCD, LTCD, and BTCD obtained in CG-Sim analysis. Based on the results obtained in the BTCD, we have designed two types of cooling systems to understand the effective cooling system on crystal quality. The CG-Sim 2d view of the Cz-process is displayed in Figure 6. In addition, Figure 7 depicts the V/G, von Mises stress, and power and defects in the Si-single crystal ingot by CG-Sim analysis, and the corresponding results are listed in Table 3.

The V/G ratio, which represents the defect-free area was measured. The DTCD-based crystal showed that the DTCD has a lower V/G compared to BTCD, in which the reduction percentage was approximately 8.29%. The V/G was only decreased by 1% compared to LTCD and BTCD. Moreover, the V/G values have approached the defect-free area using DTCD design. This is due to an increase in crystal CP, which slightly decreases the V/G values, thereby approaching the defect-free area, which has been advantageous in obtaining production yield and quality crystals. Similarly, the stress effect was also compared for the DTCD, LTCD, and BTCD designs wherein the decreasing rate of stress was almost 10% for DTCD than BTCD design and 3% for LTCD than BTCD. This is due to the improvement in the cooling efficiency as the distance between the DTCD and the crystal is less, which led to a fast reduction of heat from inside the crystal. As the stress decreases, the stress required to destroy the crystal also decreases with decrease in temperature gradient. However, the power consumption was increased by approximately 3% for DTCD compared to BTCD (Figure 7) because more power has been consumed in order to maintain the hot zone temperature of the molten state from the effect of cooling system inside the chamber.

The obtained results indicated that each CP has improved with decreasing V/G and von Mises stress. The process was considered as a mass production since the CZ-process has required high cost and longtime operation. Therefore, the higher height of CP has the faster cooling effect, and the lower value of stress can be observed from the results. The CP3000 has been selected as a best test target since it has the best cooling effect.

To further analyze the effect of cooling systems design on crystal height, the temperature distribution was analyzed at CP3000 using CG-Sim for all the cooling designs. Figure 8 represents the thermal distribution and heat flux inside the chamber and crystal using BTCD, LTCD, and DTCD cooling systems at CP3000. The lowest temperature distribution was observed for DTCD design compared to the BTCD and LTCD, which was due to the increase in the cooling capacity and influence of some region nearby. Moreover, when the temperature gradient decreases inside the crystal, the crystal growth rate would increase, simultaneously.

The rapid crystal growth rate is important to produce large crystal ingot with improved productivity. Therefore, the simulations were conducted at pulling speeds of 1.0–2.0 mm/min at CP3000. The corresponding comparative plots are displayed in Figure 9.

Table 4. CG-Sim analysis results at various pulling speeds between 1.0 and 2.0 mm/min at CP3000.

Pulling Speed	V/G (cm2/min/K)				von_Mises_Stress (MPa)
	BTCD	LTCD	DTCD	Decreasing Rate %	BTCD	LTCD	DTCD	Decreasing Rate %
1.0	0.00238	0.00237	0.00219	8.67%	14.71	14.49	13.59	8.24%
1.1	0.00256	0.00255	0.00236	8.47%	16.87	16.64	15.09	11.79%
1.2	0.00274	0.00272	0.00253	8.30%	19.32	18.66	17.13	12.78%
1.3	0.00290	0.00288	0.00269	7.80%	21.9	21.36	18.79	16.55%
1.4	0.00308	0.00305	0.00285	8.07%	25.18	24.5	21.33	18.04%
1.5	0.00325	0.00323	0.00301	7.97%	28.78	27.79	24.39	17.99%
1.6	0.00341	0.00338	0.00316	7.91%	32.13	31.59	27.33	17.56%
1.7	0.00359	0.00356	0.00332	8.13%	36.79	35.75	30.81	19.40%
1.8	0.00377	0.00374	0.00347	8.64%	40.71	39.42	35.27	15.42%
1.9	0.00396	0.00392	0.00364	8.79%	44.28	43.1	38.79	14.15%
2	0.00413	0.00411	0.00381	8.39%	48.03	46.48	42.3	13.54%

shows the difference between V/G and von Mises stress, which are expressed as the quality of the LTCD and DTCD crystals according to the pulling speed.

The above results illustrate that the increased pulling speed would move the V/G values towards defect area from defect-free area. However, the V/G values were lower for DTCD than BTCD and LTCD designs, suggesting a possibility to improve the crystal production by modifying the cooling system of DTCD. Furthermore, von Mises stress defects were also greatly decreased for DTCD design compared to the other designs at each pulling speed. The average stress differences between BTCD and LTCD were 1–3%, whereas the comparison of average stress difference between BTCD and DTCD was 15%. As a result, the experimental data has shown the possibility of satisfying the high crystal growth rate while improving the stress, by reducing the temperature gradient inside the crystal from the cooling system. The efficient difference between LTCD and DTCD were significantly indicated from the comparison results with BTCD. In other words, the effect of crystal growth was lower than that of the double design due to the cooling effects, and the pulling speed could be improved. The V/G value of DTCD was approximately 7.5% better than that of LTCD. Moreover, the internal thermal stress of the crystal was improved by approximately 12% because longer residence time from the cooling zone during crystal growth and deeper crystal would result in lower internal temperature gradient of the crystal. With this experience, methods to enhance the pulling speed and growth rate have been developed. This study was aimed at increasing the pulling rate of a single crystal silicon ingot by the Czochralski method using a Crystal Growth Simulation.

4. Conclusions

The crystal growth rate of BTCD and LTCD does not have a large margin of error, a difference of about 0.2 mm/min between the DTCD and the range of velocities could be confirmed through this experiment. That is, the cooling effect had a great influence on the crystal growth rate, and the growth rate could be improved. From this study, the difference between the V/G and von Mises stress was shown for each cooling system design, and the defect-free area V/G value of DTCD was approximately 8% better than that of BTCD. In addition, a longer residence time and faster cooling effect from the cooling zone during crystal growth and a lower internal temperature gradient resulted in an approximately 15% decrease in the internal thermal stress of the crystal. Through this experience, one of the growth enhancement solutions was applied to achieve an improved pulling speed. The purpose of this study was to achieve a high rate of increase in single crystal silicon ingots by the Czochralski method using Crystal Growth Simulation. Although the power was slightly different from the BTCD, the crystal growth rate was improved by 18.1% or more by the modified DTCD. In other words, increasing the crystal growth rate of the single crystal silicon ingot reduces the production time of the process more than the amount of electric power consumed, thereby improving productivity. In addition, efficiencies were achieved from an economic perspective. As a result of this research, further progress was able to be completed in the development of the Czochralski process and step to secure the source technology for achieving economic efficiency and improved productivity of single crystal silicon ingots was implemented.