Lumped Parameter Model for Silicon Crystal Growth from Granulate Crucible

In the present paper, a lumped parameter model for the novel Silicon Granulate Crucible (SiGC) method is proposed, which is the basis for a future model‐based control system for the process. The model is analytically deduced based on the hydromechanical, geometrical, and thermal conditions of the process. Experiments are conducted to identify unknown model parameters and to validate the model. The physical consistency of the model is verified using simulation studies and a prediction error of below 2% is reached.


Introduction
Monocrystalline silicon is an essential material for modern technologies today and in the near future. The continuous improvement of the manufacturing processes for monocrystalline silicon is, therefore, key to meet the demand for ever-increasing device performance and reduced production costs.
The novel SiGC method, invented at Leibniz-Institut für Kristallzüchtung (IKZ), has the potential to produce high-quality monocrystalline silicon comparable to the quality of crystals grown by the Floating Zone (FZ) method. [1] Low-cost raw material in form of polysilicon granules from the fluidized bed reactor is used. Thus, the production costs are expected to be much lower compared to the FZ method with low oxygen contents. [2] In the SiGC process, a melt pool contained by a layer of solid silicon ("self-crucible") is generated by inductive heating. The "self-crucible" stabilizes itself in a bed of silicon granules. Then, a crystal is pulled through a central hole in the inductor.
The melt volume during SiGC growth is kept constant as in the continuous Czochralski (CZ) and FZ process. This is achieved by a continuous replenishment and melting of silicon granules.
Contact between the melt and containment is avoided, and graphite heaters can be avoided during the crystal growth.

