3D Numerical Analysis of the Asymmetric Three-Phase Line of Floating Zone for Silicon Crystal Growth

: A numerical simulation has been carried out to study the asymmetric heat transfer, ﬂuid ﬂow, and three-phase line to explain the phenomenon of the spillage of the melt in ﬂoating zone (FZ) silicon growth. A three-dimensional high-frequency electromagnetic (EM) ﬁeld is coupled with the heat transfer in the melt and crystal calculation domains. The current density along the three-phase line is investigated to demonstrate the inhomogeneous heating along the three-phase line. The asymmetric heating is found to a ﬀ ect the ﬂow pattern and temperature distribution of the melt. The three-dimensional solid–liquid interface results show that, below the current supplies, the interface is deﬂected due to strong heating below the current supplies. The calculated asymmetric three-phase line shows a similar trend as the experimentally observed results. The results indicate that the re-melting and spillage phenomenon could occur below the current supplies.


Introduction
With the expansion of the power device market, highly cost effective floating zone (FZ) silicon has been developed. A needle-eye inductor is used to provide high-frequency electromagnetic (EM) heating in FZ silicon with a diameter larger than 100 mm. Because the needle-eye inductor is a one-turn coil, asymmetric heating becomes more obvious when the diameter of silicon is larger according to the previous study [1]. The asymmetric heating can not only affect the melting speed of the feed rod but also affect the shape of the solid-liquid interface between the melt and the single crystal. Inhomogeneous heat generation at the three-phase line could lead to spillage of the melt or to the formation of a bulge at the grown single crystal.
According to a previous study [1], the FZ experiment for silicon crystal was conducted. During the process, the inhomogeneous three-phase line and spillage phenomena have been observed. The three-phase line is deflected downward under the main slit of the inductor. This phenomenon needs to be analyzed by numerical studies. The first 2D numerical study on the FZ was published by Mühlbauer et al. [2]. The 2D global numerical studies were later conducted and were found to provide a good fit with experimental results [3][4][5][6][7]. However, the physical phenomena along the azimuthal direction cannot be analyzed with a 2D numerical model. Therefore, 3D numerical simulation results for the EM field, heat transfer, and impurities have been conducted [8][9][10][11][12]. The asymmetric three-phase line analysis by numerical simulation has not yet been reported. The asymmetric three-phase line, therefore, requires investigation. In a previous study, we proposed a 3D global model to calculate the asymmetric interface [13]. The proposed model could calculate the shape of crystallization interface in 3D by using enthalpy method to track the crystallization interface. The asymmetric 3D deflection of the interface is useful to evaluate the operating conditions in the experiment. Because serious deflection of crystallization interface periphery causes spillage of melt. In the present paper, on the basis of our previous model, the EM field model is improved using the impedance boundary condition offered by the commercial software COMSOL [1]. The three-phase line was calculated to investigate the experimental phenomena.

Numerical Models
Induction heating in the FZ is typically conducted at a high frequency. The typical frequency used in FZ is 3 MHz. The skin depth in the silicon melt is 0.26 mm. The dimension of the skin layer is far too small compared with the calculation domain. Therefore, an impedance boundary condition should be imposed to model such a high-frequency EM field. Figure 1 shows the model constructed with COMSOL.
Crystals 2020, 10, x FOR PEER REVIEW 2 of 10 shape of crystallization interface in 3D by using enthalpy method to track the crystallization interface. The asymmetric 3D deflection of the interface is useful to evaluate the operating conditions in the experiment. Because serious deflection of crystallization interface periphery causes spillage of melt.
In the present paper, on the basis of our previous model, the EM field model is improved using the impedance boundary condition offered by the commercial software COMSOL [1]. The three-phase line was calculated to investigate the experimental phenomena.

