Numerical Simulation of CdTe Crystal Growth Using the Vertical Gradient Freeze Technique Assisted by Axial Low-Frequency Oscillations of the Melt

Abstract

The problem of intensification of the melt crystal growth process has been analyzed using CdTe as an actual material. Numerical simulation of 100 mm diameter CdTe crystal growth using the VGF technique has been carried out. The heat–mass transfer was controlled by introducing low-frequency oscillating baffle into the melt, which is a so-called axial vibrational control (AVC) technique. The baffle configuration has been optimized to destroy solid “tails”, which were formed near the crucible walls at high cooling rates due to the low thermoconductivity and the corresponding latent heat. Analysis of CdTe homogeneity range showed that during fast crystal cooling, Te micro precipitations were formed, resulting from the decay of oversaturated Cd-rich nonstoichiometric solid solution during the Bridgman crystal growth technique. After full crystallization, a VGF-grown CdTe crystal stays inside the phase field of the high-temperature wurtzite polymorph. This makes it possible to go through the polymorph transition without Te micro-precipitating using the advantages of the VGF-specific feature of very slow cooling.

1. Introduction

The problem of increasing the efficiency of industrial crystal growing processes has been of concern to both scientists and businesses for the last 30 years due to the growing volume of single crystal production. Moreover, over 90% of industrial crystals are grown from the liquid phase. In turn, more than 50% of these crystals are obtained via directional crystallization of melts (Figure 1a). This applies to 90% of semiconductor crystals, which are mainly obtained using the Czochralski, Bridgman, VGF, and TNM methods and their variations (Figure 1b).

The impressive achievements of the VGF method for production of perfectly structured GaAs crystals with a diameter of up to 200 mm are associated with the implementation of a multi-crucible configuration of the growth installation [1]. In this installation, despite the low growth rate of the VGF process, it was possible to implement a highly profitable industrial process [2]. However, these achievements are inextricably linked with an increase in the chemical purity of the starting material to a purity of 7–8N [3,4,5].

Figure 1. Distribution by mass of grown crystals of the main methods for the most popular crystals (a) and various technologies for growing semiconductor crystals (b). The results are based on an analysis of the report: “SEMICONDUCTOR CRYSTAL GROWTH SYSTEM MARKET—Global Outlook and Forecast 2023–2030” [6].

Figure 1. Distribution by mass of grown crystals of the main methods for the most popular crystals (a) and various technologies for growing semiconductor crystals (b). The results are based on an analysis of the report: “SEMICONDUCTOR CRYSTAL GROWTH SYSTEM MARKET—Global Outlook and Forecast 2023–2030” [6].

Semiconductor crystals of cadmium telluride (CdTe) and its solid solutions (Cd_xZn_1−xTe) are the best materials for room-temperature X-ray and gamma-ray detection systems [7], which are widely used in medicine, security systems, and scientific technology. Therefore, interest in these materials has a steady growth trend [8].

To date, the industrial production of high-quality CdZnTe crystals with a diameter of up to 100 mm has been carried out using the THM method [9,10]. This method made it possible to overcome the problem of high-temperature polymorph transition and reduce the concentration of non-stoichiometric defects to the required level [11]. Lowered growth temperatures, low segregation of Zn, and the ability to use the seed crystal are advantages of THM. The achieved mobility–lifetime product, leakage current, and resistivity of the growth crystals are compared with the best available CdZnTe crystals [12]. However, from a commercial point of view, the THM method remains low-profit due to the low growth rate and the complexity of preparing a polycrystalline seed of variable composition.

The growth of structurally perfect CdTe and CdZnTe crystals with a diameter of 100 mm due to the low thermal conductivity of the material provided by the Bridgman method is accompanied by a number of problems, one of which is a concave front at a relatively high crystallization rate [13,14]. It is formed, on the one hand, due to the extremely low thermal conductivity of cadmium telluride, which entails a very slow radial heat transfer. On the other hand, during the CdTe crystallization, a large amount of latent heat is released, and taking into account natural convection in the melt, this leads to its uneven distribution along the crystallization front, which also leads to its sagging.

