Numerical simulation of the oxygen distribution in silicon melt for different argon gas flow rates during Czochralski silicon crystal growth process

The effects of argon gas flow rate on the oxygen concentration in Czochralski (CZ) grown silicon crystal were examined. To analyze the influence of the argon gas flow rate in CZ growth process, a 200 mm length silicon single crystal was grown. Different argon gas flow rates are considered. The melt flow pattern, temperature and oxygen concentration distributions in the melt and crystal-melt interface are calculated. The results show that the transport of oxygen impurity is quite dependent on the flow motion in the melt. As the argon gas flow rate increases, there is no fundamental change in flow motion of the melt and the oxygen concentration decreases to a minimum value. When the argon gas flow rate increases further, the flow pattern under the melt-crystal interface starting changes and the oxygen concentration has increased after. Therefore, there is an optimum value for the argon gas flow rate for obtaining the lowest oxygen concentration in the melt.


Introduction
Most of the fabrication of microelectronics use mono-crystalline or silicon single crystal as the base materials.The production of commercial silicon wafers has dominated by Czochralski (CZ) growth process, one of the popular methods for growing silicon single crystal.The popularity of CZ method are by its ability to meet the stringent requirements for high purity, easy of doping, suitable electrical and mechanical properties, and satisfactory crystallographic perfection [1].Heat and mass transport by melt convection play an important role in the quality of crystal growth.During the last decades, much effort has been made in describing global heat transfer in CZ system, impurity transport, and defect in the silicon crystal.In the thermal regime of silicon crystal, much attention has been given to the problem of global exchange by radiation and conduction at the same time the inert gas effect on the melt flow was neglected.
The incorporation mechanism of oxygen and how to control its concentration in CZ process are important for both the ultra-large-scale integration (ULSI) device and wafer manufacturing technologies.The various technique has been applied to control the oxygen concentration in CZ silicon wafers.The technique includes source reduction by controlling the lower crucible temperature and lesser melt contact area between melt-crucible interface or gas-melt interfaces; control the evaporation rate of SiO or melt flow pattern by changing the crystal and crucible rotation, or control the argon gas flow velocity and heat shield shape; and the application of various kinds of magnetic fields.
However, controlling the oxygen distribution can be more effective by knowing the effect of the argon gas flow although there are still some challenges to reveal the gas flow effect on the heat transport, melt convection, and melt-crystal interface geometry in CZ system.The aid of the radiation guide a shear stress introduced by the argon gas flow influences the heat and oxygen transport in the melt.The melt flow reconstruction then affected by the interface gas shear stress on the melt free surface which is influenced by the argon gas flow.
Many kinds of research have been carried out on the argon flow effect on the impurities transport in a silicon crystal growth process for silicon single crystal.Machida et al. [2][3] experimentally observed the melt flow changes which is responsible for preventing the oxygen concentration through the free surface by the argon gas flow change.Evstratov et al. [4] showed oxygen concentration in the crystal enhances along the increasing argon flow rate reported that the radial motion of flow pattern near the crystal rotation has more pronounced flow and isolated from the peripheral motion as the argon gas flow has increased.Bingyan et al. [5] reported that the oxygen concentration of the single crystal is reduced along with the increment of the argon flow velocity.On the contrary, when the argon flow velocity enlarges to some extent, the oxygen concentration increases.Kalaev et al. [6] examined the effect of gas flow on melt convection and shows that there is no fundamental changes are seen in the melt flow pattern and turbulence production even for the highest argon gas flow rate.
The aim of this study is to better understand the influence of the argon flow rate in the interface geometry which is related to the flow pattern construction during the grown crystal.The current research focuses on providing the change of the melt flow pattern under crystal-circulation by increasing argon flow rate will possible has the optimum value for the oxygen distribution along the crystal growth process.

