# MODELLING OF THE PEDESTAL GROWTH OF SILICON CRYSTALS

K. Surovovs a, A. Kravtsov b, and J. Virbulis a

a University of Latvia, Jelgavas iela 3, 1004 Riga, Latvia

b KEPP EU, 5 Carnikavas Street, 1034 Riga, Latvia
Email: kirils.surovovs@lu.lv

Received 18 December 2020; accepted 13 May 2021

The pedestal method is an alternative to the well-known floating zone method, both of which are performed with high-frequency electromagnetic heating. Unlike the floating zone method, in the pedestal method a single crystal is pulled upwards from the melt. It allows one to lower feed rod quality requirements and simplify the process control due to the absence of open melting front. As the pedestal method has not been widely used in industry for silicon crystals, its development requires extensive numerical modelling. The present work describes application of the previously created mathematical model for crystals with diameters higher than it is currently possible in the experimental setup, as well as for the cone growth phase. Supplementary free surface heating, that prevents melt centre freezing during the seeding phase, has been added at the beginning of cone phase. After multiple sets of simulations, an optimal scheme of heating control for cone growth was proposed.

Keywords: single silicon crystal growth, numerical modelling, pedestal method

## 1. Introduction

### 1.1. Silicon crystals

Single silicon crystals have been widely used in industry for decades and are continuing to be demanded despite advances in research of other semiconductors. Resistivity of the  grown crystal can be controlled by dopants that are added during the crystal growth. To enable precise resistivity control and successful operation of produced devices, the content of other impurities should be reduced to the lowest possible level.

To measure the impurity level in polycrystalline rods (e.g. with Fourier-transform infrared spectroscopy) they must be melted and recrystallized as single crystals. To ensure that the addition of dopants is as low as possible a crucible-free method, that eliminates liquid silicon contact with any other material, should be used.

### 1.2. Pedestal method

The pedestal method (PM), proposed in 1950s by Dash [1], is an alternative to the  well-known floating zone (FZ) method. In both methods, polycrystalline silicon rods are being melted by using a  high-frequency (HF) electromagnetic one-turn inductor. The  main difference between the two methods is that in PM the upper surface of polycrystalline rod is molten, thus the melt is situated on a  silicon ‘pedestal’, and the  single crystal is being pulled upwards from the melt (see Fig. 1, right). A complicated feature of the FZ method – an open melting front, i.e. a lower part of a feed rod where the silicon is being melted and flows downwards – is absent in PM. Therefore, feed rod quality requirements for PM can be lower (and the process control can be simpler) than in the FZ method [2].

PM has a  limitation on an inductor shape in comparison with FZ that prohibits ‘needle-eye’ inductors: the diameter of the grown crystal cannot exceed the  inner diameter of the  inductor. It happens because in the cylindrical stage of crystal growth (when the shapes of the molten zone and the grown crystal diameter do not change) 11° meniscus angle (the angle between the free melt surface and the crystal side surface) is required [3], and the diameter of the molten zone cannot be smaller than the crystal diameter. Due to a large inductor diameter, freezing of the molten zone in the centre becomes more likely when the  crystal diameter increases, and an additional medium-frequency (MF) inductor needs to be used for pedestal heating (see Fig. 1, left).

### 1.3. Numerical modelling

The experimental results are of limited quality and often do not provide an insight into all principal features of the growth process. To get this kind of information, numerical modelling is widely used. It allows one not only to reduce the number of expensive experiments, but also to supplement existing experiments by describing process features that are very hard to measure. Examples include distributions of various physical fields inside solid, melt and surrounding atmosphere: temperature, velocity, electromagnetic field, phase boundaries, etc. Numerical modelling is particularly useful in the early development stages when parameter ranges of successful experiments are not known yet.

### 1.4. Aims of the study

The used inductor shape was optimized via the gradient method [4] to increase the  melt height HM for the  cylindrical phase of a  large (diameter DC  =  100  mm) crystal growth from the  pedestal with a diameter of 200 mm. However, the experiments demonstrate [5] that at the  seeding phase even the addition of a MF inductor is not enough to ensure that the  free surface (FS) is completely molten (see Fig.  2). This difficulty arises because of large heat losses from the  central part. It can be compensated by additional FS heating QFS, e.g. with infrared lamps (schematically shown in Fig. 1, right).

