**1. Introduction**

The floating zone (FZ) method is used to produce silicon (Si) single crystals with high purity. In this method, a feed material is supplied as a polycrystalline rod, which is melted by a high-frequency inductor, and the molten Si flows along the open melting front (see the scheme in Figure 1). This forms a volume of melt, which cools and begins to crystallize at some distance from the inductor, thereby, creating the crystallization interface. The single crystal is then pulled downwards. Currently, needle-eye inductors (with the inner diameter smaller than the feed rod diameter) are used for the FZ process because they allow the use of a larger feed rod and achieve a larger crystal diameter.

To obtain the desired electrical resistivity of silicon crystals grown with the FZ method, it is necessary to precisely control the concentration of dopants—typically boron (B) or phosphorus (P) [1]. The level and distribution of various impurities, e.g., carbon (C), is also important in applications, such as the kilogram definition project [2]. A crucial phenomenon that defines the transport of such species is the segregation at the crystallization interface.

Let us consider a one-dimensional distribution of species concentration *C*(*z*), where *z* is perpendicular to the crystallization interface, which is located at *z* = 0; see Figure 2, left. The equilibrium segregation coefficient *<sup>k</sup>*<sup>0</sup> is theoretically defined as *Cs*(0) *Cm*(0), where *Cs*(0) is the species concentration in solid silicon and *Cm*(0) is in molten silicon at the interface.

However, crystal growth experiments with the same species (same *k*0) can lead to different axial concentration distributions if other process parameters differ [3]. The effective segregation coefficient *<sup>k</sup>*eff <sup>≡</sup> *Cs*(0) *Cm*(*<sup>z</sup>* <sup>0</sup>) is more practical than *<sup>k</sup>*0, because it considers the concentration outside the thin diffusion boundary layer *Cm*(*z* 0), which can be measured or evaluated more easily than *Cm*(0). Therefore, *k*eff also depends on the melt flow. An analytical expression for *k*eff has been derived in the classical work by Burton, Prim and Slichter (BPS) [4]:

$$k\_{\rm eff} = \frac{k\_0}{k\_0 + (1 - k\_0)e^{-\frac{\eta l}{\Theta}}}, \quad d \sim 1.6 D^{\frac{1}{3}} \eta^{\frac{1}{6}} \omega^{-\frac{1}{2}},\tag{1}$$

where *v* is the crystallization velocity, *d* is the thickness of the diffusion layer in the melt, *D* is the diffusion coefficient of the impurity in the melt, *η* is the kinematic viscosity, and *ω* is the crystal rotation rate. The shape of function *k*eff(*d*) for three different species (carbon, phosphorus and boron) is shown in the right part of Figure 2. According to the BPS model, *k*eff ≥ *k*0, and, when the thickness of the diffusion layer increases (for example, due to a decrease in crystal rotation rate), *k*eff increases and saturates at *k*eff = 1 for large *d*. This saturation happens at smaller *d* for the species with a smaller diffusion coefficient.

**Figure 1.** Axially symmetrical scheme of the FZ method, with designations of the most important parts and the used coordinate system.

**Figure 2. Left**—a simplified scheme of species distribution in the crystal (*z* < 0) and the molten zone (*z* > 0) near the interface (*z* = 0). **Right**—the dependence of the effective segregation coefficient *k*eff on diffusion layer thickness *d* (e.g., [4]) for different species.

Figure 2 shows that *k*eff(*d*) can vary in a large interval, particularly for carbon. Moreover, the BPS model uses simplified assumptions of melt flow pattern and does not predict *k*eff correctly in some cases (e.g., when melt flow is radially converging [5] or when microsegregation is taken into account [6]). Therefore, numerical flow simulations are necessary to obtain a more precise estimation of the effective segregation coefficient for FZ Si crystal growth.

Dopant segregation in the FZ process has been already studied in the literature—an axially symmetrical model of quasi-stationary dopant transport in the melt was described as early as 1997 [7]. Due to the increased availability of computational power, transient 3D simulations became possible around 2010 [8], and these kinds of simulations are still relevant now, e.g., for modelling radial resistivity distribution in 200 mm diameter FZ crystals [9]. However, all the aforementioned publications considered only the radial distribution of dopant concentration, without considering the axial distribution or the effective segregation coefficient.