Numerical Models
Induction heating in the FZ is typically conducted at a high frequency. The typical frequency used in FZ is 3 MHz. The skin depth in the silicon melt is 0.26 mm. The dimension of the skin layer is far too small compared with the calculation domain. Therefore, an impedance boundary condition should be imposed to model such a high-frequency EM field. Figure 1 shows the model constructed with COMSOL. The alternating current frequency is high, so the EM field is calculated on the basis of the timeharmonic assumption: where is the magnetic vector potential, is electrical conductivity, 0 is the permeability, is the electric potential, is imaginary unit, is the angular frequency, and is generated current density. Table 1. shows physical properties and calculation conditions. The hydrodynamic and temperature fields are calculated using the 3D model ( Figure 2) with the finite volume method (FVM). The model was constructed using the open-source mesh software The alternating current frequency is high, so the EM field is calculated on the basis of the time-harmonic assumption: where A is the magnetic vector potential, σ is electrical conductivity, µ 0 is the permeability, V is the electric potential, j is imaginary unit, ω is the angular frequency, and J e is generated current density. Table 1 shows physical properties and calculation conditions. The hydrodynamic and temperature fields are calculated using the 3D model ( Figure 2) with the finite volume method (FVM). The model was constructed using the open-source mesh software SALOME, and the calculation was performed with FzGlobalFOAM [13], which was developed on the basis of the open-source library OpenFOAM.
Crystals 2020, 10, x FOR PEER REVIEW 3 of 10 SALOME, and the calculation was performed with FzGlobalFOAM [13], which was developed on the basis of the open-source library OpenFOAM. In case of steady-state calculation, the melt flow is solved in the following equations: where is the velocity vector in the melt, is the melt density, is viscosity, is the gravitational acceleration constant, and is the pressure.. The temperature field is calculated by the following equation: where ℎ is enthalpy and is thermal diffusivity. The thermal diffusivities differ at different temperatures in the melt and in the crystal as a result of the thermal conductivity of silicon being temperature dependent. The latent heat is calculated as follows: where is the local crystal growth rate [14] and is the latent heat of silicon. The surface current density distribution from the EM model is coupled with the heat generation at the free surface of the melt in the heat transfer model in the following equation [15]: where and are the surface power density and the current density at the melt free surface, respectively, and skin depth is calculated by the following equation: where f is frequency. The rotation of crystal is considered at the crystallization interface between crystal and melt. The crystallization interface is calculated according to the temperature. The shapes of melting interface and free surface are determined using the data from previous study [13]. The phase boundary of 3D free surface is imposed as a fixed wall. For the temperature field, first-type boundary condition is applied to the melting interface. The temperature at melting interface is melting point. Marangoni and EM forces are imposed as second-type boundary condition at free surface [15]: In case of steady-state calculation, the melt flow is solved in the following equations: where U is the velocity vector in the melt, ρ is the melt density, µ is viscosity, g is the gravitational acceleration constant, and p is the pressure. The temperature field is calculated by the following equation: where h is enthalpy and a is thermal diffusivity. The thermal diffusivities differ at different temperatures in the melt and in the crystal as a result of the thermal conductivity of silicon being temperature dependent. The latent heat Q L is calculated as follows: where V g is the local crystal growth rate [14] and L is the latent heat of silicon. The surface current density distribution from the EM model is coupled with the heat generation at the free surface of the melt in the heat transfer model in the following equation [15]: where Q EM and J s are the surface power density and the current density at the melt free surface, respectively, and skin depth δ is calculated by the following equation: where f is frequency. The rotation of crystal is considered at the crystallization interface between crystal and melt. The crystallization interface is calculated according to the temperature. The shapes of melting interface and free surface are determined using the data from previous study [13]. The phase boundary of 3D free surface is imposed as a fixed wall. For the temperature field, first-type boundary condition is applied to the melting interface. The temperature at melting interface is melting point. Marangoni and EM forces are imposed as second-type boundary condition at free surface [15]: where γ is surface tension. Our previous study gives more explanations of the model in detail [13]. Table 2 shows the process parameters and physical properties. In the calculations, we assume that the shape of free surface is not changing with the temperature and flow fields. The geometry of high-frequency EM model is not updated with fluid flow and heat transfer model. The high-frequency EM model only provides surface current result. The surface current result is used to calculate the EM force F EM and EM heat generation Q EM . For the calculation of fluid flow, EM force F EM is imposed as a fixed gradient boundary condition at free surface of melt. For the calculation of heat transfer, Q EM is imposed as a fixed gradient boundary condition at free surface of melt. EM force F EM and EM heat generation Q EM are updated iteratively if the mean temperature of triple points is not equal to the melting point.

Current Density Analysis
The high-frequency current density is obtained from the EM calculation ( Figure 3). The current density is low below the side slits. The current density is high below the tips of side slits. Compared with the current density below the three side slits, that in the vicinity of the main slit is lower. However, near the three-phase line, the current density below the current supplies is higher. Figure 4 shows the surface current density distribution along the three-phase line. The surface current density distribution below the main slit is steep and high. Except for the main slit, the surface current density below the three side slits is also higher than in the other areas because of the strong magnetic field below the tips of the side slits. This asymmetric EM field at the free surface is used to calculate the EM heat generation in the heat transfer model.