The present work describes the calculations of phase boundaries for different crystal diameters, performed to find appropriate proportions of HF, MF and additional FS heating during the  cone growth.

## 2. Modelled system

The used parameters are listed in Table 1. For a more detailed description of the growth apparatus see Ref. [5].

The MF inductor was made of a 10 mm copper tube and consisted of 3 turns, located approximately 10 mm from the pedestal side surface, as shown in Fig. 1. To separate MF and HF electromagnetic fields, a copper shield was installed between them.

## 3. Mathematical model

The shape of phase boundaries was obtained in the  axially symmetrical approximation, using the  previously developed software [6], based on the principles described in [7]. Melt flow has not been taken into account. To save computational time the  growth of the  crystal cone was approximated by running quasi-stationary simulation for different cone diameters.

To calculate the induced heat of the HF inductor, the high-frequency approximation, as given in [7], was used due to a relatively small skin layer depth (δ = 1.4  mm for the  frequency of 2.6  MHz). MF induced heat, however, cannot be approximated by surface heat density only, thus the vector potential →A was calculated in all volume domains,

$\nabla \text{\hspace{0.17em}}×\text{\hspace{0.17em}}\nabla \text{\hspace{0.17em}}×\text{\hspace{0.17em}}\stackrel{\to }{A}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{\mu }_{0}\sigma \left(-\text{i}\omega \stackrel{\to }{A}\right)\text{\hspace{0.17em}}+\text{\hspace{0.17em}}{\mu }_{0}\stackrel{\to }{J},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left(1\right)$
$q\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{{j}^{2}}{2\sigma }\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{\sigma {E}^{2}}{2}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{\sigma }{2}\text{\hspace{0.17em}}{\omega }^{2}\text{\hspace{0.17em}}{A}^{2},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left(2\right)$

where q, j, E and A are the magnitudes of induced heat density, current density, electric field and vector potential, respectively; µ0 = 4π · 10−7 H/m, ω = 2πfMF = 6.3 · 105 Hz is angular frequency, $\stackrel{\to }{J}$ is current density in the MF inductor and σ is the Si electric conductivity. Only the  azimuthal component Aφ was taken into account due to axial symmetry (disregarding current suppliers). Boundary conditions were the following: Aφ = 0 on the axis and $\frac{\partial {A}_{\phi }}{\partial n}$ on the outer boundaries of the domain.

MF induced heat calculation is coupled with the calculation of phase boundaries. After the quasistationary shape of phase boundaries has been reached, QMF is found as the integral power over all Si domains.

Table 1. Main parameters of the modelled system
Crystal diameter DC 10–100 mm
Pedestal diameter DP 200 mm
Crystal pulling rate vp 2 mm/min
MF heating power QMF. 8.6–9.4 kW
HF heating power QHF 9.0–10.0 kW
Free surface heating QFS 0–0.45 kW
MF inductor frequency fMF 100 kHz
HF inductor frequency fHF 2.8 MHz
Distance from MF inductor to HF inductor 33 mm
Distance from ETP to HF inductor Hind 9 mm
Distance from ITP to FS heating region 5 mm
HF inductor parameters:
internal diameter 122 mm
distance to ETP 11 mm
cross-section length 40 mm
cross-section inclination angle
thickness 6 mm

## 4. Results

During the  study, the  inductor shape remained constant, and the  inductor current was adjusted using a  proportional-integral-derivative (PID) algorithm to keep the vertical distance to the external triple point (ETP) constant: Hind = 9 mm (see Fig. 1, right) for all of considered DC. In this way, the total induced HF heat QHF has been obtained for each calculation. Other integral heat fluxes were user-defined. QFS was defined directly, QMF was defined indirectly – by defining MF inductor current. Multiple values were tested, and only some of calculations converged. For example, the  melt centre was freezing when QMF was set too low, and part of FS crystallized near the internal triple point (ITP) if QFS was set too low. In this paper, the calculations with the lowest possible values of additional heat fluxes QFS and QMF are summarized, as they are beneficial for cost-effectiveness and equipment design. The optimized change of integral heat fluxes during the cone phase is given in Fig. 3, the corresponding phase boundaries and silicon temperature are presented in Fig. 4.