The literature on *k*eff of silicon impurities in crucible-free, inductively heated FZ processes was found to be scarce; therefore, we analyse other Si growth processes. An experiment has been reported that considered boron segregation during the Czochralski (CZ) process [10]. In that study, *k*eff was found to be between 0.75 and 0.98, and it decreased when the melt flow became more intensive. This seemed to agree with the simplified BPS model (1): when the crystal rotation rate increased from 1 to 13 rpm, the diffusion layer became thinner, thus, diminishing *k*eff.

A study where the CZ crystal rotation was 16 rpm [11] reported a slightly smaller *k*eff = 0.73. For crystals grown from a copper crucible using electron beam heating, *k*eff was found to be 0.8 [12]. In yet another CZ experiment [13], boron *k*eff approached 1 as the axial magnetic field increased. However, there are studies of oxygen transport in CZ Si that, on the contrary, indicate a lower *k*eff when magnetic heaters are used [14]. The application of cusp magnetic field also lowered the oxygen content in grown crystals; however, it can be explained not only by the change in *k*eff—oxygen evaporation at the free surface may have played a large role [15].

In this study, we investigated the influence of the melt flow on *k*eff in the Si FZ growth process and obtained the species distribution *C*(*r*, *z*) in the entire grown crystal including the start cone and end cone growth stages, i.e., not only the radial profile *C*(*r*) but also the axial profile *C*(*z*). For the numerical simulations of the phase boundaries, electromagnetic field, melt flow, and species transport in the melt, we used previously developed numerical models.

The *k*eff was then evaluated from the species concentration distribution in the melt. We created a new program that calculates the evolution of the average species concentration in the melt over time, which influences the species distribution along the crystal axis. To cover a wide range of *k*0, we simulated two species: carbon (*k*<sup>0</sup> = 0.07) and boron (*k*<sup>0</sup> = 0.8).

#### **2. Numerical Model**

The present section describes the system of numerical models used for the modelling of the FZ crystal growth process. For the previously published models, only a brief summary and references are given in the following subsections.

This system of models is used to describe the following aspects of the FZ process: the shape of phase boundaries, electromagnetic (EM) field, melt flow, species transport in the melt and, finally, species distribution in the grown crystal. Section 2.1 presents a graphical and textual overview of the used models and data. The section explains how the simulated phase boundaries and EM field are used for the 3D melt flow simulations.

The melt flow simulations, in turn, are necessary for obtaining the effective segregation coefficient *k*eff. When *k*eff is known, we can simulate the evolution of the average species concentration in the melt, which translates to the axial species distribution in the grown crystal. Sections 2.2–2.4 then focus on three main parts of the system of models: the phase boundaries, melt flow and species distribution in the crystal.

#### *2.1. Overview of Modelling Scheme*

The overall scheme of the models and their interactions is shown in Figure 3. The used data (or the obtained results) are shown in grey boxes, while the models are shown in orange boxes. In the beginning, one needs to specify the process parameters (see grey box No. 1): the zone height, inductor frequency, system geometry, etc. For the modelling of phase boundaries (orange box 2), one of two programs can be used: a quasi-stationary program FZone [16], which is suitable only for the cylindrical phase, or a transient program FZoneT [17], which describes the entire process: start cone, cylinder and end cone.

Both programs assume the axially symmetrical (also denoted by 2D) approximation of the shape of silicon parts. FZone and FZoneT programs are used to calculate melt shape (box 3). FZoneT also allows obtaining the melt volume dynamics (box 8): the history of the change in crystal diameter, melt volume, etc. More information about these programs is given in Section 2.2.

**Figure 3.** The overall scheme of the used numerical models (**orange boxes**) and the data that is being exchanged between them (**grey boxes**)[16–19].

Using the obtained melt shape, a 3D finite volume mesh is created for OpenFOAM [18] hydrodynamic simulation (box 6). This includes transient simulation of the melt flow, heat transfer and species transport; see Section 2.3. To obtain boundary conditions for the melt velocity and temperature at the free surface, a 3D inductor shape is used to calculate the induced heat and Lorentz forces [19] (boxes 4 and 5).