Mathematical Model
The schematic diagram in Figure . 1 Shows an industrial Czochralski furnace used in this study.The crucible has an inner diameter 480 mm with thickness of 9 mm and a height of 322 mm.The total amount of silicon in the crucible is 85 kg, which is massive enough to grow a 900 mm crystal in length.The crystal diameter selected to be 8 in.During the growth process, the furnace pressure is preserved at 2666 Pa and the crystal pulling rate is 51.5 mm/h.The furnace is axi-symmetric and the quasi-steady growth process is assumed.The silicon melt is assumed to be a Newtonian fluid, while argon gas is considered to be an ideal gas.The deformation of the free surface is neglected and the surface tension is assumed to be a linear function of temperature.The differential equations governing heat and mass transports are as follows: where  ⃗⃗    ⃗  are the density, reference density, heat capacity, velocity vector, temperature, thermal conductivity, pressure, stress tensor, gravitational acceleration, impurity concentration and  is the diffusion coefficient of the impurity, respectively.The subscript  maybe or  , where and  indicate the argon gas and silicon melt respectively.The subscript  is O or SiO, where O and SiO indicate the oxygen in the melt and silicon oxide in the argon, respectively.
In the heater: where  is the heat generation from the heater.
The normal velocity components of a wall are zero, while the tangential velocity component of a wall is equal to the wall velocity.Along the free surface, the normal velocity component is equal to zero and tangential velocity component must satisfy the shear-stress balance condition where is the dynamic viscosity,  and are the normal and tangential directions of the free surface, respectively, and is the surface tension of liquid silicon.
The thermal conditions on the interface between two opaque surfaces are (8) The heat fluxes along the interface between the opaque surface and gas can be expressed as follows where  is the Stefan-Boltzmann constant,  is the emissivity, and  in is the incoming radiation heat absorbed by the opaque surface.The view factor methoddeveloped by Dupret et al. [8] is used to compute the radiative heat exchange in the furnace.The temperature at the melt/crystal interface is equal to the melting temperature of silicon, and the energy should satisfy the Stefan condition   (10) where  is the local pulling rate normal to the melt-crystal interface and ∆ is the latent heat.
The segregation effect is taken into account for oxygen impurity.The segregation phenomenon at a melt/crystal interface is described by (12) MATEC Web of Conferences 204, 05013 (2018) https://doi.org/10.1051/matecconf/201820405013IMIEC 2018 Where  is segregation coefficient, and  and  are the oxygen concentration in the melt and crystal, respectively.
Following Matsuo et al. [9], the boundary condition for the oxygen concentration at the crucible wall immersed by silicon melt can be expressed as Following Smirnov and Kalaev [10], the boundary condition for the oxygen at free surface can be expressed as where  O , , , and  are the molar masses of oxygen, silicon, argon and silicon dioxide, respectively, and  is the furnace pressure.The conservation of the concentration flux of oxygen in the melt and SiO in argon at the free surface is The SiO concentration at inlet is zero the flux in the SiO concentration at the gas/solid interface is zero because the effect of deposition is ignored.The Reynolds average Navier-Stoke equation (RANS) is employed with the one equation model to simulate the turbulent motion inside the melt.The flow module of the CGSim package based on Finite Volume Method (FVM) is employed to solve the continuity, momentum, energy and species equations with proper boundary conditions.the equations for oxygen transport are approximated by the fourth order accuracy for the convective and diffusive terms, while the equations for momentum and energy conservations are discretized by second order accuracy for the convective and diffusive terms [9].SIMPLEC is introduced for pressure correction.To capture the sharp oxygen concentration boundary layer, the grids near the wall are refined.The residuals for temperature, velocity and concentration in each control volume are represented by  i , where  may be the temperature, velocity or concentration.The average residual is defined by  ∑  , where is the number of the total control volume.In the present computation, the average residual is restricted to be smaller than 10 -4 .