Melt Flow Analysis
3D heat transfer and fluid flow calculation were carried out, and the temperature distribution is shown in Figure 5. The pattern of temperature distribution is similar with that of current density distribution ( Figure 3). The pattern is asymmetric due to the rotation of crystal. The results in the melt are shown in Figure 6. In the cross section paralleled with current supplies (X plane in Figure 6a), the maximum temperature occurs below the side slit tip. The temperature is lower below the main slit than below the side slit. The position of the solid-liquid interface is higher in the left side (below the main slit). In the plane perpendicular to the main slit (Y plane in Figure 6b), the temperature distribution is more symmetric than that of the X plane. The temperature is not strictly symmetric because the crystal rotation shifts the temperature distribution. The flow field is shown in Figure 7. The maximum velocity occurs at the inner triple point (ITP) between silicon melt, polycrystalline feed rod and gas because the temperature gradient at the inner triple point is very large. The Marangoni force at the ITP is strong. According to the asymmetric temperature distribution in Figure 6a, the velocity is asymmetric in X plane (Figure 7a). The velocity at ITP is higher in the right part than that in the left part due to larger temperature gradient below side slit than that below main slit. Ratnieks et al. also mentioned this phenomenon [8]. Figure 7b shows the flow field in the Y plane. The presence

Melt Flow Analysis
3D heat transfer and fluid flow calculation were carried out, and the temperature distribution is shown in Figure 5. The pattern of temperature distribution is similar with that of current density distribution (Figure 3). The pattern is asymmetric due to the rotation of crystal. The results in the melt are shown in Figure 6. In the cross section paralleled with current supplies (X plane in Figure 6a), the maximum temperature occurs below the side slit tip. The temperature is lower below the main slit than below the side slit. The position of the solid-liquid interface is higher in the left side (below the main slit). In the plane perpendicular to the main slit (Y plane in Figure 6b), the temperature distribution is more symmetric than that of the X plane. The temperature is not strictly symmetric because the crystal rotation shifts the temperature distribution. The flow field is shown in Figure 7. The maximum velocity occurs at the inner triple point (ITP) between silicon melt, polycrystalline feed rod and gas because the temperature gradient at the inner triple point is very large. The Marangoni force at the ITP is strong. According to the asymmetric temperature distribution in Figure 6a, the velocity is asymmetric in X plane (Figure 7a). The velocity at ITP is higher in the right part than that in the left part due to larger temperature gradient below side slit than that below main slit. Ratnieks et al. also mentioned this phenomenon [8]. Figure 7b shows the flow field in the Y plane. The presence

Melt Flow Analysis
3D heat transfer and fluid flow calculation were carried out, and the temperature distribution is shown in Figure 5. The pattern of temperature distribution is similar with that of current density distribution (Figure 3). The pattern is asymmetric due to the rotation of crystal. The results in the melt are shown in Figure 6. In the cross section paralleled with current supplies (X plane in Figure 6a), the maximum temperature occurs below the side slit tip. The temperature is lower below the main slit than below the side slit. The position of the solid-liquid interface is higher in the left side (below the main slit). In the plane perpendicular to the main slit (Y plane in Figure 6b), the temperature distribution is more symmetric than that of the X plane. The temperature is not strictly symmetric because the crystal rotation shifts the temperature distribution. The flow field is shown in Figure 7. The maximum velocity occurs at the inner triple point (ITP) between silicon melt, polycrystalline feed rod and gas because the temperature gradient at the inner triple point is very large. The Marangoni force at the ITP is strong. According to the asymmetric temperature distribution in Figure 6a, the velocity is asymmetric in X plane (Figure 7a). The velocity at ITP is higher in the right part than that in the left part due to larger temperature gradient below side slit than that below main slit. Ratnieks et al. also mentioned this phenomenon [8]. Figure 7b shows the flow field in the Y plane. The presence of current supplies does not strongly affect the symmetric flow pattern in the Y plane. The melt converges and flows downward due to strong EM force below the three side slits. However, below the main slit, there is no similar flow pattern driven by EM force because the EM force is weak compared to the side slits. In the actual growth process, with the rotation of the crystal, the flow pattern repeatedly varies from symmetric to asymmetric. The flow separation point, which is the origin of accumulated dopant, also changes repeatedly. These calculation results can explain the asymmetric distribution of the measured radial resistivity distribution in the as-grown crystal [4].
Crystals 2020, 10, x FOR PEER REVIEW 6 of 10 of current supplies does not strongly affect the symmetric flow pattern in the Y plane. The melt converges and flows downward due to strong EM force below the three side slits. However, below the main slit, there is no similar flow pattern driven by EM force because the EM force is weak compared to the side slits. In the actual growth process, with the rotation of the crystal, the flow pattern repeatedly varies from symmetric to asymmetric. The flow separation point, which is the origin of accumulated dopant, also changes repeatedly. These calculation results can explain the asymmetric distribution of the measured radial resistivity distribution in the as-grown crystal [4].   Crystals 2020, 10, x FOR PEER REVIEW 6 of 10 of current supplies does not strongly affect the symmetric flow pattern in the Y plane. The melt converges and flows downward due to strong EM force below the three side slits. However, below the main slit, there is no similar flow pattern driven by EM force because the EM force is weak compared to the side slits. In the actual growth process, with the rotation of the crystal, the flow pattern repeatedly varies from symmetric to asymmetric. The flow separation point, which is the origin of accumulated dopant, also changes repeatedly. These calculation results can explain the asymmetric distribution of the measured radial resistivity distribution in the as-grown crystal [4].

