Smart Design of Cz-Ge Crystal Growth Furnace and Process

Abstract

The aim of this study was to evaluate the potential of the machine learning technique of decision trees to understand the relationships among furnace design, process parameters, crystal quality, and yield in the case of the Czochralski growth of germanium. The ultimate goal was to provide the range of optimal values of 13 input parameters and the ranking of their importance in relation to their impact on three output parameters relevant to process economy and crystal quality. Training data were provided by CFD modelling. The variety of data was ensured by the Design of Experiments method. The results showed that the process parameters, particularly the pulling rate, had a substantially greater impact on the crystal quality and yield than the design parameters of the furnace hot zone. Of the latter, only the crucible size, the axial position of the side heater, and the material properties of the radiation shield were relevant.

1. Introduction

Today, a century after its invention, Czochralski (Cz) growth is the dominant method in industry to produce large crystalline materials for modern electronic and photonic technologies [1,2,3]. Out of many parameters influencing the crystal quality and yield, the furnace design (e.g., geometry and the material properties of the hot zone parts) and process parameters (e.g., pulling and rotational rates, heating power, etc.) are essential. Their optimization has been the subject of extensive experimental and numerical research, especially in the case of the Cz growth of silicon (Cz-Si), e.g., [3,4,5,6,7,8,9], followed by germanium, gallium arsenide, and various oxides [10,11,12,13].

The furnace design studies focused mainly on optimizing the geometry of the radiation shield [4,11,14,15,16,17,18] and rarely on the position of the side heater [14], the geometry of the cooler [16], the flow guide [4], and the insulation materials [4,17]. Concerning process parameters, numerous papers were devoted to the study of the influences of the crystal and crucible rotational rates on the interface shape [9,19] and the limitations of the growth rates [20,21]. In magnetically driven Cz-Si growth, optimization of the magnetic parameters with respect to oxygen and carbon transport was a topic of numerous studies, e.g., [22,23,24,25,26].The sensitivity of the Cz-Ge process dynamics to the pulling rate and heating power at the beginning and towards the end of the growth was investigated in [10].

Computational fluid dynamics (CFD) is a traditional numerical tool for solving the governing differential equations describing the transport phenomena during crystal growth that helped understand the crucial process steps and factors determining the crystal growth. The main challenge of this approach is the fact that CFD simulations are laborious, expensive, and time consuming without the ability to generalize.

One of the most exciting cutting-edge tools that entered the field of crystal growth over the last decade is machine learning (ML). It has a great potential to revolutionize the optimization of the crystal growth processes and equipment, process control, and interpretability of the results in general [27].

Based on the added value achieved by combining ML and crystal growth research fields, we aimed to assess the potential of a selected interpretable ML technique for understanding the relationships among furnace design, process parameters, crystal quality, and yield in the case of Cz growth of germanium (Cz-Ge) to improve on the latter two. This study was motivated by the fact that Ge, a pioneering crystalline material in the history of electronics, is making a comeback as a next-generation semiconductor for applications in infrared optics, gamma-radiation detection, and photovoltaics [28,29,30,31]. The ultimate success in this endeavor depends on the quality and yield of the Ge crystals and their relationship to actual Si values [8].

One of the first attempts to develop the geometry of a bulk crystal growth furnace with ML support goes back to the work of Yu et al. [32], which was focused on the top-seeded solution growth of silicon carbide (TSSiC). The authors derived an artificial neural network (ANN) from 500 axisymmetric CFD simulations varying seven geometrical furnace parameters and optimizing the geometry using a genetic algorithm. Further ML studies were devoted to the applications of ANN and convolutional neural networks (CNN) to optimize process parameters in Cz-Si, Cz-YAG, DS-Si, and TSSiC growth [33,34,35,36,37].

The ANN and CNN are the most popular ML methods, characterized by a high accuracy of predictions but also a black-box nature, long computational times, and the demand for a large amount of training data. Our ML focus was on using the decision trees (DTs) method [38,39], which is suitable for applications with a small amount of data and high demands on the interpretability of the results (i.e., on understanding cause–effect relationships), which are common in crystal growth. Recently, we successfully applied regression and classification trees to study the influence of process parameters on interface deflection in VGF-GaAs growth [40].

In this study, Cz-Ge furnace and process design considerations included the hot zone geometry parameters, radiation shield material, heating power, rotation, and pulling rates.

Altogether, 13 input and 3 output parameters were considered. In total, 77 training data sets were provided by CFD modelling.

The pros and cons of the here proposed data-driven DT approach for Cz-Ge furnace and process design will be discussed and results compared with the common CFD and data-driven ANN approach.

