MODELLING OF THE PEDESTAL GROWTH OF SILICON CRYSTALS

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 H_M for the  cylindrical phase of a  large (diameter D_C  =  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 Q_FS, 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,

where q, j, E and A are the magnitudes of induced heat density, current density, electric field and vector potential, respectively; µ₀ = 4π · 10⁻⁷ H/m, ω = 2πf_MF = 6.3 · 10⁵ Hz is angular frequency, $\overrightarrow{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, Q_MF is found as the integral power over all Si domains.

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: H_ind = 9 mm (see Fig. 1, right) for all of considered D_C. In this way, the total induced HF heat Q_HF has been obtained for each calculation. Other integral heat fluxes were user-defined. Q_FS was defined directly, Q_MF 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 Q_MF was set too low, and part of FS crystallized near the internal triple point (ITP) if Q_FS was set too low. In this paper, the calculations with the lowest possible values of additional heat fluxes Q_FS and Q_MF 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 H_M value and FS shape are satisfactory, i.e. do not threaten a  stable growth process, for all of crystal diameters except the  largest D_C = 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 Q_FS (solid lines: from grey to black) or Q_HF (dashed lines: from grey to black). To illustrate this idea, Fig.  5 demonstrates the  influence of Q_FS on the phase boundaries, while D_C = 10 mm and other parameters are constant. The  embedded graph also shows the molten zone height H_Z, 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 D_C are summarized, a correlation between H_Z 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 H_ind, 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 H_M, 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

SILICIO KRISTALŲ AUGINIMO PJEDESTALO BŪDU MODELIAVIMAS