Both HM value and FS shape are satisfactory, i.e. do not threaten a  stable growth process, for all of crystal diameters except the  largest DC = 100 mm. In the latter case, FS is bulged outwards near the  ETP, thus creating a  high risk of melt spilling. We define α as the  angle between the  FS and the  vertical line, with α  >  0 indicating that the surface goes ‘outwards’ from the ETP. Figure 4 shows that α is increasing when the induced heat in the vicinity of the crystal increases: either QFS (solid lines: from grey to black) or QHF (dashed lines: from grey to black). To illustrate this idea, Fig.  5 demonstrates the  influence of QFS on the phase boundaries, while DC = 10 mm and other parameters are constant. The  embedded graph also shows the molten zone height HZ, defined as a vertical distance between the ETP and the ITP.

From Figs  4 and 5 it can be concluded that the  increase of heat flux near the  crystal shifts the ITP upwards, increasing the zone height, which increases α. This conclusion is affirmed by Fig. 6: when all considered DC are summarized, a correlation between HZ and α is very strong. It means that to mitigate the risk of melt spilling over the ETP, the heat flux should be reduced in the vicinity of ITP.

A possible solution is to increase the  inner diameter of the  HF inductor. However, it would mean that induced heat is concentrated too close to the ETP, and low HF current will be sufficient to maintain the  defined Hind, which will increase the  risk of melt centre freezing. The  risk could be prevented by increasing the  distance between the  ETP and induced heat maximum, e.g. by increasing the  pedestal diameter. It should also be noted that the  currently used inductor shape has been optimized to improve only the  melt height HM, and another target function that includes α may be useful.

## 5. Conclusions

The performed calculations indicate that both MF and HF heating must increase during the  cone phase, in order to a) compensate the  decrease of FS heating and b) compensate radiative heat losses from the crystal. The scheme of this increase is proposed in Fig. 3. The maximal FS heating, which is required for crystal diameters smaller than 20 mm, is less than 500 W.

The risk of melt spilling occurs only for the largest crystal diameter, i.e. in the cylindrical phase. To decrease the  risk, the  free surface angle could be improved by increasing the pedestal diameter or by modifying the target function in the inductor optimization method to take this angle into account.

## Acknowledgements

The authors express their gratitude for the financial support provided by the  European Regional Development Fund Project ‘Development of Competence Centre of Mechanical Engineering’, Contract No.  1.2.1.1/18/A/008, Research No.  4.5 ‘Technological Research and Silicon Production with Diameters up to 100 mm for Low-current and Highpower Microelectronic Solid State Devices’.

## References

[1] W.C. Dash, Silicon crystals free of dislocations, J. Appl. Phys. 29, 736 (1958).

[2] W.C. Dash, Improvements on the pedestal method of growing silicon and germanium crystals, J. Appl. Phys. 31, 736 (1960).

[3] M. Wünscher, A. Lüdge, and H. Riemann, Growth angle and melt meniscus of the RF-heated floating zone in silicon crystal growth, J. Cryst. Growth 314, 43–47 (2011).

[4] K.  Surovovs, A.  Kravtsov, and J.  Virbulis, Optimization of the shape of high-frequency inductor for the  pedestal growth of silicon crystals, Magnetohydrodynamics 55(3), 353–366 (2019).

[5] A.  Kravtsov, K.  Surovovs, and J.  Virbulis, Float zone single crystals for testing rods, pulled under electron beam heating, IOP Conf. Ser. Mater. Sci. Eng. 503, 012022 (2019).

[6] K. Surovovs, M. Plate, and J. Virbulis, Modelling of phase boundaries and melt flow in cruciblefree silicon crystal growth using high-frequency heating, Magnetohydrodynamics 53, 715–722 (2017).

[7] G.  Ratnieks, A.  Muižnieks, and A.  Mühlbauer, Modelling of phase boundaries for large industrial FZ silicon crystal growth with the needle-eye technique, J. Cryst. Growth 255, 227–240 (2003).

#### SILICIO KRISTALŲ AUGINIMO PJEDESTALO BŪDU MODELIAVIMAS

K. Surovovs a, A. Kravtsov b, J. Virbulis a

a Latvijos universitetas, Ryga, Latvija

b KEPP EU, Ryga, Latvija