2. Models and Methodology

2.1. Generation of Training Data by CFD Modelling

A sketch of the Cz-Ge furnace geometry with a definition of the main geometric and process parameters used in the simulations is given in Figure 1. The furnace was equipped with two graphite resistance heaters positioned on the side and bottom of the crucible. The quartz crucible was charged with 6.7 kg of Ge. The radiation shield was made of different materials and varied in length but had a constant slope. Its geometry guaranteed the visibility of the Ge triple point.

The transport phenomena taking place during the Cz-Ge growth were described by a common two-dimensional axisymmetric model. Governing equations for CFD modelling included equations of continuity, Navier-Stokes with the Boussinesq approximation, and energy Equations (1)–(5).

$$\frac{\partial \rho }{\partial t}+\nabla \left(\rho ·\overrightarrow{u}\right)=0$$

(1)

$$\frac{\partial \left(\rho \overrightarrow{u}\right)}{\partial t}+\left(\overrightarrow{u}·\nabla \right)\rho \overrightarrow{u}=-\nabla p+\nabla \tau +(\rho -{\rho }_0)\overrightarrow{g}$$

(2)

$$\frac{\partial \left(\rho c_{p}T\right)}{\partial t}+\nabla ·\left(\rho c_p\overrightarrow{u}T\right)=\nabla ·\left(\lambda \nabla T\right)$$

(3)

The Stefan and isothermal conditions must be fulfilled at the melt–crystal interface:

$$\left[{\left(-{\lambda }_l\nabla T\right)}_l+{\left({\lambda }_s\nabla T\right)}_s\right]·{\overrightarrow{n}}_{s,l}=∆H_{s}v$$

(4)

$$T=T_m$$

(5)

Details and theoretical backgrounds can be found elsewhere [41].

Due to the small size of the crucibles and crystals, Reynolds numbers were expected to be below the range typical for the turbulent flow. Against this background, melt convection was described by the laminar flow model.

The Ge material properties used in this study are given in

Table 1. Material properties of Ge melt and crystal.

Property Unit	Tm (K)	ρ (kg m−3)	μ (Pas)	λ (W m−1K−1)	cp (J kg−1K−1)	α (K−1)	ε (-)	Pr (-)
Ge melt	1211	5670	7.42 × 10−4	39	394	1.11 × 10−4	0.5	0.008
Ge crystal		5323	-	17.3	404		0.5

. The crystals were grown in Ar atmosphere under atmospheric pressure. The CFD simulations were performed using the commercial code CGSim.

In order to obtain maximum information from a minimum number of CFD simulations, i.e., to cover relationships among 13 inputs and 3 outputs with good quality, a Design of Experiments (DoE) method [42] was used. If all input variables were varied at three levels, the theoretical number of CFD simulations (i.e., virtual experiments) based on a full factorial design would be 3¹³ = 1,594,323, which is impossible to perform. In this study, we applied a D-optimal design for eight variables at three levels. However, nonsensical geometrical combinations and constrained process variables were discarded from DoE considerations.

For the generation of data sets for ML simulations, 77 combinations of seven geometric, one radiation shield material choice, and five crystal growth process parameters were used as input parameters for the CFD simulations. From the obtained axisymmetric CFD results, the interface deflection, crystal growth rate, and ratio of growth rate to the vertical temperature gradient (v/G) were extracted as output parameters for the ML training database. The interface deflection y1 was selected as the main measure of the crystal quality due to its well-known correlation with the development of thermal stress-induced crystal dislocations [43].

The output v/G was included as an additional measure of crystal quality based on the following considerations. The so-called Voronkov criterion is one of the key parameters that defines the critical value of v/G for the growth of defect-free perfect silicon crystals for device processing [44]. For v/G > 1.34 × 10⁻³ cm²/K min, the silicon crystal is vacancy rich, while for v/G < 1.34 × 10⁻³ cm²/K min, the crystal is self-interstitial rich. Recent studies showed that application of the Voronkov criterion to a Cz-Ge single crystal is also possible when Ge is doped with Ga and Sb [45,46,47,48]. Otherwise, the vacancy is always the dominant intrinsic point defect for Cz-Ge [48]. Finally, we selected the crystal growth rate as a crucial parameter for the economics of the process.

The common optimization goal in industrial crystal growth is to find ranges of input parameters in which the economy and crystal quality criteria are met simultaneously. In our study, this multi-objective optimization goal corresponded to the cases where s/l interface was flat, the crystal growth rate was maximal, and the Voronkov criterion was satisfied. In other words, we looked for the minimum of the objective (fitness) function y₁₂₃ that was defined as a combined normalized output:

