The preparation of Czochralski (CZ) monocrystalline silicon is always accompanied by chemical reactions under high-temperature conditions, which produces oxygen and carbon impurities, as well as the deposition and etching of their compounds^[1-4]. Impurities such as oxygen atoms, carbon atoms, and their compounds are unfavorable for growing high-quality monocrystalline silicon. Thus, while growing the crystals, it is necessary to control and reduce these impurities and contamination on the thermo-field components in the crystal-melt and the furnace. It is vital to explore a suppressing method for oxygen and carbon impurities deposited in the thermal system when preparing high-quality large monocrystalline silicon. In the past few decades, oxygen and carbon transporting in the growth of CZ monocrystalline silicon has been extensively explored^[5-7]. Gao et al.^[8] reported the development of a coupled transport model for SiO and CO in argon gas and C and O in silicon melt. Lan et al.^[9] pointed out that silicide depositions on heat shield and SiO particles deposited on the outer surface of the shield fall into the melt, which results in the decreased quality of the monocrystal. Subsequently, Vorob’ev et al.^[10] used advanced chemical models to capture the main trends in parasitic sediment formation observed in the experiments. However, research of the evolution and distribution of deposits (i.e., SiO₂) during crystal growing are seldom reported because monocrystalline silicon grows in a high-temperature, vacuum, and closed environment, and the temperature changes in the thermo-field components cannot be detected when the crystal is growing. In addition, the chemical reactions during the growth will consume graphite elements, which reduces the life of the thermo-field components. To seek an effective suppression method, our study explores the formation of carbon etching and the evolution of silicide deposition on the shield surface and heater when the crystal is growing. The finite element method is applied to numerically predict chemical deposition reactions during crystal growth. Meanwhile, the predicated results have been verified by experiments followed. Moreover, theoretical research and numerical simulation are combined for structure optimization of the thermal system in CZ monocrystalline silicon. Furthermore, we also analyzed chemical reactions on the shield and heater surface, which changes because of argon flow rate and pressure varied in the furnace.

Oxygen impurities in CZ monocrystal furnace come from the dissolving quartz crucible. The dissolved oxygen and silicon atoms are combined at the gas/melt interface, leading to gas-phase volatilization of silicon monoxide (SiO), which is discharged via the argon gas flow. Several experimental observations have confirmed that during this process, the heat shield and heater of graphite component in thermal system react with SiO to form deposits such as CO, SiC, and SiO₂. However, part of the CO is transported back to the melt surface, which is then decomposed into carbon and oxygen atoms, entering the silicon melt and contaminating silicon material while affecting the quality of monocrystalline silicon. According to crystal growth, a numerical model of oxygen and carbon transport has been set up by the finite element method. In the gas region, SiO gas from the surface of free melt flows with argon and reacts with graphite elements and the heat shields nearby, in the temperature range of 1300–1800 K, and CO and SiC are generated^[11-13]:

$ {\rm{SiO}}\left( {\rm{g}} \right) + {\rm{C}}\left( {\rm{s}} \right) \leftrightarrow {\rm{CO}}\left( {\rm{g}} \right) + {\rm{SiC}}\left( {\rm{g}} \right). $  (1)

In the same temperature range, the highly-concentrated CO generated in argon will further react with silicon vapor to produce SiC and SiO:

$ \mathrm{C}\mathrm{O}\left(\mathrm{g}\right)+2\mathrm{S}\mathrm{i}\left(\mathrm{g}\right)\leftrightarrow \mathrm{S}\mathrm{i}\mathrm{C}\left(\mathrm{s}\right)+\mathrm{S}\mathrm{i}\mathrm{O}\left(\mathrm{g}\right). $ (2)

Meanwhile, under the low-temperature condition (below 1600 K), SiO is unstable and decomposes in SiO₂ and Si:

$ 2\mathrm{S}\mathrm{i}\mathrm{O}\left(\mathrm{g}\right)\leftrightarrow {\mathrm{S}\mathrm{i}\mathrm{O}}_{2}\left(\mathrm{s}\right)+\mathrm{S}\mathrm{i}\left(\mathrm{g}\right)\text{,} $ (3)