Three-Phase Line Analysis
The solid-liquid interface is affected by the asymmetric heating. Figure 8 shows the 3D shape of the solid-liquid interface. The color legend indicates the vertical height from the bottom of the interface. The symmetric interface center demonstrates that the asymmetric heating at the melt free surface does not strongly influence the center of the interface. However, the interface periphery is obviously inhomogeneous because the inhomogeneous heat is generated at the three-phase line ( Figure 4). The three-phase line is calculated and compared with experimentally observed results [1] in Figure 9. Below the current supplies, the position is low because of the high current density and high heat power below the current supplies. The experiment was conducted in IKZ (Leibniz Institute for Crystal Growth). From the experimental observation (Figure 9a), the three-phase line is deflected below the current supplies. The deflection is regarded as the reason for spillage down the silicon melt. In the calculation results (Figure 9b), along the rotation direction, the three-phase line descends and ascends below the current supplies. In the experimental observation, the deflection is smaller than description in Figure 9a and still visible to the naked eye [1].

Three-Phase Line Analysis
The solid-liquid interface is affected by the asymmetric heating. Figure 8 shows the 3D shape of the solid-liquid interface. The color legend indicates the vertical height from the bottom of the interface. The symmetric interface center demonstrates that the asymmetric heating at the melt free surface does not strongly influence the center of the interface. However, the interface periphery is obviously inhomogeneous because the inhomogeneous heat is generated at the three-phase line (Figure 4).

Three-Phase Line Analysis
The solid-liquid interface is affected by the asymmetric heating. Figure 8 shows the 3D shape of the solid-liquid interface. The color legend indicates the vertical height from the bottom of the interface. The symmetric interface center demonstrates that the asymmetric heating at the melt free surface does not strongly influence the center of the interface. However, the interface periphery is obviously inhomogeneous because the inhomogeneous heat is generated at the three-phase line (Figure 4). The three-phase line is calculated and compared with experimentally observed results [1] in Figure 9. Below the current supplies, the position is low because of the high current density and high heat power below the current supplies. The experiment was conducted in IKZ (Leibniz Institute for Crystal Growth). From the experimental observation (Figure 9a), the three-phase line is deflected below the current supplies. The deflection is regarded as the reason for spillage down the silicon melt. In the calculation results (Figure 9b), along the rotation direction, the three-phase line descends and ascends below the current supplies. In the experimental observation, the deflection is smaller than description in Figure 9a and still visible to the naked eye [1]. The three-phase line is calculated and compared with experimentally observed results [1] in Figure 9. Below the current supplies, the position is low because of the high current density and high heat power below the current supplies. The experiment was conducted in IKZ (Leibniz Institute for Crystal Growth). From the experimental observation (Figure 9a), the three-phase line is deflected below the current supplies. The deflection is regarded as the reason for spillage down the silicon melt. In the calculation results (Figure 9b), along the rotation direction, the three-phase line descends and ascends below the current supplies. In the experimental observation, the deflection is smaller than description in Figure 9a and still visible to the naked eye [1]. Through the qualitatively comparison, the calculated deflection shows the similar trend with the experimental observation. Because of the inhomogeneity of the interface along the azimuthal positions, the local crystal growth rate is inhomogeneous at the three-phase line. Along the crystal rotation direction, the local crystal growth rate is lower than the pulling rate when the three-phase line descends and the local crystal growth rate is higher than the pulling rate when the three-phase line ascends. The latent heat, which is dependent on the local crystal growth rate, is considered using the method suggested in a previous study [14]. The inhomogeneous local growth rate causes the asymmetric three-phase line deflection below the current supplies. Moreover, when the local growth rate is lower than the pulling rate, the re-melting phenomenon occurs. Figure 10 shows the three-phase line position along the azimuthal direction. The three-phase line is asymmetric because the rotation shifts the fluid and temperature field. The three-phase line is deflected downward in the vicinity of the main slit and side slits as a result of high temperature in the vicinity of the slits ( Figure 5). The position of three-phase line is higher between −90°and 90° because the temperature is lower in the vicinity of the main slit. Below the main slit, the defection is about 1 mm. It is demonstrated that even though the crystal is rotated at 5 RPM, the deflection of three-phase line still exists.