$$y_{123}=\left\{\left|y_1\right|\right\}+\left\{-y_2\right\}+\{y_3-1.34\times {10}^{-3}\}$$

(6)

The data were pre-processed by normalization rescaling so that they ended up ranging between 0 and 1:

$$\left\{y_i\right\}=\frac{y_i-y_{i,min}}{y_{i,max}-y_{i,min}}$$

(7)

The fitness function y₁₂₃ can also be defined in different ways, emphasizing the influence of certain outputs by adding weighting factors to Equation (6), as suggested in the literature [49]. In this study, we assumed that all outputs were equally important.

The interface deflection y₁ was measured at the melt symmetry axis with respect to the three-phase junction (melt/crystal/gas) and varied between detrimental concave (y₁ > 0) and favorable slightly convex (y₁ < 0), as shown in Figure 1c. Process parameters were selected from the typical range of values given in the literature [28], e.g., pulling rate, crystal, and crucible rotational rates were varied in the intervals 0.16–1.5 mm/min, 10–60 rpm and, −1 to −5 rpm, respectively. Although the radiation shield can be made from multiple materials in a stacked arrangement, for simplicity we used a single material shield made of either quartz, ceramics, or graphite. All material properties of the shield materials were taken into account in the CFD simulations. In the ML simulations, only the emissivity of the radiation shield was considered as an input parameter.

2.2. Correlation Coefficients

Because of its simplicity, we used correlation coefficients as a starting point for data analysis. The Pearson correlation coefficient r of any input x_i and output y_j pairs of random variables is a measure of their linear dependence [50]. If each variable has N observations, then the Pearson correlation coefficient is defined as

$$r\left(x_i,y_j\right)=\frac{1}{N}{\displaystyle \sum }_{n=1}^N\left(\frac{x_{i,n}-{\overline{x}}_i}{{\sigma }_{x_i}}\right)\left(\frac{y_{j,n}-{\overline{y}}_j}{{\sigma }_{y_j}}\right)$$

(8)

where ${\overline{x}}_i$ and ${\sigma }_{x_i}$ are the empirical mean and standard deviation of x_i, respectively, and ${\overline{y}}_j$ and ${\sigma }_{y_j}$ are the empirical mean and standard deviation of y_j.

Two variables, x_i and y_j, are positively correlated (r > 0) if their values change in the same direction, negatively correlated (r < 0) if they change in the opposite direction, or may exhibit zero correlation (r = 0) if there is no relationship between the change of the variables. When interpreting the r = 0 results, the limitations of this simple approach must be taken into account, e.g., results may indicate zero correlation where non-linear dependences exist. Additionally, an observed correlation between the variables does not imply causation.

2.3. Regression Trees

The dependency between input and output parameters was modelled using regression trees (RTs), a non-linear regression model suitable for small data sets and high interpretability requirements. Regression trees are a type of decision trees where the output variables are real numbers [39].

The method enables us to analyze how the combination of input variables affects the outputs and to predict them individually. It has a tree structure, starting from the topmost root node with the initial data set, undergoing recursive binary splitting down in each node until the terminal node is reached that either is in a given maximal depth or has at most a given minimum size. The nodes represent the input variables: each branch is a decision and both nodes below a branch are the more homogeneous subsets of the data set above the branch. The splits are based on the sum of squared errors (SSE) with respect to the means of the y_j in the two subsets, $S_1 {\mathrm{and}S}_2,$ forming the split

$$SSE\left(S_1,S_2\right)={\displaystyle \sum }_{x_i\in S_1}{\left(y_i-\frac{1}{\left|S_1\right|}{\displaystyle \sum }_{x_i\in S_1}{y}_i\right)}^2+{\displaystyle \sum }_{x_i\in S_2}{\left(y_i-\frac{1}{\left|S_2\right|}{\displaystyle \sum }_{x_i\in S_2}{y}_i\right)}^2$$

(9)

From all possible splits of the data set S above the branch, a split $S_1^{*},S_2^{*}$ was chosen that led to the minimal sum of squared errors.

$$SSE\left(S_1{}^{*},S_2{}^{*}\right)=\min \left\{SSE\left(S_1,S_2\right);\left(S_1,S_2\right) \mathrm{is}\mathrm{some}\mathrm{split}\mathrm{of} S\right\}$$

(10)

The fitted model can be used to predict y_j values as a function of x_i variables. In this study, the maximum tree depth was set at 10.