where g indicates gas, and s refers to solid particle. As shown in reactions (1) to (3), while SiO vapor is carried away by argon gas stream, different reactions are at different positions of graphite components, such as heat shields and heaters, according to the changes of thermo-field temperature and SiO concentrations. Therefore, the highly-concentrated CO gas, SiC particle precipitation, and SiO₂ particle precipitation are formed. Fig. 1 illustrates the furnace inside: chemical reactions are in the argon circuit area; the left-hand side reveals temperature distributions, and the material settings are on the right-hand side.

The in-furnace global simulation calculation has two parts: (1) concentration equilibrium equation for both carbon and oxygen atoms in the melt; and (2) concentration equilibrium equation for CO and SiO in argon. The mathematical model of oxygen–carbon transport includes oxygen–carbon transport in silicon melt, SiO and CO transport in argon, and boundary conditions of the equation.

For the melt, the governing equations for oxygen and carbon transport are as follows^[14]:

$ \frac{\partial \left({\rho }_{\mathrm{S}\mathrm{i}}{c}_{\mathrm{O}}\right)}{\partial t}+\frac{\partial }{\partial {x}_{i}}\left({\rho }_{\mathrm{S}\mathrm{i}}{{u}_{\mathrm{i}}c}_{\mathrm{O}}\right)={D}_{\mathrm{e}\mathrm{f}\mathrm{f}}\frac{{\partial }^{2}{c}_{\mathrm{O}}}{\partial {x}_{i}^{2}}\text{,} $ (4)

$ \frac{\partial \left({\rho }_{\mathrm{S}\mathrm{i}}{c}_{\mathrm{C}}\right)}{\partial t}+\frac{\partial }{\partial {x}_{i}}\left({\rho }_{\mathrm{S}\mathrm{i}}{{u}_{\mathrm{i}}c}_{\mathrm{C}}\right)={D}_{\mathrm{e}\mathrm{f}\mathrm{f}}\frac{{\partial }^{2}{c}_{\mathrm{C}}}{\partial {x}_{i}^{2}}\text{,} $ (5)

where represents the molar concentration of carbon atoms in the melt; represents the molar concentration of oxygen atoms in the melt; is the flow vector of the silicon melt, and is the effective dynamic-diffusion coefficient of oxygen and carbon atoms in silicon raw material.

For the argon gas flow, the governing equations for SiO and CO transport are as follows^{[15, 16]}:

$ \frac{\partial \left({\rho }_{\rm{Ar}}{\omega}_{\rm{SiO}}\right)}{\partial t}+\frac{\partial }{\partial {x}_{i}}\left({\rho }_{\rm{Ar}}{{u}_{\rm i}{\omega }}_{\rm{SiO}}\right)={D}_{\rm{SiO}}\frac{{\partial }^{2}{\rho }_{\rm{Ar}}{\omega}_{\rm{SiO}}}{\partial {x}_{i}^{2}}\text{,} $ (6)

$ \frac{\partial \left({\rho }_{\rm{Ar}}{{\omega }}_{\rm{CO}}\right)}{\partial t}+\frac{\partial }{\partial {x}_{i}}\left({\rho }_{\rm{Ar}}{{u}_{\rm i}{\omega }}_{\rm{CO}}\right)={D}_{\rm{CO}}\frac{{\partial }^{2}{\rho }_{\mathrm{Ar}}{{\omega }}_{\rm{CO}}}{\partial {x}_{i}^{2}}\text{,} $ (7)

where and are mass fractions of SiO and CO in the mixed gas, and and are diffusivities of SiO and CO, respectively.

Temperature linear functions of SiO, CO, and SiC is supportive for getting Gibbs free energy of the chemical reaction Eq. (3):

$ {\Delta G}^\circ =-81300+3.02T\ \ \mathrm{J}/\mathrm{m}\mathrm{o}\mathrm{l},\quad T \lt 1640\mathrm{ }\mathrm{K}, $ (8)