Conclusions
Three-dimensional numerical simulation of the EM field and heat transfer were carried out to analyze the influence of asymmetric heating on the asymmetric temperature distribution and three- Through the qualitatively comparison, the calculated deflection shows the similar trend with the experimental observation. Because of the inhomogeneity of the interface along the azimuthal positions, the local crystal growth rate is inhomogeneous at the three-phase line. Along the crystal rotation direction, the local crystal growth rate is lower than the pulling rate when the three-phase line descends and the local crystal growth rate is higher than the pulling rate when the three-phase line ascends. The latent heat, which is dependent on the local crystal growth rate, is considered using the method suggested in a previous study [14]. The inhomogeneous local growth rate causes the asymmetric three-phase line deflection below the current supplies. Moreover, when the local growth rate is lower than the pulling rate, the re-melting phenomenon occurs. Figure 10 shows the three-phase line position along the azimuthal direction. The three-phase line is asymmetric because the rotation shifts the fluid and temperature field. The three-phase line is deflected downward in the vicinity of the main slit and side slits as a result of high temperature in the vicinity of the slits ( Figure 5). The position of three-phase line is higher between −90 • and 90 • because the temperature is lower in the vicinity of the main slit. Below the main slit, the defection is about 1 mm. It is demonstrated that even though the crystal is rotated at 5 RPM, the deflection of three-phase line still exists. Through the qualitatively comparison, the calculated deflection shows the similar trend with the experimental observation. Because of the inhomogeneity of the interface along the azimuthal positions, the local crystal growth rate is inhomogeneous at the three-phase line. Along the crystal rotation direction, the local crystal growth rate is lower than the pulling rate when the three-phase line descends and the local crystal growth rate is higher than the pulling rate when the three-phase line ascends. The latent heat, which is dependent on the local crystal growth rate, is considered using the method suggested in a previous study [14]. The inhomogeneous local growth rate causes the asymmetric three-phase line deflection below the current supplies. Moreover, when the local growth rate is lower than the pulling rate, the re-melting phenomenon occurs. Figure 10 shows the three-phase line position along the azimuthal direction. The three-phase line is asymmetric because the rotation shifts the fluid and temperature field. The three-phase line is deflected downward in the vicinity of the main slit and side slits as a result of high temperature in the vicinity of the slits ( Figure 5). The position of three-phase line is higher between −90°and 90° because the temperature is lower in the vicinity of the main slit. Below the main slit, the defection is about 1 mm. It is demonstrated that even though the crystal is rotated at 5 RPM, the deflection of three-phase line still exists.

Conclusions
Three-dimensional numerical simulation of the EM field and heat transfer were carried out to analyze the influence of asymmetric heating on the asymmetric temperature distribution and three-

Conclusions
Three-dimensional numerical simulation of the EM field and heat transfer were carried out to analyze the influence of asymmetric heating on the asymmetric temperature distribution and three-phase line deflection. The induction heating at the free surface of the melt is asymmetric because of the asymmetric inductor design. Using the three-dimensional induction heating calculation result, the asymmetric results of temperature distribution and flow field were calculated. The calculation results of the three-phase line show a similar trend as the experimental results. We confirmed that the spillage phenomenon can be caused by asymmetric heat generation at the three-phase line. According to the calculation result, the position of maximum heat generation is close to the three-phase line. In order to reduce the risk of spillage of melt, the position of the maximum heat generation should be designed farther from the three-phase line. For example, decreasing the length of the side slits in high-frequency inductor is beneficial to prevent spillage phenomenon.