From the regression trees for the crystal growth rate, the interface deflection, and the v/G ratio, the importance of each input variable (also known as feature importance) was also derived. Feature importance is the sum of the mean square error (RMSE) reductions across all those interior nodes in which this input variable resulted in the highest RMSE reduction in the node.

3. Results and Discussion

3.1. CFD Modelling

The CFD results provided 77 data sets in the form of 16-dimensional vectors, each with 13 input and 3 output parameters. All data are presented in parallel coordinates in Figure 2. Each line in the diagram corresponds to one data tuple (x₁ … x₁₃, y₁ … y₃). The generated database was used for ML training and analysis. The line color in Figure 2a,b matches the value of the interface deflection y₁ and v/G ratio y₃, respectively. In both cases, the visual inspection of the results indicated the strong influence of the pulling rate x₅, as known from experimental investigations, i.e., the higher the pulling rate x₅, the lower the crystal quality due to higher concavity of the interface (y₁ > 0) and the v/G value (y₃). Similarly, as the crystal radius x₁ increases, so does the interface concavity y₁ increase, which is the well-known challenge of upscaling.

Examples of axisymmetric quasi-steady state CFD simulation results for various furnace designs and growth recipes in a form of temperature and melt velocity distributions in the Cz-furnace are shown in Figure 3. The melt flow patterns reflected the typical superposition of buoyancy-driven convection and forced convection due to crystal and crucible rotation. Interface deflection varied among concave (Figure 3a), convex (Figure 3b), and nearly flat (Figure 3c) depending on the used growth recipe, hot zone design, and materials’ allocation.

3.2. Data Mining and Machine Learning

The resulting values for the correlation coefficients for all inputs and outputs are shown in Figure 4 and

Table A1. Correlation coefficients for input and output variables x₁–x₁₃, y₁–y₃.

	x1	x2	x3	x4	x5	x6	x7	x8	x9	x10	x11	x12	x13
x1	1.00	0.77	−0.24	−0.20	−0.19	0.10	0.00	0.07	0.21	0.16	−0.10	0.09	−0.95
x2	0.77	1.00	−0.20	−0.17	−0.48	−0.01	−0.01	0.14	0.38	0.21	−0.05	−0.55	−0.92
x3	−0.24	−0.20	1.00	0.05	0.00	−0.02	0.01	−0.06	0.23	−0.02	0.02	0.00	0.22
x4	−0.20	−0.17	0.05	1.00	−0.06	0.03	0.01	0.04	−0.02	0.25	0.30	0.02	0.20
x5	−0.19	−0.48	0.00	−0.06	1.00	0.13	0.07	−0.15	−0.26	−0.28	−0.21	0.49	0.33
x6	0.10	−0.01	−0.02	0.03	0.13	1.00	−0.04	0.09	0.08	−0.14	−0.03	0.15	−0.06
x7	0.00	−0.01	0.01	0.01	0.07	−0.04	1.00	−0.07	0.15	−0.02	0.00	0.02	0.02
x8	0.07	0.14	−0.06	0.04	−0.15	0.09	−0.07	1.00	−0.75	−0.08	0.07	−0.11	−0.09
x9	0.21	0.38	0.23	−0.02	−0.26	0.08	0.15	−0.75	1.00	0.15	−0.07	−0.29	−0.31
x10	0.16	0.21	−0.02	0.25	−0.28	−0.14	−0.02	−0.08	0.15	1.00	0.57	−0.13	−0.21
x11	−0.10	−0.05	0.02	0.30	−0.21	−0.03	0.00	0.07	−0.07	0.57	1.00	−0.07	0.07
x12	0.09	−0.55	0.00	0.02	0.49	0.15	0.02	−0.11	−0.29	−0.13	−0.07	1.00	0.20
x13	−0.95	−0.92	0.22	0.20	0.33	−0.06	0.02	−0.09	−0.31	−0.21	0.07	0.20	1.00
y1	0.30	−0.01	0.00	0.12	0.59	0.20	0.25	−0.29	0.14	0.09	−0.01	0.38	−0.19
y2	−0.19	−0.48	0.00	−0.06	1.00	0.13	0.07	−0.15	−0.26	−0.28	−0.21	0.49	0.34
y3	−0.06	−0.41	0.03	−0.01	0.98	0.16	0.08	−0.13	−0.24	−0.25	−0.20	0.55	0.22

.