$ \Delta {G^\circ } = - 22100 - 33.1T\ \ {\rm{J}}/{\rm{mol,}}\quad 1640\;{\rm{K}} \lt T \lt 870\;{\rm{K}}, $  (9)

$ {\Delta G}^\circ =-72100-3.44T \ \ \mathrm{J}/\mathrm{m}\mathrm{o}\mathrm{l},\quad T \gt 1687\mathrm{ }\mathrm{K}. $ (10)

Hence, equilibrium constant of the reaction Eq. (3) is:

$ K={\mathrm{e}}^{-\Delta G/RT}. $  (11)

Coupling boundary conditions on the surface of graphite components can be expressed as:

$ {c_{{\rm{Ar}}}}{D_{{\rm{SiO}}}}\nabla \left( {\frac{{{c_{{\rm{SiO}}}}}}{{{c_{{\rm{Ar}}}}}}} \right) = - {c_{{\rm{Ar}}}}{D_{{\rm{CO}}}}\nabla \left( {\frac{{{c_{{\rm{CO}}}}}}{{{c_{{\rm{Ar}}}}}}} \right), $  (12)

$ {c_{{\rm{Ar}}}}{D_{{\rm{CO}}}}\nabla \left( {\frac{{{c_{{\rm{CO}}}}}}{{{c_{{\rm{Ar}}}}}}} \right) = {c_{{\rm{Si}}}}{D_{\rm{C}}}\nabla \left( {\frac{{{c_{\rm{C}}}}}{{{c_{{\rm{Si}}}}}}} \right). $  (13)

These simultaneous equations help us to obtain impurity distributions of oxygen and carbon, as well as their compounds in the monocrystal furnace. It requires extremely fine meshes on both sides of melt-gas to ensure a stable numerical solution.

To clarify the mechanism of oxygen-carbon impurity transport, we conduct a global calculation to simulate CZ-Si crystal growth and impurities distributed in crystals, melts, and gas streams. The simulation selects the equal-diameter growth length of 600 mm as quasi-steady-state modeling to study the chemical deposition, and graphite etching evolution on the surface of heat shield and heater. The parameter settings in calculation are shown in Table 1.

Fig. 2 displays the temperature on heat shield surface, deposits of silicon and its compounds, together with carbon etching distribution. In Fig. 2, when the heat shield temperature rises, the temperature of thermo-field continues to decrease. In the top area where lays a cover, there are no chemical reactions because of the lower temperature. The black curve in the figure reflects simulation calculation of the silicon deposition distributed on shield surface. Thus, it can be seen that, on the surface, chemical reactions took place in a large area, resulting in obvious carbon etching and silicon depositions. At the shield’s upper part, silicon deposition rapidly drops to zero, which exhibits a narrow peak distribution. Meanwhile, carbon etching doubles, which indicates the (2)–(3) chemical reactions and SiO₂ deposition. The heat shield repeatedly used in the experiment is shown in Fig. 2 (left-hand side). The yellow-in-white bright ring represents significant erosion caused by silicide particles deposited on the shield surface.

Under high-temperature conditions, as argon flows, a large amount of unreacted SiO and CO passes through both sides of graphite heater, and react with each other, generating Si and SiC deposits. Fig. 3 shows the deposition rates of silicon and silicide inside the heater and the etching rate of carbon distribution. It has been found that various silicon deposits (white) are formed near the top of the heater. At the same time, carbon etching has also increased linearly. Silicide deposition of higher distribution rate appears at the upper center, while the bottom of the heater is relatively less damaged. By observing the physical map of the heater used for growing monocrystalline silicon, the yellow or dim white is also found on the silicon-rich areas.

To sum up, by comparison, it has been found that the model predicting chemical deposition reactions quite corresponds to the experimental deposition results, which fully demonstrates that the constructed model could be applied to the resolution of practical problems. In the actual CZ crystal preparation, furnace structure design, argon flow rate, and in-furnace pressure are usually changed to suppress chemical reactions in the furnace.