OpenFOAM hydrodynamic simulation produces the radial distribution of the species concentration *C*(*r*) in the grown crystal and the effective segregation coefficient *k*eff. This data, together with the melt volume dynamics, is then used in the transient 0D-segregation program (box 9), which predicts the temporal evolution of the average species concentration in the melt [20]. The complete description of the mathematical model for 0D-segregation is presented in Section 2.4. Finally, the results of 0D-segregation program are interpolated on the point grid that is created in crystal domain (box 10); this process is schematically shown in Figure 4.

**Figure 4.** The scheme of species concentration interpolation from OpenFOAM melt flow simulations of the cone stage and cylinder stage (**a**) and 0D-segregation simulations (**b**) on the grown crystal mesh (**c**).

### *2.2. Phase Boundaries*

The description of the program FZone, which has been used for the calculation of phase boundaries, is provided in [16]. The program considers axially symmetrical approximations of phase boundaries. The melting interface, crystallization interface and open melting front are moved according to the heat flux balance; the influence of the melt flow is neglected. The program FZone operates in quasi-stationary approximation: it describes the cylindrical phase of the process, i.e., it calculates the shape of the phase boundaries for a process stage when they do not change in time.

The transient version of this program that simulates growth process dynamics (e.g., a change in the crystal radius during the cone phase) is called FZoneT [17]. The evolution of the crystal diameter and melt volume is saved in a file, which can be later read by the 0D-segregation program. The species concentration inflow from the melting feed rod (the integral from Equation (3)) is also saved in a file. In the present work, the transient FZoneT program was used to simulate the entire crystal growth (the start cone, cylinder and end cone), and the results from three time instants from the start cone stage were selected for melt flow simulations.

#### 2.2.1. Electromagnetic Field

When quasi-stationary FZone is used, the EM field is calculated in 3D using boundary elements and assuming negligible skin depth [19]. The EM field is iteratively recalculated until the shape of the phase boundaries converges. However, in the transient FZoneT version, only an axially symmetrical (also denoted by 2D) EM field is simulated due to limited computational resources. In the 2D EM simulations, inductor slits are modelled by setting an artificial magnetic field source surface density [16], which allows part of the magnetic field lines to penetrate the inner part of the inductor.

To calculate the 2D EM field distribution, axially symmetrical inductor shape was constructed in a way that ensures the best correspondence to the results obtained with a 3D non-symmetrical shape. Quasi-stationary simulations of the corresponding 3D system were performed to calibrate 2D inductor characteristics so that the differences between the induced power on the silicon surface are the smallest (the results are presented in Section 3.2.2). A similar comparison of different EM field models is described in [21].

The initial simulations were performed using real slit parameters in the 2D model. Then, the width and length of the inductor side slits as well as the width of the main slit were changed to achieve better agreement between 2D and 3D systems.

#### *2.3. Species Transport in Melt*

The species transport in melt is simulated with the OpenFOAM hydrodynamic solver described in [22]. This assumes fixed phase boundaries and uses the transient, incompressible, laminar Navier–Stokes equation for melt velocity, with the Boussinesq approximation for the description of thermal convection. Boundary conditions for velocity are as follows:


A standard convection-diffusion equation is used for simulation of the melt temperature with the EM-induced heat and radiation on the free melt surface and fixed temperature on the solid–liquid interfaces. The equation that describes species transport is also of the convection-diffusion type:

$$\frac{\partial \mathcal{C}}{\partial t} + (\vec{v} \,\nabla)\mathcal{C} = D\Delta \mathcal{C},\tag{2}$$

where *C* is the species concentration, *t* is the time, *v* is the melt velocity, and *D* is the species diffusion coefficient. The initial conditions were uniform: *C* = 1 arb.u. (arbitrary unit). The boundary conditions for concentration are as follows:


#### *2.4. 2D Species Distribution in Crystal*

To simulate the axial distribution of species in the crystal, the time-dependence of the average species concentration in the melt *Cm*(*t*) was considered. It was assumed (and later supported by simulation results; see Section 3.3.1) that the *v* and *Cm* in the melt reach a quasi-stationary state sufficiently fast in comparison to the characteristic time of axial changes of concentration in the grown crystal. Therefore, we assumed that the average concentration of species in the crystallizing silicon at any given time *t* is directly proportional to *Cm*(*t*). Thus, the following equation of species mass conservation was proposed:

$$
\Delta(\mathbb{C}\_m V\_m) = \int \mathbb{C}\_F v\_F \,\Delta t \,\mathrm{d}S - k\_{\mathrm{eff}} \mathbb{C}\_m \Delta V\_{\mathrm{out}} \tag{3}
$$

where *Vm* is the melt volume, *CF* is the species concentration in the feed rod at the melting interface, *vF* is the feed rod pushing rate, Δ*t* is the simulation time step, d*S* is the surface element of melting interface and open melting front, *k*eff = *CC*/*Cm* is the effective segregation coefficient, Δ*V*out is the volume of the crystal that crystallizes during the time step, and *CC* is the average species concentration in the layer that crystallizes during the time step. The value of *<sup>k</sup>*eff is obtained from OpenFOAM simulations (see Section 3.3.1) and can depend on the process parameters, e.g., the crystal diameter *DC*.

Equation (3) can be rearranged to express Δ*Cm*, which then is used to calculate *Cm* iteratively with the time step Δ*t*:

$$
\Delta \mathbb{C}\_m \, V\_m = -\mathbb{C}\_m \, \Delta V\_m + \int \mathbb{C}\_F \upsilon\_F \, \Delta t \, \mathrm{d}S - k\_{\mathrm{eff}} \mathbb{C}\_m \Delta V\_{\mathrm{out}}.
$$

$$
\Delta V\_{\rm in} = \int \upsilon\_{\rm F} \,\Delta t \,\mathrm{d}S
$$

$$
\Delta V\_{\rm m} = \Delta V\_{\rm in} - \Delta V\_{\rm out}
$$

$$
\mathsf{C}\_{\rm m}(t + \Delta t) = \mathsf{C}\_{\rm m} + \Delta \mathsf{C}\_{\rm m} = \mathsf{C}\_{\rm m} + \frac{\mathsf{C}\_{\rm m}(\Delta V\_{\rm out} - \Delta V\_{\rm in}) + \int \mathsf{C}\_{\rm F} \upsilon\_{\rm F} \Delta t \,\mathrm{d}S - k\_{\rm eff} \mathsf{C}\_{\rm m} \Delta V\_{\rm out}}{V\_{\rm m}}.\tag{4}
$$

Equation (4) is implemented using Python language. This program is called 0D-segregation, because it describes transient species segregation without spatial dimensions, i.e., disregarding spatial concentration distribution in the melt. The program source code is published in [20]. Two functionalities of the program are available:

	- Due to the assumption of constant pulling velocity, *L* ∝ *t*.
	- Cone surfaces are approximated as having constant slope, and thus *DC*(*t*) ∝ *t*. • The free surface height above the external triple point is assumed to be constant even during the cone phases because it is impossible to predict its evolution for an
	- arbitrary crystal shape (without experimental data); therefore, *Vm*(*t*) ∝ *DC*(*t*)2. • The crystallized volume is proportional to the crystal cross-section: Δ*V*out = *vCSC*Δ*t*, where *vC* is the crystal pulling velocity, and *SC* = *<sup>π</sup>* <sup>4</sup> *DC*(*t*)<sup>2</sup> is the crystal cross-section area. Therefore, Δ*V*out ∝ *DC*(*t*)2.
	- Due to silicon mass conservation, Δ*V*in = Δ*Vm* + Δ*V*out.

After *Cm*(*t*) has been calculated, it is combined with *C*(*r*), which was obtained during species transport simulations in melt; see the scheme in Figure 4. The final result is *C*(*r*, *z*) distribution in the grown crystal.

### **3. Results**

#### *3.1. Description of the Experiment*

The present work uses process data from a 4" FZ crystal growth experiment performed at the IKZ (Leibniz Institute for Crystal Growth, Berlin) . The main parameters are listed in Table 1. A standard one-turn inductor with three side slits is used; see Figure 5 and our earlier work [25,26]. The shape of the grown crystal as well as the height of the molten zone can be seen in Section 3.2.3, where a comparison to transient simulation results is performed. Process photographs with the visible parts of the phase boundaries were used to validate the simulation of the cone phase in particular. Note that characterization of the grown crystals is not discussed in the present study and will be addressed in further publications.



#### *3.2. Phase Boundaries*

#### 3.2.1. Quasi-Stationary Simulations

First, quasi-stationary simulations are performed with FZone program to obtain the shape of phase boundaries for the system described in Section 3.1, which are later used in hydrodynamics simulations in Section 3.3. 3D EM simulations are performed to precisely describe the EM field created by the high-frequency inductor. An example of the results of 3D EM simulations (induced heat on silicon free surface) is shown in Figure 5.

**Figure 5.** 3D high-frequency electromagnetically induced heat sources on silicon surfaces used for 3D melt flow simulations. The one-turn inductor is shown in grey.

#### 3.2.2. Influence of Three-Dimensionality of the EM Field

Three-dimensional melt flow simulations require 3D distributions of induced heat and Lorentz force, which can only be achieved using a 3D model of the high-frequency inductor shape. On the other hand, transient simulations can be performed only with a 2D EM field model, which approximates the 3D features of the asymmetrical EM field (see Section 2.2.1).

Induced current on the open melting front and the free melt surface of silicon are compared between different EM models in Figure 6. For the 3D model, the azimuthally averaged radial distribution is shown. In the case of the 2D model, the effect of inductor slits can be considered only approximately, which results in a significantly different distribution. Various parameter studies were conducted to obtain the optimal inductor parameters (modified slit width and length).

**Figure 6.** Comparison of induced currents on silicon surfaces in the case of 2D and 3D simulations.

The resulting induced electrical current distribution from the optimized 2D model is similar to the 3D model results with the same defined inductor current. The main differences can be observed on the free melt surface and in the outer part of the open melting front. This results in slightly different shapes of the phase boundaries, which are shown in Figure 7. The internal triple point radius is slightly smaller, which results

in a different melting interface and free melt surface shape. However, the 2D model is sufficiently precise to describe the melt volume dynamics, because the differences are large only in the central part of the melt.

**Figure 7.** Comparison between shapes of the phase boundaries for optimized 2D and 3D EM simulations.

#### 3.2.3. Transient Simulations

While the quasi-stationary approach used in FZone was applied in the present work for the cylindrical growth stage, the transient model should be considered for the initial and final stages of the process when the crystal or feed rod is short and the diameter is changing. Therefore, a transient simulation of the temperature field and phase boundaries has been conducted using FZoneT, starting from a small cone with a diameter of 20 mm and a length of 12 mm at time *t* = 0 min. A rather good agreement with the experiment was obtained for the shape of the phase boundaries during the cone and cylindrical growth stages; see Figure 8.

**Figure 8.** Comparison of the shapes of phase boundaries in the transient FZoneT simulation with experimental photos during the cone growth (0–30 min) and cylindrical growth (40 min) stages.

**Figure 9.** Comparison of transient FZoneT simulation results to the experiment during the entire growth process.

The time-dependent inductor current (instead of the power used in the experiment) and feed rod velocity were adjusted to achieve good agreement with the experimentally measured crystal diameter and zone height as demonstrated in Figure 9. A slightly imperfect prediction of the zone height by FZoneT could be caused by differences in input data (inductor current and feed velocity) or by model limitations (e.g., a fluid film model describing the melting of a macroscopically smooth open melting front).

The feed diameter *DF* was not explicitly measured in the experiment; therefore, only the curve from FZoneT is plotted in Figure 9. It was verified that *DF* during the cylindrical growth stage was 90 mm both in the experiment (from process photos) and FZoneT simulation in accordance with the steady-state mass conservation *DF* = *DC* <sup>√</sup>*vC*/*vF*.