The strongest linear correlation between interface deflection y₁ and one of the inputs was observed for pulling rate x₅. Their correlation coefficient with value r_x5,y1 = 0.5877 showed medium-strong positive correlation. This result is in agreement with the experimental observations that higher pulling rates generate more latent heat at the crystallization front and finally more concave interface. Please note that in this study convex interface deflection had a negative value (y₁ < 0) and a negative correlation was beneficial for the crystal quality. The strongest linear correlation between the growth rate y₂ and one of the inputs was again observed for the pulling rate x₅, with a correlation coefficient of r_x5,y2 = 0.9999. This finding is supported by crystal growth experiments and is one of the reasons why pulling rates and crystal growth rates are often treated as interchangeable terms. Finally, the strongest positive linear correlation between y₃ and other inputs was observed for x₅ (the pulling rate) and for x₁₂ (distance from the radiation shield to the melt, i.e., melt gap) with correlation coefficients of r_x5,y3 = 0.9781 and r_x12,y3 = 0.5517, respectively. Similar results were reported in the literature for Cz-Si growth [5,7]. Additionally, the strong linear correlations were observed between inputs x₁ (crystal radius) and x₂ (crucible radius) from one side and x₁₃ (melt height) from the other side. The large negative values r_x1,x13 = −0.9512 and r_x2,x13 = −0.9213 were a consequence of the assumption of constant Ge loading in this study.

The correlation coefficients provided the first simple estimate of the linear dependencies between the examined parameters. However, many dependencies in crystal growth are strongly non-linear and need to be further investigated. A viable solution for studying complex relationships in small data sets is provided by regression trees.

The regression trees for outputs y₁–y_3, y_1,2,3 are given in Figure 5, Figure 6, Figure 7 and Figure 8 and summarized in

Table 2. The most decisive inputs and their optimal values for the Cz-Ge growth with flat or slightly convex s/l interface (y₁ ≤ 0 mm), derived from RT analysis. The data ranges correspond to the white nodes in Figure 5.

Mean	Decisive Inputs
y1	x4	x5	x9	x10	x13
−1 −1.15	0.55 < ∗ < 0.825 <0.55	<22.8 <22.8	1.25 < ∗ < 9.4 1.25 < ∗ < 9.4	12.5 < ∗ < 52.5 12.5 < ∗ < 52.5	>41.5 >41.5

,

Table 3. The most decisive inputs and their optimal values for the Cz-Ge growth with high growth rate y₂, derived from RT analysis. The data ranges correspond to the white nodes in Figure 6.

Mean	Decisive Inputs
y2	x4	x5	x7
89.6 89.2	>0.575 <0.575	>82 >82	<104 <104

,

Table 4. The most decisive inputs and their optimal values for the favorable y₃ (cm²/K min) value corresponding to the Voronkov criteria, derived from RT analysis. The data ranges correspond to the white nodes in Figure 7.

Mean	Decisive Inputs
y3	x2	x5	x9	x10	x11
0.00126 0.00141	>3.18 <3.18	12.5 < ∗ < 22.8 12.5 < ∗ < 22.8	<11.3 <11.3	>15 >15	<−2.5 <−2.5

and

Table 5. The most decisive inputs and their optimal values for the favorable y₁₂₃ value, derived from RT analysis. The data ranges correspond to the white nodes in Figure 8.

Mean	Decisive Inputs
y123	x2	x4	x7
0.858	>2.73	>0.55	>121

. The value in the root node at the top of each regression tree represents the mean of the examined output y_i in all the data in the database. The values in interior nodes and leaf nodes correspond to the mean value of output y_i in the set of data remaining over the node after the last splitting. The root mean square error (RMSE) for each split is given in

Table A2. The RMSE values for each split in regression trees for y₁, y₂, y₃ and y₁₂₃. The nodes are numbered from the top of the tree down and from the left to the right at each tree depth.

	y1	y2	y3	y123		y1	y3		y1	y3