Analysis suggests that an excessive temperature gradient is the main reason for deposition and etching on the heat shield. Therefore, optimizing the furnace structure and changing the thermo-field temperature distribution could be used to explore reaction changes in the furnace. As shown in Fig. 4 (right-hand side; i.e., the optimized furnace structure), a top insulation structure has been added onto the cover.

After optimizing the furnace structure, simulation calculation has been carried out, while setting the reference values both for the in-furnace pressure and flow rate of argon gas at 20 Torr and 100 slm, respectively. Hence, we obtained temperature distribution of heat shield surface and material deposition rate before and after the optimization, as shown in Fig. 5. The surface temperature of the optimized heat shield increases. In particular, the temperature of the top area increases by more than 500 K, as shown in Fig. 5(a). The etchings on the shield surface and cover plate are presented in Figs. 5(b)–5(d). The silicide deposition rate decreases rapidly, the position of chemical reaction changes, and carbon etching and silicon deposition rates decrease. New chemical reactions are taking place in the top cover area.

Two differences before and after the structural optimization:

(1) After optimizing the structure, the deposition rates of silicon and silicide (mainly SiO₂) have been reduced, thereby reducing the etching rate of carbon accordingly. Although the chemical reaction area expands, the sediment thickness has been cut and impurities are difficult to shed off.

(2) After optimizing the structure, the SiO₂ deposition layer moves from the upper part of position 3 to position 2 (Fig. 1 (left)). Meanwhile, new carbon etching and silicon deposition are formed at position 1 on the top cover. In brief, reaction position shifts from the shield to the top cover area, which avoids wear of the heat shield.

The increasing temperature in the furnace causes SiO gas to rise with argon gas flow, and to react at the cover. However, the temperature at the top cover is relatively low and the chemical reaction rate decreases, which is beneficial for suppressing the reaction. As for material loss, etch deposition transfers to the top cover area of the shield, which prolongs service life of the shield. Also, the replacement of the cover is easier and more economical than replacing the heat shield.

After optimizing the furnace structure, the deposition distribution on the heater surface changes accordingly. As shown in Fig. 6, the deposition rate on the surface decreases. However, compared to the original structure, new silicide deposition and carbon etching are found on the heater surface after the optimization. Moreover, the deposition area is relatively large, which greatly shortens the service life of the heater.

In-furnace pressure is the key parameter to study the SiO volatilization above the free melt interface. To explore the effects of the furnace pressure on chemical deposition etching, a series of simulation calculations with pressures ranging from 18 to 25 Torr have been performed at a fixed flow rate of 100 slm. The position is selected where the maximum deposition rate occurs. The comparison results before and after the structure optimization have been given; i.e., the original structure takes 355 mm as the observation distance. Meanwhile, the optimized structure adopts 424 mm for the same purpose, as shown in the positions marked by the red dots in Fig. 4.

Fig. 7 shows silicide distribution and carbon etching on the shield under different furnace pressures of the original structure. As shown in Fig. 7, when furnace pressure increases, silicide deposition rate decreases rapidly. Meanwhile, a significant decrease has been observed for carbon etching rate. This indicates that reaction (1) balance shifts to the left, CO diffusion is suppressed, carbon consumption is reduced, and carbon pollution in the melt is fundamentally reduced. In addition, calculations show that when the furnace pressure is greater than 22 Torr, in-furnace reaction maintains equilibrium, and the rates of carbon etching and silicide deposition remain unchanged.

Fig. 8 shows silicide distribution and carbon etching on the shield under different furnace pressures after the structure optimizing. It can be found that when compared with the original structure, as the furnace pressure increases, the silicide deposition rate and the carbon etching rate change little and there is almost no change. The main reason is that temperature rise at upper cover of heat shield and the partial pressure of the atmosphere decreases after the optimization of the structure, which slows down the reaction. Therefore, the change in pressure has little influence on the silicide deposition and carbon etching rate of the upper cover of the heat shield under the optimization of the structure.