Results and discussion
A series of numerical simulations is carried out for different argon gas flow rates.The crystal and crucible rotation rates are fixed at 13 and -10 rpm, respectively.The distribution of temperature and oxygen concentration together flow pattern in the melt for different argon gas flow rates is described in Figure .2 for the 200 mm length crystal.As can be seen from the flow pattern, under the crystal-melt interface generated by crystal rotation produces a pair of vortices.Due to the velocity change of the melt around the center of the crucible emerges the movement of the melt under the crystal upward to the crystal-melt interface and then move down near the crucible wall.A buoyancythermocapillary vortex emerges near the crucible wall consequences of the temperature gradient in the melt.The temperature gradient induces the density gradient in the melt lead the flow to move upward along the crucible wall to reach the melt-free surface and then move down to the crucible wall.Therefore, a secondary vortex appears between both of the two vortices.When the argon gas flow rate increases, the temperature difference between all cases is quite significant.The melt flow structure changes slightly when the argon gas flow rate reached 60 SLM (Figure 2d).The secondary vortex moves upwards near the free surface.The formation of these vortices might prevent the movement of oxygen impurities towards the crystal-melt interface from the crucible wall.This helps explain why the lowest oxygen concentration is in these argon gas flow rate case.In the latter cases, the strength and size of these vortices changes are not significant.
We can see that the iso-lines (right Figure ) showing near the crucible wall and the gasmelt interface, the oxygen concentration are much denser.The oxygen impurities originating from the silica crucible wall are moved to the gas-melt interface by the flow motion of buoyancy-thermocapillary vortex and vaporize to form SiO. However, oxygen concentration under the crystal is not the same as oxygen concentration on the melt surface near the crystal, and is strongly dependent on the melt pattern.The oxygen concentration on the melt surface increases as the argon gas flow rate increases if a strong Marangoni effect acts on the surface Li et al. [14].The radial oxygen concentration at the different argon gas flow rate at the melt-crystal interface is also plotted in Figure 3. On increasing the argon gas flow rate, oxygen concentration along the crystal-melt is reduced significantly, but when the argon gas flow increases further the oxygen concentration enhances for the maximum value of oxygen concentration.For further increasing argon gas flow rate from 70 SLM to 100 SLM, oxygen concentration in the crystal-melt interface decreases.This mechanism is difficult to understand and might be explained by consider the flow velocity above the free surface, and also the furnace pressure.
The evaporation at the gas-melt interface takes the SiO away from the melt-free surface by argon gas flow motion.The distribution of the oxygen concentration in the melt during the solidification process is strongly affected by the flow pattern in the melt and gas side.Figure 4 show the distribution of argon flow velocity and silicon oxide above the free melt surface at the different argon gas flow rate.The flow velocity of argon gas is one of the key factor of the change in the melt pattern.In radial direction, the dominant melt flow near the free surface is from the crucible wall to the crystal because of the strong thermal convection.By increasing the argon gas flow rate, the movement of the melt near the free surface could be suppressed by the shear stress of argon gas.For higher argon gas flowrate, it drives a flow motion, large enough to reduce the evaporation of SiO (Figure 4d).In this case, after reach the minimum oxygen concentration, the melt from the bottom of the crucible flow underneath the growth interface without releasing a large amount of SiO at the free surface of the silicon melt, and the oxygen concentration of the grown crystal increases as a result (Fig. 4e).In Figure 4d, can be seen that the SiO evaporation above the free surface reaches the lowest value compare to the other cases.It is agree with the results in Figure .2d, where the lowest value for oxygen distribution along the crystal-melt interface is obtained.

Conclusions
Numerical simulations have been performed to investigate the effects of the argon gas flow rate on melt flow pattern and oxygen concentration at the crystal-melt interface.The results show that the transport of oxygen impurity is quite dependent on the flow motion in the melt.By increasing the argon gas flow rate, the melt flow pattern is slightly altered.The oxygen concentration along the crystal-melt interface is reduce significantly and then increase in maximum value.However, these analysis have been conducted for 200 mm length crystal.The simulations will be further extended in order for higher crystal length by considering the vertical gradient of melt velocity and diffusion layer for SiO transportation..

Fig. 1 .
Fig. 1.Schematic diagram of the major components of the CZ furnace.