Node	RMSE	RMSE	RMSE	RMSE	Node	RMSE	RMSE	Node	RMSE	RMSE
1	7.59 × 1000	2.87 × 1001	2.63 × 10−03	2.52 × 10−01	42	1.80 × 10−01	1.09 × 10−04	83	0.00 × 1000	0.00 × 1000
2	6.47 × 1000	1.18 × 1001	1.11 × 10−03	6.62 × 10−02	43	5.89 × 10−01	0.00 × 1000	84	4.44 × 10−16	0.00 × 1000
3	6.24 × 1000	1.44 × 1001	1.63 × 10−03	2.08 × 10−01	44	6.24 × 10−01	9.88 × 10−05	85	0.00 × 1000	0.00 × 1000
4	3.47 × 1000	3.91 × 1000	5.92 × 10−04	6.26 × 10−02	45	6.24 × 10−01	4.34 × 10−19	86	8.88 × 10−16	0.00 × 1000
5	4.62 × 1000	7.39 × 1000	2.64 × 10−04	5.67 × 10−02	46	3.55 × 10−15	1.81 × 10−04	87	0.00 × 1000	0.00 × 1000
6	1.87 × 1000	7.72 × 1000	1.06 × 10−03	1.37 × 10−01	47	4.08 × 10−01	8.67 × 10−19	88	6.10 × 10−01	0.00 × 1000
7	5.15 × 1000	1.39 × 10−01	8.89 × 10−04	1.68 × 10−01	48	2.18 × 1000	6.78 × 10−05	89	3.75 × 10−01	6.90 × 10−05
8	0.00 × 1000	2.82 × 1000	3.58 × 10−04	5.55 × 10−02	49	1.17 × 1000	0.00 × 1000	90	2.16 × 10−01	2.06 × 10−05
9	2.15 × 1000	1.58 × 1000	2.65 × 10−04	4.27 × 10−02	50	0.00 × 1000	0.00 × 1000	91	1.00 × 10−01	5.37 × 10−05
10	2.92 × 1000	1.41 × 1000	2.23 × 10−04	4.97 × 10−02	51	0.00 × 1000	2.49 × 10−04	92	0.00 × 1000	6.60 × 10−05
11	2.33 × 1000	9.17 × 10−01	1.95 × 10−04	6.00 × 10−02	52	5.00 × 10−01	8.38 × 10−05	93	0.00 × 1000	1.78 × 10−05
12	1.67 × 1000	2.37 × 1000	6.49 × 10−04	1.41 × 10−01	53	8.88 × 10−16	5.03 × 10−05	94	4.44 × 10−16	3.45 × 10−05
13	1.15 × 1000	1.61 × 1000	8.15 × 10−04	1.08 × 10−01	54	0.00 × 1000	2.96 × 10−04	95	0.00 × 1000	2.64 × 10−05
14	1.61 × 1000	2.05 × 10−01	5.17 × 10−04	7.30 × 10−02	55	0.00 × 1000	1.75 × 10−04	96	2.22 × 10−16	7.01 × 10−05
15	2.64 × 1000	0.00 × 1000	5.49 × 10−04	1.19 × 10−01	56	2.50 × 10−01	1.19 × 10−04	97	0.00 × 1000	7.62 × 10−05
16	1.12 × 1000	0.00 × 1000	2.94 × 10−04	3.60 × 10−02	57	2.22 × 10−16	0.00 × 1000	98	8.88 × 10−16	0.00 × 1000
17	2.05 × 1000	0.00 × 1000	2.64 × 10−04	4.16 × 10−02	58	0.00 × 1000	9.60 × 10−05	99	0.00 × 1000	3.40 × 10−06
18	2.03 × 1000	1.12 × 10−01	1.55 × 10−04	8.84 × 10−02	59	7.07 × 10−01	9.02 × 10−05	100	0.00 × 1000	0.00 × 1000
19	2.12 × 1000	0.00 × 1000	1.14 × 10−04	7.98 × 10−02	60	1.14 × 1000	0.00 × 1000	101	0.00 × 1000	1.65 × 10−05
20	7.07 × 10−01	6.28 × 10−01	2.68 × 10−05	1.03 × 10−01	61	6.62 × 10−01	4.33 × 10−05	102	0.00 × 1000	2.68 × 10−05
21	1.15 × 1000	3.55 × 10−15	4.34 × 10−19	1.55 × 10−02	62	0.00 × 1000	1.32 × 10−04	103	0.00 × 1000	8.12 × 10−05
22	6.61 × 10−01	7.11 × 10−15	1.54 × 10−04	2.46 × 10−02	63	1.78 × 10−15	1.50 × 10−05	104	8.16 × 10−01	6.55 × 10−05
23	0.00 × 1000	0.00 × 1000	0.00 × 1000	6.93 × 10−03	64	0.00 × 1000	0.00 × 1000	105	5.00 × 10−01	0.00 × 1000
24	8.54 × 10−01	0.00 × 1000	3.62 × 10−04	4.81 × 10−02	65	0.00 × 1000	0.00 × 1000	106	0.00 × 1000	0.00 × 1000
25	0.00 × 1000	6.20 × 10−01	2.73 × 10−04	2.87 × 10−02	66	0.00 × 1000	0.00 × 1000	107	0.00 × 1000	0.00 × 1000
26	0.00 × 1000	0.00 × 1000	1.72 × 10−04	4.37 × 10−02	67	0.00 × 1000	0.00 × 1000	108	0.00 × 1000	0.00 × 1000
27	1.14 × 1000	1.26 × 10−01	5.86 × 10−04	2.81 × 10−02	68	1.00 × 10−01	0.00 × 1000	109	4.90 × 10−01	0.00 × 1000
28	2.15 × 1000	9.43 × 10−02	0.00 × 1000	6.82 × 10−02	69	1.00 × 10−01	4.90 × 10−05	110	0.00 × 1000	3.06 × 10−05
29	9.43 × 10−01	0.00 × 1000	2.81 × 10−04	2.34 × 10−02	70	2.00 × 10−01	8.04 × 10−05	111	0.00 × 1000	0.00 × 1000
30	0.00 × 1000	4.33 × 10−02	3.62 × 10−04	9.16 × 10−02	71	0.00 × 1000	4.34 × 10−19	112	8.16 × 10−02	3.86 × 10−05
31	8.16 × 10−01	7.45 × 10−02	3.39 × 10−04	8.80 × 10−02	72	2.50 × 10−01	2.02 × 10−06	113	0.00 × 1000	1.08 × 10−19
32	9.43 × 10−01	7.11 × 10−15	0.00 × 1000	1.36 × 10−02	73	0.00 × 1000	5.84 × 10−06	114	0.00 × 1000	0.00 × 1000
33	8.88 × 10−16	0.00 × 1000	1.60 × 10−04	1.72 × 10−02	74	2.50 × 10−01	8.67 × 10−19	115	0.00 × 1000	0.00 × 1000
34	6.24 × 10−01	1.36 × 10−01	5.59 × 10−05	6.25 × 10−02	75	0.00 × 1000	1.39 × 10−04	116	5.00 × 10−02	0.00 × 1000
35	1.06 × 1000	5.00 × 10−01	1.19 × 10−04	1.32 × 10−01	76	3.55 × 10−15	0.00 × 1000	117	0.00 × 1000	4.82 × 10−06
36	0.00 × 1000	0.00 × 1000	5.36 × 10−05	3.71 × 10−02	77	2.50 × 10−01	8.67 × 10−19	118	0.00 × 1000	0.00 × 1000
37	1.26 × 1000	8.66 × 10−02	5.60 × 10−05	4.67 × 10−02	78	1.85 × 1000	0.00 × 1000	119	0.00 × 1000	1.43 × 10−05
38	0.00 × 1000	0.00 × 1000	8.66 × 10−05	3.77 × 10−02	79	4.71 × 10−01	0.00 × 1000	120	0.00 × 1000	0.00 × 1000
39	4.71 × 10−01	7.28 × 10−02	0.00 × 1000	4.32 × 10−02	80	7.45 × 10−01	8.67 × 10−19	121	5.00 × 10−01	0.00 × 1000
40	2.50 × 10−01	7.11 × 10−15	0.00 × 1000	2.77 × 10−02	81	5.00 × 10−01	8.67 × 10−19	122	0.00 × 1000	2.55 × 10−05
41	7.50 × 10−01	0.00 × 1000	0.00 × 1000	2.99 × 10−02	82	0.00 × 1000	8.67 × 10−19	123	0.00 × 1000	2.17 × 10−19
								124	4.71 × 10−01	0.00 × 1000
								125	1.78 × 10−15	0.00 × 1000
								126	5.00 × 10−02	6.52 × 10−06
								127	0.00 × 1000	0.00 × 1000
								128	0.00 × 1000	0.00 × 1000
								129	0.00 × 1000	0.00 × 1000
								130	1.78 × 10−15	0.00 × 1000
								131	0.00 × 1000	0.00 × 1000
								132	5.00 × 10−01	0.00 × 1000
								133	0.00 × 1000	0.00 × 1000
								134		1.08 × 10−19
								135		0.00 × 1000
								136		0.00 × 1000
								137		0.00 × 1000