DOI: 10.1002/crat.202000044
Thus, the SiGC method combines the advantages and avoids the disadvantages of the two industrially-established processes for monocrystalline silicon crystal growth.
Manual control of the SiGC process is complicated because the operator has to stay in a small window of suitable growth parameters. This windows is even smaller for growth of crystals with larger diameter. Thus, the automation of the SiGC method will be a key issue in the further process development and would be the next step towards its industrial application.
The performance of linear control concepts using, for example, Proportional-Integral-Derivative (PID) controllers is limited for crystal growth. [3] As seen in the CZ as well as in the FZ method, crystal growth processes typically show a nonlinear behavior. [4,5] Hence, the parameters of linear controllers have to be determined individually for every phase of the growth process and have to be set up again if the process setup is modified. This expensive procedure can be avoided by using a control concept based on the mathematical model of the system proposed in this paper. Compared to linear control concepts, better performance can be achieved using a model-based control system.
For the CZ process many models and control methods have been proposed. For example, Zheng et al. have recently developed a first-principle model for the industrial CZ silicon production process for predicting crystal radius and crystal growth rate. [6] A sketch of fundamental challenges in the automation of the CZ method and a review of selected control systems can be found in a paper by Winkler. [7] The methods of feed-forward control, feedback control and state estimation are discussed ibidem.
A recent approach to control the CZ process is presented by Winkler and Neubert in a series of papers. [4,[8][9][10] The crystal growth rate and diameter is controlled using a combination of a nonlinear model-based controller and conventional PID controllers. The nonlinear model-based controller is used for parts of the process for which the model is known with sufficient precision. System states not measured in the process are reconstructed utilizing a nonlinear observer.
Werner developed a model-predictive control (MPC) system for the FZ method without additional PID controllers. [5] For this purpose, a measurement system based on visual image processing is used. A nonlinear low-order model is derived considering the geometrical and thermodynamical conditions of the process. Based on this model, an extended Kalman filter (EKF) is www.advancedsciencenews.com www.crt-journal.org implemented. Using these techniques, a high-control precision was achieved.
The crystal growth using the SiGC method is investigated by Menzel using a 2D transient numerical model of the process implemented in Comsol Multiphysics. [2] In this work, suitable process parameters for stabilizing the SiGC growth are identified by simulating the crystal shape under varying growth conditions.
The remainder of this paper is organized as follows: In Section 2 the mathematical model for the SiGC method is derived and a parameter identification is conduced. Section 3 presents simulation results of test cases. Section 4 shows the results from a auto-and cross-validation and conclusions are drawn in Section 5.

Mathematical Model
In the following, a lumped parameter model for the SiGC method is proposed. A set of ordinary differential equations (ODEs) in the form of a nonlinear state-space model is deduced, describing the dynamic behavior of the the crystal radius r c , the crystal slope angle c , the melt filling level h l , the effective power of the inductor P eff and the crystal growth rate v c . The adjustable inputs are the pull rate v p , the silicon replenishment rateṀ in and the generator anode voltage U a . The output is the measurable crystal radius.
Relevant quantities are depicted in Figure 1 and required constants used in the calculations can be seen in Table 1.
The first part of the model deals with the hydromechanical and geometrical conditions of the SiGC method. Three ODEs describing the crystal radius, the crystal slope angle and the melt filling level are derived.

Crystal Radius
An ODE describing the crystal radius r c can be derived from the geometric relations at the triple-point. [12] A rotationally symmetric crystal and a planar crystallization interface is assumed. The crystal growth rate v c is described by the temporal rate of change of the crystal length l with the pull rate v p , the temporal rate of change of the meniscus heightḣ m and the melt filling levelḣ l . The crystal slope angle c is defined as where the melt angle is the angle at the triple-point line between the melt surface and a vertical line. 0 = 11 • is the constant growth angle of silicon. [13]

Melt Filling Level
The melt filling level h l is defined as the depth of the melt pool in the center. It must be kept constant during SiGC growth. Strong fluctuations of the melt filling level due to an improper silicon replenishment rateṀ in will lead to process instabilities and must be avoided by a control system. Thus, calculation of the dynamic behavior of the melt filling level is essential to determine the correct amount of silicon granulate replenishment during the whole process. The ODE describing the melt filling level can be derived from a mass balance. Hence, the time derivatives of the crystal mass, the melt mass, and the meniscus mass are calculated.

Crystal Mass
The rate of change of the crystal massṀ c is calculated aṡ with the density of solid silicon s . A rotationally-symmetric crystal and a planar crystallization interface is assumed.

Melt Mass
The shape of the self-crucible, as well as the melt volume V l it contains, is determined by the temperature field. This field is not modeled in this paper. Consequently, a geometric approximation of the shape of the crucible is derived. It is assumed that the crucible is rotationally symmetric and that the outer radius r l of the free melt surface is constant in time. To model the shape of the self-crucible, the melt height h is approximated as a cubic function depending on the radius r. Unknown constants in this approximation are fitted to experimental data. The melt mass M l and its time derivative are calculated by a volume integral resulting in This result is verified by comparison to a lateral cut of a solidified residual melt of a SiGC crystal growth experiment.

Meniscus Mass
The meniscus mass M m depends on the shape of the meniscus. The shape of the meniscus is governed mainly by the balance of surface tension and hydrostatic force. This is described by the Young-Laplace equation, neglecting the impact of the electromagnetic force on the free melt surface shape (cf. [14,15]). To avoid the computational cost of a numerical solution, an analytical approximation derived by Boucher [16] is used.
The capillary constant of silicon a is calculated by where is the surface tension of the melt, l is the density of liquid silicon and g is the gravitational constant. The force F m from the meniscus acting on the crystal consists of the vertical component of the surface tension at the triple-point line and the hydrostatic pressure drop caused by the melt being elevated over the free melt level in the crucible. [17,18] Using Newton's second law in scalar form (F = mg), an approximation for the meniscus mass M m and its time derivative is calculated utilizing Equation (8). It has been shown by Johansen [19] that this approximation is reasonable by comparing its results to the solutions from solving the Young-Laplace equation.

Mass Balance
The mass balance of the process is written as follows, assuming that no silicon can flow off and all replenished silicon is molten in the melt pool Substituting the left-hand side of the mass balance by the Equations (5), (4) and (9) and solving forḣ l yields the ODE needed to describe the melt filling level.

Crystal Slope Angle
The ODE describing the crystal slope anglėc is derived by substituting the right-hand side of Equation (2) by the Equations (6) and (11) and solving foṙc.

Thermal Model
In the following, ODEs describing the crystal growth rate and the effective power of the inductor based on the thermal behavior of the process are presented. Unknown model parameters Θ i are introduced.
The primary heat source in the process is a one-turn pancakeshaped inductor. A high frequency (f ≈ 2 MHz) of the inductor current is used. The effective power of the inductor P eff is not known in SiGC growth experiments. Hence, it is estimated from the adjustable generator anode voltage U a using a first-order delay elemenṫ where the gain factor is set to K p = 1 W∕V. The relevant heat fluxes considered in the model are shown in Figure 2. Heat sources of the crystal, assumed to be of high relevance for the dynamics of the thermal model, are: • heat conduction from the melt P c , • latent heat released at the crystallization interface P c,l , Cryst. Res. Technol. 2020, 55, 2000044 www.advancedsciencenews.com www.crt-journal.org • induction into the crystal P c,ind .
Hence, a quasi-steady power balance for the crystal can be written as (cf. [5]) where P c,loss is the power loss of the crystal. The dependence on the heat convection is neglected in this model, since it depends to a high degree on the rotational rates of the crucible and the crystal which are not considered. The heat balance for the melt and the heat loss through the porous silicon were not explicitly considered in the model. However, the influence of the melt on the heat balance in the crystal is effectively described by the P c .

Crystallization
A silicon mass M releases an amount of heat Q l , when it solidifies: Where q 0 is the latent heat of fusion of silicon. The power released due to crystallization P c,l is calculated as the time derivative of Q l .
The mass flowṀ is the rate of change of the crystal massṀ c , which is calculated by Equation (4). By substituting Equation (4) into Equation (16), P c,l is expressed as

Power Loss
The power loss of the crystal P c,loss can be analytically approximated as shown by Werner, [5] based on an analytical equation for the surface temperature of the crystal derived by Billig. [20] The power loss of the crystal is calculated as P c,loss = Θ 0 lost r 3 2 c .
The power loss constant lost is calculated by where is the Stefan-Boltzmann constant, T 0 is the melting point of silicon and 0 is the emissivity of solid silicon at T 0 .

Heat Conduction
The power introduced into the crystal by heat conduction from the melt P c is assumed to increase with higher effective power of the inductor P eff or larger crystal radius. Moreover, it is assumed that less power is introduced into the crystal if the meniscus height h m is larger than the equilibrium meniscus height h m,eq = h m (r c , c = 0). This is caused by a smaller temperature gradient at the crystallization interface. For this purpose, the variable Δh m = h m − h m,eq is introduced. Hence, P c is calculated as follows where U a,0 is a reference value for the generator anode voltage.

Induction
The induced power into the crystal P c,ind is assumed to depend on the crystal radius and the effective power of the inductor and is described as:

Crystal Growth Rate
The power released due to crystallization P c,l depends on the crystal growth rate v c (see Equation (17)). Therefore, the power balance for the crystal (Equation (14)) can be solved for v c after substituting the right-hand side by the Equations (18), (20), (17), and (21). The ODE describing the crystal growth rate is obtained by taking the time derivativė

Length-Dependent Mathematical Model
A length-dependent model is derived from the time-dependent model for more accurate parameter identification. The crystal slope angle c and the crystal growth rate v c can be calculated in the length-dependent domain and these reference values are used for the identification of the unknown model parameters. To www.advancedsciencenews.com www.crt-journal.org convert the model to a length-dependent domain, the ODEs are divided by the crystal growth rate v c (cf. [4]) The length-dependent ODEs describing the crystal radiuṡr c and the crystal slope anglėc are calculated as: Both equations are not explicitly dependent on any of the four remaining ODEs. Thus, the model is reduced to these two ODEs for the following derivation.
Flatness: Following the definition given in ref. [21], the reduced model (Equation (25)) is differentially flat. The flat output y flat = r c can be calculated from the reduced states x reduced = [r c , c ]. The new input v z and the reduced state vector x reduced can be calculated from y flat and its first and second derivative without integration.
Using this method, c can be calculated from Equation (26) and the crystal growth rate v c by solving Equation (24) for v c since the pull rate v p and the replenishment rateṀ in are known from experimental data and the crystal radius r c is measured.
In the derivation of the Equations (26) and (24), only the hydromechanical-geometrical model is utilized. This is advantageous since this part of the model is very accurate and deduced without introducing model parameters. Thus, it is reasonable to use the calculated values of c and v c for the identification of unknown model parameters introduced in the thermal model.

Parameter Identification
In order to identify the unknown model parameters Θ and validate their accuracy, two separate SiGC growth experiments were conducted. In the first experiment, steps in the pull rate v p , the generator anode voltage U a and the silicon replenishment ratė M in were realized. The step size, the direction of the steps, and their sequencing was chosen such that the crystal radius was varied within the known stability range of around 30 mm. For this  radius, suitable growth parameters were known from previous experiments and a stable process could be ensured. The grown crystal can be seen in Figure 3.
The whole crystal length, except for the the thin neck and the cone, was used for the parameter identification and the autovalidation. To generate particularly suitable data for the crossvalidation, ramp experiments were conducted in the pull rate and the generator anode voltage in the second growth experiment.
The SiGC growth setup comprises an on-line measurement system based on visual image processing for time-dependent measurement of the crystal radius. Additionally, a method to measure the crystal radius from grown crystals in dependence of the crystal length l was developed. This was done utilizing image analysis techniques.
Using the least squares method, values for the model parameter Θ are found. The simulated crystal radius is compared to the measurement of the crystal radius. The calculated crystal slope angle and crystal growth rate are compared to the reference values calculated from Equations (26) and (24) using the flatness of the reduced model.
The identified values of the model parameters are shown in Table 2.

Simulation of Test Cases
To verify the physical consistency of the mathematical model and to evaluate the influence of each input, step responses are calculated.
In Figure 4, a results for a step in the silicon replenishment rateṀ in are shown.
By decreasing the silicon replenishment rateṀ in , the melt mass and volume decreases. The melt filling level decreases, which increases the meniscus height. This causes the triple-point to move inward, and the crystal radius decreases. At the end of this step response, the crystal radius, the crystal slope angle, and   the meniscus height return to a new equilibrium state with a reduced crystal radius.
In the special case of a permanently depleting melt pool, the mathematical model of the SiGC process shows similar behavior as expected for CZ growth. During CZ growth no raw material is added and, hence, the melt level is constantly decreasing.
For the simulated steps in the pull rate ( Figure 5) and in the generator anode voltage (Figure 6), the rate of change of the silicon replenishment rate was adjusted to match the rate of change of the crystal mass calculated by Equation (4) (Ṁ in =Ṁ c = s r 2 c v c ). Increasing the pull rate yields an increasing meniscus height and a decreasing melt angle and the crystal radius decreases. Adjusting the generator anode voltage changes the temperature field in the process. The temperature in the melt increases slowly, which temporarily reduces the crystal growth rate. After the crystal growth rate reaches its minimal value, the conditions are comparable to the ones during increased pull rate. In both simulations the system returns to an equilibrium state with a reduced crystal radius after the crystal growth rate is equal to the pull rate.
The results seen in all three step responses are physically consistent and agree with experimental observations.

Auto-and Cross-Validation
The model predictions are compared to experimental data from the first crystal growth experiment, which are involved in the parameter identification (auto-validation) and to data from the second experiment, which are not involved in the parameter identification (cross-validation). To quantify the accuracy of the model prediction in comparison with the experimental results, the mean absolute percentage error is calculated for the crystal radius and can be seen in Table 3. In the second experiment, which is used for the cross-validation, only the pull rate and the generator anode voltage were adjusted, thus no error value for the silicon replenishment rate is given. The maximum error of 1.74% for the auto-validation is small. The values for the cross-validation are in similar range as for the auto-validation. This good agreement is seen in the plots of the experimental and simulation data, shown in Figure 7. Hence, the mathematical model provides a general description of the SiGC method is suitable for use in a future model-based control system for the process.

Conclusions
For the first time, a lumped parameter model was derived and validated for SiGC growth. The coupled ODEs have the form of a nonlinear state space model and were analytically deduced from the fundamental physical properties of the SiGC method. The physical consistency of the model was verified and the influences of the inputs were evaluated by simulation studies. The validity of the model was successfully shown by auto-and crossvalidation. A high accuracy, with a prediction error below 2%, was calculated for the crystal radius. This is an important result since the control of the crystal radius is the main goal of a future automation of the SiGC method. A reliable prediction of the crystal radius, which can be obtained by the mathematical model is, thus, essential.
Using the derived lumped parameter model, the behavior of the relevant quantities of the SiGC method: • the crystal radius r c , • the crystal slope angle c , • the melt filling level h l , • the effective power of the inductor P eff , • the crystal growth rate v c , can be predicted with high accuracy for given initial conditions and sequence of inputs.
The model is, thus, suitable to be used as the basis for a future model-based control system for the SiGC process.