An example of the calculated global temperature field and shape of phase boundaries during the entire process is shown in Figure A1 of Appendix A.

### *3.3. Species Transport in Melt*

Using the simulated phase boundaries, that were obtained using FZone or FZoneT and presented in the previous section, 3D melt flow simulations were run with OpenFOAM using the model described in Section 2.3. The most important physical and numerical parameters are given in Table 2. An example of the simulation results for the carbon transport during the cylindrical phase is shown in Figure 10. Despite the non-symmetrical high-frequency inductor, the carbon concentration field is mostly axially symmetrical.


**Table 2.** The physical [27] and numerical parameters that were used for OpenFOAM simulations.

**Figure 10.** The time-averaged meridional melt velocity in the *xz* plane, temperature in the *yz* plane and carbon concentration on the crystallization interface, simulated for the cylindrical phase (*DC* = 102 mm). The direction of the main inductor current suppliers corresponds to the *x* axis.

The OpenFOAM simulations can be performed only using a fixed melt shape (stationary phase boundaries). The FZone calculation of the cylindrical growth stage with a diameter of 102 mm was used. In addition, the results from FZoneT with several cone diameters were selected to simulate the melt flow in a cone phase: 80, 60 and 45 mm.

The simulated melt velocity in the vertical plane parallel to the main inductor slit is shown in Figure 11, left. Comparison between the EM and Marangoni force densities on the free surface shows that the EM forces dominate; see Figure 12. The EM force is creating an inwards-directed flow at the outer part of the melt. The melt flow is more intensive for the smaller *DC* due to larger EM field gradients.