in appendix. The path from the root to the leaf at the bottom of the tree indicates the influence of a given input x_i on the output y_i under study, with the highest relevance at the top and lowest at the bottom.

The resulting RT for interface deflection y₁ reveals x₅, followed by x₁₀ and x₉ as the most decisive inputs for the favorable flat or slightly convex interface shape (interface deflection y₁ ≤ 0, the nodes marked white in the tree graph in Figure 5). Their importance decreased in the order mentioned above. The most decisive input was the pulling rate x₅, which had a deteriorating effect on the interface deflection, i.e., x₅ should be below 22.8 mm/min to strongly reduce y₁ from average y₁ = 5.84 mm to average y₁ = 1.4 mm. All decisive inputs and ranges of their optimal values that assured the Cz-Ge growth with nearly flat s/l interface, derived from RT analysis, are given in Table 2. From these results, it was also possible to derive the relative feature importance for y₁, as shown in Figure 9a. Again, the most important input concerning its influence on interface deflection had a pulling rate x₅ (32.7%), followed by crystal rotational rate x₁₀ (30.91%), crystal radius x₁ (22.31%), and power of bottom heater x₉ (10.34%). The importance of the input x₁ was high but not conducive to interface flattening. Therefore, x₁ did not appear in the branches marked in white, leading to the optimal value of the interface deflection in Figure 5.

