High-pressure liquid-encapsulated Czochralski growth (HPLEC) of InP

Unlike the conventional Czochralski process, modeling of high-pressure liquid-encapsulated Czochralski (HPLEC) method is particularly complicated, because account must be taken of turbulent gas and melt convection as well as of the encapsulant preventing phosphorus evaporation from the melt. The heat transfer interactions between the recirculating gas, the B2O3 layer and the crystal is very complex and requires especially accurate simulation.

To study the influence of the gas flow on the thermal field in the crystal, we used a 3D unsteady model that can directly resolve turbulent eddies and their contribution to the heat exchange. We considered gas convection together with heat transfer in the crystal and encapsulant flow in an industrial HPLEC InP set-up, paying special attention to the temperature gradients in the crystal as they determine the thermal stresses.

Computational domain for preliminary 2D steady-state simulation of global heat transfer with the computed temperature distribution and flow pattern

Figure 1. Computational domain for preliminary 2D steady-state simulation of global heat transfer with the computed temperature distribution and flow pattern

Computational domain for detailed 3D unsteady simulations

Figure 2. Computational domain for detailed 3D unsteady simulations

The computational procedure included preliminary simulation of global heat transfer in the whole set-up (Figure 1) with the account of all the heat exchange modes — convective, radiative and conductive — performed using a 2D axisymmetric steady-state model. The resulting thermal field was used to set boundary conditions for the 3D computational domain, see Figure 2.

The results demonstrate an essentially non-axisymmetrical flow structure with multi-scale vortices in different areas of the computational domain. The solution is generally characterized by an upward flow along the pulling rod, a downward flow along the cool external wall and two re-circulating zones in the area between the crystal and the crucible, with the most intensive unsteady convection observed at the center of the growth set-up, around the crystal and the pulling rod, see Figures 3 (a, b).

Figure 3 (a). Instantaneous flow patterns and temperature distributions in the gas domain at two different time moments.

Figure 3 (b). Instantaneous flow patterns and temperature distributions in the gas domain at two different time moments.

Strong temperature fluctuations in the gas result in pronounced variations of the temperature and temperature gradients at the crystal periphery. One can see in Figures 4 and 5 that the largest absolute gradients take place at the melt/crystal/encapsulant tri-junction point, while the area, where a large radial gradient occurs, is located at the crystal/encapsulant/gas tri-junction point. Besides, the fluctuations of the radial temperature gradients in the conic part of the crystal can give rise to multiplication of dislocations.

Figure 4. Instantaneous distributions of the absolute temperature gradients at the two time moments.

Figure 5. Instantaneous distributions of the radial temperature gradients at the two time moments.

To estimate the contribution of gas convection to the heat loss from the crystal surface, we compared time-averaged heat fluxes across the crystal and the encapsulant in different regimes, Figures 6, 7. It can be seen that the heat loss from the crystal and the encapsulant increases dramatically with the gas pressure due to the effect of turbulent gas convection. Similar computations made within a 2D quasi-steady axisymmetrical approach, using two various turbulence models, significantly underestimated the heat loss from the crystal and the encapsulant. Summarizing, the nitrogen flow around the crystal greatly influences the axial temperature gradient in the crystal and, furthermore, produces observable unsteady fluctuations of the radial temperature gradients in the conic part of the crystal, which can result in significantly higher dislocation density. It is therefore important to take accurate account of unsteady gas convection in an analysis of defect formation and optimization of the growth process. Since 2D models usually strongly underestimate the gas cooling of the crystal, the 3D unsteady analysis is necessary for the adequate description of the convective gas effect.

Figure 6. Average heat fluxes across the crystal surface computed with the 3D unsteady (solid lines) and 2D steady (dashed and dash-dotted lines) models.

Figure 7. Average heat fluxes across the encapsulant surface, computed with the 3D unsteady (solid lines) and 2D steady (dashed and dash-dotted lines) models.

References

[1] “Prediction of the melt/crystal interface geometry in liquid encapsulated Czochralski growth of InP bulk crystals”, E.N. Bystrova, V.V. Kalaev, O.V. Smirnova, E.V. Yakovlev, Yu.N. Makarov, J. Crystal Growth, 250 (2003) pp. 189–194.

[2] “3D unsteady analysis of gas turbulent convection during HPLEC InP growth”, E.N. Bystrova, V.V. Kalaev, J. Crystal Growth, 287 (2006) pp. 275–280.