In Ref. [15], it was shown that the use of the AVC (axial vibrational control) technique made it possible to solve the above problems by creating controlled forced laminar flows in the melt using a submerged inert oscillating disk. However, the authors also declared that in a real experiment, the use of a cylindrical disk eventually led to the formation of a gas bubble along the surface, especially in viscous liquids, which, over time, caused the turbulence of the flow due to the chaotic movement of the gas bubble. It is assumed that this is due to the compressibility effect during crystallization. In the case of free access of the melt to the surface of the forming crystal, hydrodynamic equilibrium is established, and bubble formation does not occur. However, if access of the melt to the surface is limited, which occurs with the introduction of oscillations through a cylindrical disk due to the small distance between the periphery of the disk and the crucible wall at a high frequency of oscillations, a gradual formation of cadmium vapor occurs [16], forming a stable bubble on the surface of the oscillating disk on account of surface tension forces.

This problem can be solved by configuring an oscillating body that is capable of both effectively leveling the crystallization front and providing free access of the melt to the surface of the crystal being grown. In this work, a cylindrical ring is presented as an oscillating baffle for calculating the growth of CdTe crystal using the VGF method, which satisfies all the above requirements. In the present research, we analyzed the configuration of a low-frequency oscillating ring-baffle and its influence on heat–mass transfer during VGF crystal growth of CdTe crystals.

Since the setup of VGF is very complex (especially in multi-crucible configuration), before proceeding with actual growth experiments, we must understand whether it is in principle possible to obtain a flat interface when growing CdTe VGF or not with appropriate cooling rates characteristic of the VGF technique. This was the main purpose of this study.

2. Mathematical Model

Modeling of the crystallization front in the melt was carried out in the approximation of an incompressible linearly viscous fluid in a gravity field, which is described by the Navier–Stokes equation supplemented by the continuity equation:

$${\rho }_0\frac{\partial \overline{v}}{\partial t}+\left(\overline{v}·\nabla \right)\overline{v}=\frac{\mu }{{\rho }_0}{\nabla }^2\overline{v}-\frac{1}{{\rho }_0}\nabla p+{\overline{f}}_g+{\overline{f}}_{cr};$$

(1)

$$\nabla ·\overline{v}=0.$$

(2)

The fluid density and the gravity force in the Navier–Stokes equation was determined by the Boussinesq approximation:

$$\rho \left(T\right)={\rho }_0-{\rho }_0\beta \left(T-T_m\right);$$

(3)

$${\overline{f}}_g=-\left(\rho -{\rho }_0\right)\overline{g}=-{\rho }_0\beta \left(T-T_m\right)\overline{g}.$$

(4)

The law of energy conservation for an incompressible fluid, taking into account crystallization, is written in the enthalpy–porosity formulation [17]:

$$\frac{\partial \left(\rho H\right)}{\partial t}+\nabla \left(\rho \overline{v}H\right)=-\nabla ·\left(k\nabla T\right).$$

(5)

Neglecting kinetic energy, the enthalpy of the melt is determined as

$$H=\epsilon ∆H_m+\underset{T_{ref}}{\overset{T}{\int }}{C}_{p}dT.$$

(6)

Enthalpy–porosity formulation describes the crystallization process as a change in the “porosity” of the medium ε, or a change in the fraction of liquid, which for two-phase systems with unlimited solubility is defined as