The resulting RT for the crystal growth rate y₂ gave the pulling rate x₅, followed by the position of the side heater x₇ and the emissivity of the radiation shield x4 as the most decisive inputs for the favorable crystal yield (i.e., the maximum value of y₂, corresponding to the nodes in the tree graph in Figure 6 marked white). Their importance decreased in the order mentioned above. All decisive inputs for y₂ and ranges of their optimal values are given in Table 3. The relative inputs’ importance for y₂ is shown in Figure 9b. The pulling rate x₅ with 99.96% influenced the growth rate, followed by 0.025% and 0.005% relative influence of x₃ and x₁₃, respectively. These ML/RT results confirmed the general assumption in Cz crystal growth practice of the equivalence between the growth rate and the pulling rate. In view of this, ML results can serve as the last trusted criterion to explain various assumptions, free from human errors.

The resulting RT for the output y₃ (ratio of growth rate and axial temperature gradient in the crystal) is shown in Figure 7. The most decisive inputs for the favorable value y₃ according to the Voronkov criteria were x₅, x₁₀, x₉, x₁₁, and x₂ (the nodes in Figure 7 are marked in white). As already mentioned, the Voronkov criteria for defect-free Cz-Si crystals can only be applied to Cz-Ge if the crystal is doped with Ga and Sb. The optimal values of the most critical inputs are given in Table 4, and the relative importance of all inputs is shown in Figure 9c. The pulling rate x₅ with 95.96% influenced the y₃, followed by 1.99% and 0.71% relative influence of x₁ and x₉, respectively.

The RT for the combined output y₁₂₃ is given in Figure 8. The favorable values of the input parameters correspond to the minimal leaf mean value of the objective function y₁₂₃ = 0.878. The most decisive inputs for the favorable mean value y₁₂₃ were x₂, x₄, and x₇ (the white nodes in Figure 8). Their related data ranges are given in Table 5, and the relative importance of all inputs is shown in Figure 9d.

In summary, our data-driven study based on the decision trees and synthetic training data revealed the key furnace design and process parameters, their importance, and the ranges of their values for achieving high-quality, high-yield Cz-Ge crystals. Interestingly, the process parameters, especially the pulling rate, had a substantially larger impact on all singular outputs than the furnace hot-zone design parameters. Among the latter, only the crucible size, the axial position of the side heater, and the material properties of the radiation shield were relevant.

In the case of a combined output y₁₂₃ with equal contribution of each singular output y₁–y₃, the obtained results were different, i.e., the furnace design parameters such as crucible size, the axial position of the side heater, and the material properties of the radiation shield were more important than the pulling rate. However, the last result strongly depended on the definition of y₁₂₃ (i.e., on the weighting of contributions of individual outputs and on the justification of consideration of y₃) and could not be generalized.

4. Conclusions

The decision trees are an excellent choice for studying Cz-Ge crystal growth if one needs methods with short training times and acceptable prediction accuracy based on low-volume data (77 in this study) that are able to provide guidelines for understanding the influences of the furnace design and the process parameters on the process economy and crystal quality. The decision trees also provide ranges of 13 input parameters (7 geometric, 1 radiation shield material choice, and 5 crystal growth process parameters) where optimal values of three targeted output variables (s/l interface deflection, growth rate, and v/G ratio) can be found.

Compared to the standard approach to furnace and process development based solely on CFD simulations, e.g., [4,5], characterized by accurate but slow predictions with no ability to generalize and more novel deep learning methods that require a large amount of training data to be accurate and have a black-box nature, e.g., [32], decision trees offer high interpretability with acceptable accuracy at low cost. Against this background, decision trees are a recommendable steppingstone in the development of new processes and equipment for various grown materials and growth techniques before applying high-prediction accuracy methods such as ANN.