The flow pattern is similar between two stages with larger *DC*, where the largest velocities (exceeding 1 cm/s) appear only near the free surface and in the neck region.

The decrease in the melt flow intensity during the increase in *DC* may be related to various factors. One of them could be the EM force redistribution due to the change in the free surface shape. The influence of crystal rotation on larger melts could also play a role: when the azimuthal melt velocity is higher, it is more difficult to induce a radial flow [7]. Another possible explanation is different surface-to-volume ratios, which decrease the influence of the surface EM force for larger *DC*.

**Figure 11.** The time-averaged melt velocity (**left**) and carbon concentration (**right**) in the vertical cross-section of the melt, simulated for the cylindrical phase (*DC* = 102 mm) and cone phase (*DC* < 102 mm).

**Figure 12.** Absolute values of the EM and Marangoni force density on the free melt surface during the cylindrical phase (*DC* = 102 mm).

The carbon concentration distribution is shown in Figure 11, right. Arbitrary units are used, and normalization is made with respect to the concentration in the feed rod, i.e., the *C* = 1 boundary condition is set on the melting interface. Concentration contours, despite the averaging over 100 s, are not always smooth—possibly due to the chaotic temporal behaviour of the melt flow or due to minor numerical effects, e.g., mesh non-orthogonality.

As of segregation with small *k*0, a thin boundary layer is formed at the crystallization interface; see Figure 13. The width of this layer approximately corresponds to the values calculated using the BPS model, Equation (1): 0.020 mm for carbon and 0.024 mm for boron. The carbon concentration at the crystallization interface is shown in Figure 14. The obtained concentration distributions for large *DC* are axially symmetrical due to crystal rotation and due to partial axial symmetry of the boundary conditions.

For smaller *DC*, only minor deviations from the axial symmetry in the concentration distribution are present, and they do not significantly affect the azimuthally averaged radial distribution. There is no experimental data about such 2D distributions known to the authors.

**Figure 13.** The distribution of carbon (**left**) and boron (**right**) concentration in the melt, near the crystallization interface (*DC* = 102 mm). The normal coordinate is denoted by *n*, and different radial positions *r* are considered. Boundary layer thicknesses from the BPS model (1) are shown with vertical black lines.

From the species concentration values at the interface, as shown in Figure 14, the radial distribution of species atoms in the grown crystal can be calculated; see Figure 15. For small *DC* (60 mm and less), two distinct concentration maxima exist: one in the crystal centre and one near the rim. This effect can be explained by the intensive melt vortex, as shown in the top images in Figure 11. Such a central maximum of species concentration in small crystals (*DC* = 40 mm) was also observed in experiments [30]. When *DC* is large (*DC* = 102 mm), the species concentration maximum is obtained at *r* ≈ 20 mm. It corresponds well to the resistivity maximum that is obtained due to the distribution of phosphorus or other dopant elements in typical 4 silicon crystals [31].

**Figure 14.** The time-averaged carbon concentration on the crystallization interface, simulated for the cone phase (*DC* < 102 mm) and cylindrical phase (*DC* = 102 mm).