$$\epsilon =\left\{\begin{array}{cc}1,\hfill & \hfill T\ge T_l;\\ \frac{T-T_s}{T_l-T_s},\hfill & \hfill T_s<T<T_l;\\ 0,\hfill & \hfill T\le T_s,\end{array}\right.$$

(7)

where T represents temperature, $T_s$ represents the solidus temperature, and T_l represents the liquidus temperature. When crystallizing a phase without decomposition, T_s = T_l, and the authors [18] propose using a different formulation for ε; however, to increase the stability of the numerical model, we adopted

$$T_s=T_m-∆T, T_l=T_m+∆T, ∆T=0.01 K,$$

(8)

here, T_m is the CdTe melting temperature.

The drop in the fluid velocity at ε < 1 (transition to the crystalline phase) is due to the term ${\overline{f}}_{cr}$. The use of the modified Kozeny–Carman equation is widespread:

$${\overline{f}}_{cr}=-C\overline{v}=-B\frac{{\left(1-\epsilon \right)}^2}{{\epsilon }^3+0.001}\overline{v},$$

(9)

where coefficient B is selected manually, and the larger it is, the sharper the change in speed will be, but the less stable the numerical model becomes. In this work, for all calculations, we used $B={10}^6$.

3. Numerical Model Parameters

A numerical simulation was carried out in a two-dimensional axisymmetric approximation for CdTe melt in a crucible with an inner diameter of 100 mm with and without an oscillating body (Figure 2). A cylindrical ring with sharp edges was attached to the central rod using cylindrical spokes. Previously, it was demonstrated that the cylindrical body [18] as well as a plane disk with rounded edges [19] did not cause flows when oscillating in isothermal fluid. So, the above-described ring-baffle constructed from dense graphite was used as an oscillating body in our numerical simulation.

Grid models with an average cell size of 0.5 mm were constructed. For the near-wall layers, additional grinding was applied. For the model with an oscillating ring-baffle, the mesh was also refined with a minimum cell size of 5 μm with a smooth change in size. The example of velocity magnitude gradient distribution close to the ring-baffle sharp edge is presented in Figure 2 (see the right insertion). The calculation was carried out using the ANSYS 14.5 package. If necessary, the adaptive mesh refinement method was used to further refine the mesh during the calculation process. The specific features of the mesh model at AVC application are described in detail in Refs. [15,20].

The CdTe and dense graphite properties used for the numerical simulation are presented in

Table 1. Properties of CdTe and graphite used for numerical simulation.

Parameter	Value	Temperature Range, K	Approximation Type	Reference
	CdTe
Melting temperature, K	1365			[16]
Enthalpy of crystallization, J∙kg−1	209,200			[21]
Density CdTe melt, kg∙m−3	5663.6	1300–1400	Boussinesq	[21]
Thermal expansion, K−1	6.3952∙10−5
Dynamic viscosity, Pa∙s	$8.3215·{10}^{-3}-4.3729·{10}^{-6}·T$	1365–1400	Piecewise-linear	[22]
Heat capacity (solid), J∙kg−1∙K−1	$198.63+0.0392·T$	298.15–1365	Piecewise-linear
Heat capacity (liquid), J∙kg−1∙K−1	251.1831	1365–1400
Thermal conductivity (solid), W∙m−1∙K−1	$11.304-7.3152·{10}^{-3}·T$	1300–1365	Piecewise-linear	[23]
Thermal conductivity (liquid), W∙m−1∙K−1	$7.8909-4.6318·{10}^{-3}·T$	1365–1400
	Graphite
Density, kg∙m−3	1700
Heat capacity (solid), J∙kg−1∙K−1	2148.3
Thermal conductivity (solid), W∙m−1∙K−1	44.165

.

It should be noted that the public domain contains conflicting data on the heat capacity and thermal conductivity of molten cadmium telluride. In Ref. [23], which is devoted to the growth of cadmium telluride crystals via the Bridgman method, the thermal diffusivity coefficients α near the melting point were experimentally determined. Then, from the definition of thermal diffusivity:

$$k=\alpha \rho C_p,$$

(10)

one can determine the thermal conductivity of cadmium telluride near its melting point. The calculated thermal conductivity coefficients are presented in Table 1.

The heat capacity of the CdTe melt was taken as constant and calculated from the equation for the solid state at the melting temperature of the substance [16]. To correctly calculate the enthalpy of CdTe, the heat capacity equation in the solid state was adjusted, taking into account the change in its density in the solid state [21]:

$$C_{p\left(s\right)}^{\prime }=\frac{{\rho }_s}{{\rho }_l}{C}_{p\left(s\right)}.$$

(11)

The boundary conditions for the temperature along the crucible walls were determined using the equation

$$T\left(x,t\right)=T_0+\frac{dT}{dx}·x-\frac{dT}{dt}·t=T_0+{\tau }_{x}x-{\tau }_{t}t,$$

(12)

where T₀ is the temperature that determines the initial position of the crystallization front, τ_x is the temperature gradient along the height of the crucible, and τ_t is the temperature gradient over time (linear cooling rate).

No-slip boundary conditions apply to all solid walls. The vibration modeling of a cylindrical ring was carried out using a dynamic mesh and was specified by the equation.

$$\left\{\begin{array}{c}{v}_x\left(t\right)=A\omega \cos \left(\omega t\right);\\ v_y=0.\end{array}\right.$$

(13)

Here, A is the amplitude of oscillations, $\omega =2\pi f$ is the angular frequency of oscillations, and $f$ is the frequency of oscillations.

The equations were solved using the finite volume method. The gradients were discretized using the least-squares method. The continuity and Navier–Stokes equations were solved jointly using the Coupled scheme; pressure discretization was carried out using the PRESTO scheme. Discretization of speed and energy was carried out using second-order schemes in the calculation without a ring and using first-order schemes with a ring. Time discretization was carried out according to a first-order scheme.

We paid particular attention to relaxation coefficients. The crystallization model in the enthalpy–porosity formulation is extremely unstable, so it is necessary to introduce small relaxation coefficients for almost all variables and increase the number of internal iterations when calculating the time step. The relaxation and accuracy coefficients for the absolute discrepancy, at which constant values of the variables were achieved, are indicated in

Table 2. Relaxation factors and absolute residual criteria used for numerical simulation.

Relaxation Factors	Value	Absolute Residual Criteria	Value
CFL	9	Continuity	10−9 without the ring
Pressure	0.25		10−7 with the ring
Momentum	0.25	Velocities	10−9 without the ring
Body Forces	0.5		10−6 with the ring
Volume Fraction	0.1–0.5	Energy	10−6
Energy	1

.

4. Results and Discussion

4.1. Effect of Linear Cooling Rate

Calculations were carried out in a crucible at linear cooling rates τ_t = 0.2, 1 and 5 K/h in unsteady mode. The initial approximation was specified using a stationary calculation under boundary conditions T₀ = 1363 K, which corresponds to the position of the crystallization front before reaching the conical expansion of the crucible, i.e., the position of the seed crystal.

The value of τ_x in the VGF method is chosen to be constant and small enough to ensure minimal thermal stress, which affects the structural perfection of the crystal. It is technologically possible to provide a linear gradient τ_x = 1 K/cm, which was adopted for all calculations.

In the resulting sections of the crystallization fields (Figure 3), one can observe how, during the crystal growth process, the crystallization front bends from the crucible axis to the wall.

Two parts of the crystallization front can be distinguished: central and peripheral. Note that they have different curvature, and this is explained by the influence of the thermal effect of crystallization. Indeed, if we imagine that crystallization occurs without a thermal effect, then the curvature of the crystallization front is determined solely by the ratio of thermal conductivities in the melt and crystal (see Ref. [24]). However, taking into account the crystallization heat significantly increases the sagging of the front, and leads to the appearance of a linear “tail” part of the crystal at the periphery of the crucible. Figure 4 shows the velocity field in a logarithmic scale (on the left) and the temperature field (on the right) with the “tails” positioned at a height of 8 cm from the bottom of the crucible. Similar fields are shown in Figure 5 at a distance of 14 cm from the bottom of the crucible.

They clearly show how, at different crystallization rates, an isothermal region is established, which has been formed due to the heat of crystallization. In this region, thermal equilibrium is established, and the melt flow rate in it drops to almost zero. This means that the convective part of Equation (5) also tends to zero, and the temperature distribution, and therefore the shape of the front, is determined only by the time term and the thermal conductivity of cadmium telluride in the liquid and solid phases.

In the near-wall region, the axial temperature gradient predominates, which is determined by the boundary condition (12). Towards the center of the crucible, the axial temperature gradient decreases to zero due to the heat of crystallization. On the one hand, this leads to the formation of “tails” and the establishment of a linear front; on the other hand, it leads to a stationary thermal regime closer to the crucible axis, which determines the concave front profile due to the difference in thermal conductivity coefficients.

Consequently, to straighten the crystallization front, it is necessary to ensure the distribution of the heat of crystallization over the crucible due to the convective term in Equation (5); that is, to introduce a sufficient flow rate into the isothermal region.

The complete results of the crystallization process as movies are presented in the Supplementary Material (Videos S1–S3).

The use of the AVC technique involves the creation of a laminar controlled flow, but at high frequencies and amplitudes of vibrations a transition to a turbulent regime occurs. It can be seen that at τ_t = 0.2 K/h, the bending of the crystallization front is the smallest, the size of the “tails” is small, and the total crystallization time is approximately three days, which is acceptable for the technology of growing high-quality single crystals. Consequently, vibrational parameters will be selected for this mode, at which it will be possible to align the crystallization front.

4.2. Unsteady Calculation of Crystallization Assisted by Axial Low-Frequency Vibrations

The oscillation amplitude of a ring-baffle (Figure 2) was 0.1 mm and the oscillation frequency was 50 Hz. Calculation of oscillatory motion requires the choice of a minor step size for accurate calculation of the velocity field. In our calculation, its size was 1/400 of the oscillation period (5 × 10⁻⁵ s), optimized for the AVC technique described earlier [18]. The calculation was carried out until a stationary temperature distribution in the liquid was achieved, which took 60 s.

The oscillatory motion of a baffle in the form of a cylindrical ring in a viscous liquid in an unlimited volume leads to the formation of a flow velocity field that is sinusoidal in time and solenoidal in space. Complicating this model by taking into account the geometric dimensions of the crucible and thermal convection leads to a change in the flow velocity and deformation of the flow shape. In a viscous fluid, the AVC flows generate two regions (Figure 6).

The primary region is adjacent to the ring. The secondary region is flowing around the ring. The formation of the primary region is a consequence of viscous dissipation of velocity from the surface of the ring. Meanwhile, the formation of the secondary region is a consequence of the flow around the body (movement of the fluid mass) during its movement. The primary region is quite small, so the energy distribution in the system is determined primarily by the secondary flow.

The hot AVC flows, as expected, gradually dissolve the solidified “tail”. Note that the “blurring” of the crystallization front is a specific feature of the “mushy zone” model, which is used in FLUENT 14.5 software for solidification processing. Consequently, the exact shape of the front is plotted by the deep-blue region where the liquid fraction is equal to zero (Figure 7).

Melting of the solid “tail” occurs due to the transfer of a colder mass of the melt to a hotter region of the melt when the ring-baffle moves downwards, followed by the movement of a warmer melt to the crystallization front when the ring-baffle moves up (Figure 8).

Besides, at AVC action, a uniform redistribution of the enthalpy of crystallization occurs, followed by the expansion of the isothermal region. The position of the ring closer to the crucible wall also plays a significant role—a more intense sinusoidal axial flow occurs above the “tail” due to the presence of a narrower cylindrical split than in the center, leading to a more intense heat and mass transfer. Finally, the ring itself has a high coefficient of thermal conductivity, which allows it to serve as an additional source of heat, which is shown in the temperature distribution when the disk oscillates, resulting in further intensification of heat transfer. This is reflected in the appearance of an additional bend of the crystallization front, where the section between the “tail” and the central part of the front is parallel to the surface of the disk (Figure 9).

The introduction of axial low-frequency vibrations through the disk with sharp edges leads to the breakdown of clusters in the melt [25], improving the structural perfection of the grown crystal. The integral power of viscous dissipation was determined using the expression

$$P=\underset{V}{\overset{ }{\int }}\mu {\left(\frac{\partial v_x}{\partial y}+\frac{\partial v_y}{\partial x}\right)}^{2}dV$$

(14)

In the case of the disk with sharp edges, the maximum p value reached 8 × 10⁶ W/m³, while the momentum P_w in the whole melt volume reached 0.2 W [26].

The maximum power of viscous dissipation during our ring-baffle oscillation (Figure 10), however, does not allow the required indicators to be achieved. This value was nearly in order less (1.4 × 10⁵ W/m³) than in the case of a plane disk. This is explained by the fact that the movement of the melt when the ring-baffle oscillates is freer than when the plane disk oscillates. There is much less friction between the layers of the moving liquid, because the mass transfer of the liquid melt occurs when the disk oscillates only through the area between the melt and the wall, which leads to an increase in energy viscous dissipation.

4.3. The Role of High Temperature Phase Transition for CdTe Crystal Growth at Congruent Crystallization

To produce CdTe or CdZnTe crystals for matrix X-ray detectors, one should strongly control both nonstoichiometry and dopant concentration, as well as the structural quality of the grown crystals. The structural quality of the grown crystals depends on the lattice rearrangement upon cooling. In general, the greatest changes in the crystal lattice occur when the crystal is cooling within a polymorph transition.

The high temperature wurzite–sphalerite phase transition in CdTe has been a controversial issue for the last 30 years [27,28,29,30,31,32]. In Ref. [33], using the systematic approach to the analysis of homogeneity ranges of cadmium and zinc chalcogenides, we demonstrated that for CdTe, there is a polymorph transition at 1210–1270 K depending on nonstoichiometry (Figure 11). In the case of CdTe direct crystallization from the Cd-saturated melt in the range of polymorphic transition, we observed Te precipitation (see insertion in Figure 11), which resulted from compensation of nonstoichiometric composition passing through the polymorph transition at a fast cooling rate specific to the Bridgman technique. The quantitative analysis showed that the amount of Te precipitation within the experimental accuracy coincided with the ΔX_Te of phase transition from a high-temperature wurzite CdTe to a low-temperature sphalerite CdTe.

When comparing Bridgman and VGF processes (Figure 12, red circles) due to the very low crystallization rate of the VGF technique (Figure 12, green circles), the grown CdTe crystal after the full crystallization from point (1) will be in the field of w-CdTe (see the insertion in Figure 12). With further slow cooling, it is possible to go through the polymorph transition, achieving the equilibrium state by moving sequentially from point (2) to point (3) and then to point (4), followed by the further crystal annealing. Whereas when using the Bridgman technique, the crystallization starts from (1′) point and ends at (2′) point. At a comparatively fast cooling rate and high temperature gradient (20 K/cm [12]) we move from point (2′) to point (3′), and further crystal cooling Te micro precipitations (ΔX_Te) form, resulting from the decay of oversaturated Cd-rich nonstoichiometric solid solution.

On the other hand, we can adjust the raw ingots’ composition in relation to the green points (see Figure 11) through preliminary annealing and go through the polymorph transition without ΔX deviation very close to the stoichiometric composition. Unfortunately, to date, an exact CdTe phase diagram near the polymorph transition is in need of refinement. We hope to achieve this in the near future.

5. Conclusions

An analysis of possible application of the VGF technique to CdTe crystal growth let us conclude that in spite of the very low thermal conductivity of CdTe, it is possible to realize the crystal growth with appropriate rates specific for the conventional VGF applied to GaAs crystal growth. However, in the case of CdTe, we need to use the AVC technique to suppress the formation of solid “tails” near the crucible walls. The “tails” play a negative role by creating tensions in growing crystal, resulting in the generation of a large number of dislocations. Besides, they provoke the formation of twins.

Due to the low crystallization rate, the VGF method assisted by the AVC technique is a promising tool for growth of large-diameter CdTe crystals with adjusted nonstoichiometry. One can easily control Ar pressure up to 10–60 atm to suppress the cadmium vaporization from Cd-rich melt in a VGF growth setup. The latter let us adjust the nonstoichiometry of VGF-grown CdTe crystals.