To clarify the effect of argon flow rate influenced on the deposition of oxygen and carbon impurities in monocrystal furnace, monocrystalline silicon is growing at a fixed furnace pressure of 20 Torr and with different inlet argon rates. We analyzed the effects of furnace gas flow on the etching of heat shield and deposition distribution of silicide. The flow rate has been set to be 80, 100, and 130 slm, respectively.

Fig. 9 illustrates a series of simulation results, which explains the reverse diffusion and gas flow control of CO production at different argon flow rates. In addition, Figs. 9(a)–9(c) show the increases in CO concentration in the gas at flow rates of 80, 100, and 130 slm after optimizing the insulation structure. As shown in Fig. 9, compared with the original structure, in the optimized structure, the CO gas distribution shifts from the middle area to the top area of the shield. In addition, the position of chemical reaction shifts, which slows down the wearing of the shield. This result is also consistent with the conclusion of Section 4. Thus, as the flow rate increases, the CO concentration decreases, and a large amount of CO gas is taken out of the furnace with argon gas flow.

Fig. 10 shows a series of simulation results, which explains the increasing concentration of SiO in the gas after optimizing the insulation structure. Fig. 10 shows the increases in SiO concentration in the gas when flow rates are at 80, 100, and 130 slm after optimizing the structure. As shown in Fig. 10, after the optimization, SiO gas in the upper part of the shield decreases, and SiO concentration in the gas decreases with the continuous increasing of the flow rate.

Theoretically, partial pressure of SiO in the furnace gas and argon flow , SiO evaporation flow Q₂ at the free interface of the melt, and the pressure in the furnace satisfy the following relations:

$ {P}'=\frac{P{Q}_{2}}{{Q}_{1}+{Q}_{2}}\approx \frac{P{Q}_{2}}{{Q}_{1}}. $ (14)

When the flow rate of argon gas increases, SiO gas in the atmosphere is rapidly evacuated, resulting in decreasing partial pressure of SiO. Meanwhile, the evaporation of SiO from the silicon melt has been enhanced; therefore, the oxygen concentration in the crystal decreases. However, if the flow of argon increases to a certain extent, the continuous blowing of argon will exert two effects:

This generates shear stress, which affects natural convection in the melt, thereby indirectly affecting the doping of oxygen and carbon atoms.

As a result, it lowers the free surface temperature of the silicon melt and slows down the evaporation of SiO gas, which is not conducive to impurity discharging.

To explain the evolution of carbon etching and silicide deposition on the heat shield surface and heater during CZ-Si crystal growth, an oxygen-carbon chemical model has been constructed to predict deposition. Simulation results and experiments of the shield and heater have verified the measurement results. Then, the structure of hot area of the CZ furnace has been optimized for chemical reaction temperatures. At the same time, when we change the flow rate of argon gas and the pressure in the furnace, the studies had also been made to analyze the changes in chemical reactions on the heat shield surface and heater. In summary, we finally come to the following conclusions:

(1) After optimizing the furnace structure, the deposition rates of silicon and silicide (mainly SiO₂) on the shield decrease significantly. The position where chemical reactions take place shifts from the shield to the top cover area, thereby avoiding the wearing of the shield and reducing impurity shedding. However, as for the heater, new silicide deposits and carbon etching are found on the surface, which greatly shortens the service life of the heater.

(2) By increasing in-furnace pressure, the rate of silicide deposition decreases rapidly. At the same time, the rate of carbon etching also decreases. Fundamentally, the consumption of carbon is suppressed, and the carbon pollution in the crystal is reduced.

(3) After optimizing the insulation structure, the concentrations of SiO and SiO in the gas decrease with the increasing argon flow rate. Meanwhile, the increasing gas flow rate effectively suppresses the generation and reverse diffusion of SiO.

This work was supported by the National Natural Science Foundation of China (No. 61533014), and the Natural Science Foundation of Shaanxi Province (No. 2019JQ-734).