**Figure 15.** Normalized radial distributions of carbon (**left**) and boron (**right**) concentrations in the grown crystal (yellow line for the cylindrical phase, violet and orange lines for the cone phase), using crystal rotation rates 6 rpm (solid lines) and 12 rpm (dashed lines).

#### 3.3.1. Effective Segregation Coefficient

From OpenFOAM melt flow simulations, the effective segregation coefficient can be calculated:

$$k\_{\rm eff} = \frac{\mathbb{C}\_{\rm C}}{\mathbb{C}\_{m}} = \frac{k\_0 \, \mathbb{C}\_{m, \rm crys, inter}}{\mathbb{C}\_{m}},\tag{5}$$

where *CC* is the concentration of species in the crystal, *Cm* is the average species concentration in the melt, and *Cm*, crys. inter. is the average species concentration on the crystallization interface, calculated by averaging the concentration values from the faces of OpenFOAM calculation cells.

Values of *k*eff converged relatively rapidly—after less than 200 s, it did not change more than ±2%. This means that, until the saturation of *k*eff(*t*) curve, less than 12 mm of crystal is grown, which is a small part in comparison to total crystal length (400 mm and more). The obtained values of *k*eff are summarized in Figure 16. The effective segregation coefficient monotonously increases as the crystal diameter increases. This can be explained by considering a simplified scheme in Figure 2: when *DC* is larger, the melt motion is less intensive, the boundary layer remains undisturbed, and the drop of species concentration in the melt is sharper. Therefore, *Cm*(*<sup>z</sup>* <sup>0</sup>) becomes closer to *Cs*, and *<sup>k</sup>*eff <sup>≡</sup> *Cs Cm*(*z* 0) increases.

The obtained range of carbon *k*eff variation only partly intersects with the range calculated using the BPS model (1), where *k*eff is in the interval from 0.3 to 0.5. This could possibly be explained by the simplifications during the derivation of the BPS model, which does not capture all features of the 3D melt flow in the FZ process. It should also be noted that the used carbon diffusion coefficient *D* is rather uncertain (its values are given in [28] with ±30% precision). The range of boron *k*eff, in turn, agrees with Equation (1) rather well.

**Figure 16.** The effective segregation coefficient *k*eff of carbon (**left**) and boron (**right**) obtained for different crystal diameters *DC* and crystal rotation rates of 6 rpm (black) and 12 rpm (orange).

The obtained *k*eff value may depend on the used mesh due to the combination of extremely low *D* and low *k*<sup>0</sup> for carbon. During the research, mesh influence studies were performed, and it was shown that, in order to obtain consistent results, the cell size at the crystallization interface should not be larger than 0.03 mm in the normal direction. This criterion was ensured in all meshes that were used for the melt flow simulations described above.

3.3.2. Increased Crystal Rotation Rate

The crystal rotation rate *ω<sup>C</sup>* was increased from 6 to 12 rpm in two stages—for *DC* = 102 mm and *DC* = 60 mm—in order to investigate if it is possible to decrease carbon *k*eff and thus improve the crystal purification. The comparison of meridional velocity fields revealed only small differences; see Figure 17. Vertical motion (downwards in the centre of the melt and upwards at the middle of the radius) became slightly stronger when *ω<sup>C</sup>* was increased to 12 rpm.

Therefore, the central maximum in the radial distribution of the carbon concentration became more pronounced; see Figure 15 (dashed lines). As shown in Figure 16, an increase in *ω<sup>C</sup>* from 6 to 12 rpm led to a decrease in *k*eff from 0.36 to 0.31 for *DC* = 102 mm and decreased *k*eff by half for *DC* = 60 mm. This means that the high *ω<sup>C</sup>* could help to achieve a slightly higher average crystal purity; however, an increase in the *C*(*r*) profile maximum is predicted, as shown in Figure 15, which will make the 2D *C*(*r*, *z*) distribution in the grown crystal less homogeneous.

**Figure 17.** The time-averaged melt velocity in the vertical cross-section of the melt, simulated for the cylindrical phase with *DC* = 102 mm and crystal rotation rates of 6 rpm (**left**) and 12 rpm (**right**).
