A Coupled Approach to Compute the Dislocation Density Development during Czochralski Growth and Its Application to the Growth of High-Purity Germanium (HPGe)

Abstract

The evolution of the dislocation density during Czochralski growth is computed by the combination of global thermal calculations and local computation of the stress and dislocation density in the crystal. The global simulation was performed using the open-source software Elmer (version 8.4) and the local simulation with the open-source software MACPLAS (version of 23.1.2023). Interpolation both in space and time was used to transfer the boundary conditions from the global simulations to the local model, which uses a different mesh discretization and a considerably smaller time step. We applied this approach to the Czochralski growth of a high-purity Ge crystal. The heater power change predicted by the global model as well as the final dislocation density distribution in the crystal simulated by the local model are correlated to the experimental results.

1. Introduction

Growth of germanium (Ge) semiconductor single crystals is typically performed using the Czochralski (Cz) method. Among many applications of Ge crystals [1], we would like to highlight radiation detectors to make precise measurements of high-energy rays like gamma or X-rays [2]. A hydrogen atmosphere is needed to ensure the high purity during the growth process, but this has the drawback of forming hydrogen di-vacancies, which act as charge traps in the detector. Dislocations can reduce such defects and diminish the charge trapping effects of them. Dislocation densities in the range of ${10}^3{\mathrm{cm}}^{-2}$–${10}^4{\mathrm{cm}}^{-2}$ lead to a good performance [3]. Consequently, engineering of the dislocations is required to obtain a nearly homogeneous distribution over the crystal in the required range. It is well known that the multiplication of dislocations is influenced by the stresses in a growing crystal. For example, low thermal stresses and dislocation densities can be achieved by specially designing the hot zone to reduce the temperature gradients [4].

Such dislocation engineering would be very helpful in conjunction with the Legend project [5], with the aim of detecting the neutrinoless double-$\beta$ decay. High-purity germanium (HPGe) crystals are enriched in the isotope ${}^{76}$Ge, which is predicted to undergo a neutrinoless double-$\beta$ decay. Thus, the HPGe detector acts both as source and detector. Such experiments have been performed, e.g., in the GERDA project [6]. Further investigations will be performed within the Legend project [5], applying 200 kg of Ge detectors in the first step and 1000 kg in the second.

Numerical simulations are a valuable tool for detailed investigations of the Cz crystal growth, including such aspects as the global heat transfer, melt flow, thermal stresses, and dislocation density distribution. While a large amount of literature exists on modelling the Cz Ge crystal growth, see, e.g., [7,8,9], computing the dislocation density evolution is still a great challenge. The time evolution of the temperature field in the growing crystal has to be calculated accurately considering the growth process in the furnace. On the other hand, the dislocation density dynamics has to be modelled in a growing crystal, which in most cases has been achieved by applying the Alexander–Haasen (AH) model [10]. This describes stress relaxation through plastic deformations and the multiplication of the dislocations depending on the local temperature and the effective stress. The Alexander–Haasen model does not take into account the propagation of dislocations. Thus, it is a simplification of reality, in which different kinds of dislocations can occur, such as screw and edge dislocations. In Czochralski-grown Ge crystals, screw dislocations are dominant [11]. In the classical Alexander–Haasen model, there is also no preferential direction, like gliding planes, defined where the line defects move. A deeper insight into the dislocation dynamics could be gained by discrete dislocation dynamics, with high computational costs for small spatial domains (see, e.g., [12]). Despite all of these restrictions, the Alexander–Haasen model has been used for various materials, such as GaAs and InP [13], multi-crystalline Si [14], and SiC [15].

One of the first realistic simulations of dislocation evolution during Cz growth of semiconductor crystals was carried out by Miyazaki and Okuyama [16]. They employed the AH model and considered Si, InP, and GaAs crystals with a diameter of 8${}^{\prime \prime }$. The temperature field was calculated at several discrete crystal lengths, and interpolated on a different mesh for the dislocation density simulations. The model was later extended to take into account the crystal anisotropy [17], finding that the dislocation density for the $\langle 111\rangle$ growth direction is lower than for the $\langle 001\rangle$ direction.

A similar approach of performing separate steady-state global heat transfer simulations which are then interpolated onto the mesh for dislocation density simulations has also been used by other researchers, see, e.g., [18,19]. However, recalculation of the temperature equation in the local crystal model with a finer mesh and a much lower time step was not performed. This could potentially lead to a reduced spatial resolution of the temperature field affecting the stress distribution, and to neglected transient effects.

A fully transient calculation was performed by Artemyev et al. for the Czochralski growth of Ge [8]. They used the software package CGSim to perform a transient calculation of the crystal growth in the cylindrical part of the process. Thus, the grids in melt and crystal were deformed during the computation. However, for the computation of the dislocation density the grid was different and not changing in time [20]. This requires adding nodes at the propagating crystal/melt interface and, subsequently, setting dislocation densities in the new nodes. They used two approaches for the setting, as described in more detail in [8]. They found a reasonable agreement with the six experimental data points of etch pit density (EPD). Because they considered a Czochralski growth with very low thermal gradients, provided by a very effective thermal insulation, there was no increase in the dislocation density from the middle of the crystal to the end.

For the growth of high-purity germanium sources, where impurities have to be minimised, such insulation is not possible. Consequently, the thermal gradients are much larger; also due to the high thermal conductivity of hydrogen (H${}_2$), which is used as the gas environment. Wang et al. [21] obtained the axial temperature difference by measuring the temperature at two points using thermocouples yielding to thermal gradients of about 25 $\mathrm{K} {\mathrm{cm}}^{-1}$ in the beginning up to 80 $\mathrm{K} {\mathrm{cm}}^{-1}$ at the end of the growth process. They used the axial temperature gradient to estimate the dislocation density in the crystal in the sense of a very simple approximation.

In this paper, we present a coupled approach to computing the evolution of the dislocation density and apply it to the growth of an HPGe crystal. Details of the experimental setup and process are described in [22,23]. The coupled approach consists of an axisymmetric global simulation of the temperature field (with a larger time step) and a local calculation of the temperature, stress, and dislocation density only in the crystal (with a smaller time step). The coupling is considered in one direction, i.e., the temperature at the surface of the crystal is transferred from the output of the global simulation to the boundary conditions in the local simulation. For the global simulation, we use the open-source software Elmer, and for the local one the open-source software MACPLAS [24], which has been already applied to a parametric study of dislocations in Ge crystals [25]. The coupled approach presented here is general and can be used for any Czochralski process.

2. Numerical Model

2.1. Global Calculation Using Elmer

Global axisymmetric simulations of the growth process were performed with the open-source software Elmer (version 8.4) [26], which is a finite element program with various solvers, such as for inductive heating, heat transfer, wall-to-wall radiation, melt convection, and phase change. For the mesh creation we applied gmsh (version 4.11.1) [27] through its Python API using an object-oriented approach based on “objectgmsh” [28], which simplifies meshing of complex geometries. Recently, a Python interface for Elmer, named “pyelmer”, has been developed [29]. Based on this interface, a framework for Czochralski growth called “openCGS” [30,31] has been established. The “openCGS” framework was extended to perform quasi-transient calculations as developed earlier [9]. In the current version, geometries for Czochralski growth are described in a more systematic way (see Figure 1): (1) crystal and crystal holder (type I, movable by crystal pull velocity); (2) crucible and crucible holder (type II, movable); (3) other parts (type III, fixed). In the Czochralski growth experiments for HPGe, the crucible position is fixed, i.e., objects of type II do not move. Thus, we only need to compute the position of the crystal and its holder.

For geometry definition, openCGS has been extended with a set of functions that draw the complex furnace geometry based on the characteristic points of each body. A parameterisation allows the movable plant parts to be moved according to the current crystal length. The contact areas between the bodies are determined automatically and do not require further adjustment, e.g., when the melt height drops from the cylindrical to the bottom part of the crucible. For quasi-transient simulations, a new class was implemented in openCGS, which automatically starts the individual steady-state simulations and performs basic evaluations. This general openCGS code is available open source.

The crystal shape is given as a set of data points $(R,L)$, which might originate from a grown crystal or define a target crystal shape. The first step is to compute the corresponding crystal mass. Because the crucible is completely emptied at the end of the growth process, the mass of the crystal corresponds to the mass of the melt at the beginning of the growth process. In the model, the surface of the crystal is represented by a spline function, which means we have a continuous function which is used to compute the growth angle for every time step. From the growth angle and the contact angle (${14}^{\circ }$ for Ge) together with the crystal radius, the meniscus height $h_{\mathrm{men}}$ is computed using the approach of Hurle [32]. Further details can be found in [9]. The melt volume is computed considering the meniscus but under the assumption of a flat crystal/melt interface. The interface shape is computed in the run and cannot be considered at this step. The (numerical) melt height $h_{\mathrm{melt}}$ is defined as shown in Figure 1.

The time loop is realised by increasing the length of the crystal by $\Delta l$ in every time step. The growth velocity $v_{\mathrm{growth}}$ results from the pulling velocity $v_{\mathrm{pull}}$ and the decreasing melt height. $v_{\mathrm{growth}}$ is larger than the pulling velocity $v_{\mathrm{pull}}$, since the melt level is always decreasing. $v_{\mathrm{growth}}$ enters the computation of latent heat release in the next time step. For the calculations presented in this paper, we used $\Delta l=2\mathrm{mm}$ and $\Delta l=1\mathrm{mm}$ to check the influence of the time step in the quasi-transient calculation.

Having the geometry of melt and crystal for the global calculation, the induced heat is computed. In the thermal field calculation, the induced heat is multiplied by a factor to meet the melting point temperature at the three-phase point. The thermal field calculation comprises wall-to-wall radiation, heat conduction, and convective heat transfer in the melt and gas. The convection in the melt and the gas is obtained by directly solving the Navier–Stokes equation. To achieve a converged solution, the viscosity was artificially increased.

In the melt we consider buoyancy convection, forced convection by the crystal, and crucible rotation, as well as Marangoni convection at the free surface. For the gas convection we consider buoyancy convection and set the velocity at all boundaries to zero, except at the inlet and the outlet. For comparison with MACPLAS we also compute the elastic stress in the crystal.

In every time step, we write out the temperature along the surface of the crystal. These data are used as the boundary condition in the local calculation with MACPLAS. The correctness of this procedure was checked by comparison of the stress field obtained with Elmer and a steady-state computation with MACPLAS.

2.2. Local Calculation Using MACPLAS

Since the results of the global calculation are only available at discrete crystal lengths, a local transient model has been developed to predict the time-dependent temperature and dislocation density distributions in the crystal. The local model uses the previously developed open-source solver package MACPLAS [24] which is based on the open-source finite element library deal.II [33]. MACPLAS has recently been applied to analyse thermally induced dislocation generation in floating-zone Si crystals [34] and dislocation dynamics in a parametric model for Czochralski Ge growth [25]. For the present study, the MACPLAS crystal growth solver has been extended with an option to import external temperature boundary conditions $T(r,z,t)$.

The time-dependent temperature field in the crystal is described by the equation

$$\rho c_p\frac{\partial T}{\partial t}=\nabla ·(\lambda \nabla T),$$

(1)

where $\rho$ is the density, $c_p$—specific heat capacity, and $\lambda$—thermal conductivity. Heat transfer by convection is not included in the temperature equation, but is considered by the temperature update in each time step according to the shifts in the mesh nodes [25]. We chose to solve the local temperature Equation (1) instead of simply interpolating the whole field from the global results, in order to obtain a more accurate T distribution on a finer local mesh and take into account the transient effects (at least partially, since temperature at the crystal surface is still taken from the steady-state global results).

A zero-flux boundary condition is set on the crystal axis; Dirichlet (first-type) boundary conditions are applied at the melt/crystal interface (melting point) and crystal side surface (non-homogeneous temperature profile from the global simulation). A steady-state temperature distribution is calculated and used as the initial condition at the beginning of the simulation at $t=0$.

The crystal geometry and temperature field on the crystal surface (for Dirichlet boundary conditions) from global Elmer simulations are used in the local model. Since Elmer results are only available at certain crystal lengths, linear interpolation in time between the two closest crystal positions is employed to obtain the necessary data for the local simulation, similar to previous studies, cf. [16,18]. To correctly treat the time-varying crystal shape, the radial and axial coordinates are normalised by the crystal radius and crystal length, correspondingly.

For simplicity, the stretching mesh approach is used in the local model, i.e., without changing the number of points and without mesh regeneration. Additional aspects related to the mesh movement, such as the update of the temperature field to model the heat transfer by convection, are described in [25].

The macroscopic dislocation density dynamics is described by the Alexander–Haasen model [10]. The time evolution of the dislocation density $N_m$ depends on the dislocation velocity v and the effective stress ${\tau }_{\mathrm{eff}}$:

$$\frac{\mathrm{d}{N}_m}{\mathrm{d}t}=Kv{\tau }_{\mathrm{eff}}^l{N}_m,$$

(2)

where K and l are material constants. The dislocation velocity v is given by

$$v=k_0{\tau }_{\mathrm{eff}}^{m}exp(-\frac{Q}{k_{B}T})$$

(3)

with material constants $k_0$ and m, activation energy (Peierls potential) Q, and the Boltzmann constant $k_B$. The effective stress ${\tau }_{\mathrm{eff}}$ is computed from the stress $\tau$ and strain hardening ${\tau }^{*}$ caused by a dense structure of dislocations:

$$\begin{array}{ccc}\hfill {\tau }_{\mathrm{eff}}& =\tau -{\tau }^{*}& \mathrm{if}\tau >{\tau }^{*}\hfill \\ \hfill {\tau }_{\mathrm{eff}}& =0& \mathrm{if}\tau \le {\tau }^{*}.\hfill \end{array}$$

(4)

${\tau }^{*}$ is given by ${\tau }^{*}=D\sqrt{N_m}$ with the strain-hardening factor D. ${\tau }^{*}$ is only relevant for large dislocation densities. $\tau$ is defined by the second invariant of the stress deviator

$$\tau =\sqrt{J_2},\mathrm{with}{J}_2=\frac{1}{2}\sum _{i,j}{S}_{ij}^2.$$

(5)

$\sqrt{3J_2}$ is the von Mises stress. The stress deviator is defined via: $S_{ij}={\sigma }_{ij}-{\delta }_{ij}\frac{1}{3}{\sum }_k{\sigma }_{kk}$. Here, ${\delta }_{ij}$ is the Kronecker symbol and ${\sigma }_{ij}$ is the stress tensor.

In addition to the dislocation multiplication law (2), the Orowan relation in isotropic formulation is used in the AH model to describe the plastic strain components ${\epsilon }_{ij}^c$:

$$\frac{\mathrm{d}{\epsilon }_{ij}^c}{\mathrm{d}t}=\frac{bv{N}_m}{2\sqrt{J_2}}{S}_{ij},$$

(6)

where b is the magnitude of the Burgers vector. Non-zero homogeneous dislocation density distribution and zero plastic strain are set at $t=0$. For additional details regarding the AH model, we refer the reader to our prior works [25,34].

3. Material Parameters

3.1. Germanium

The physical properties for the melt are given in

Table 1. Physical parameters of the Ge melt.

Density	${\mathit{\rho }}_{\mathit{m}}$	$\mathrm{kg} {\mathrm{m}}^{-3}$	5534	${}^{*}$
Heat capacity	$c_p$	$\mathrm{J} \mathrm{k}{\mathrm{g}}^{-1} {\mathrm{K}}^{-1}$	358 (@ $T_m$)	[35]
Thermal expansion	${\alpha }_{\mathrm{V}}$	${\mathrm{K}}^{-1}$	$1.06\times {10}^{-4}$	[36]
Thermal conductivity	$\lambda$	$\mathrm{W} {\mathrm{m}}^{-1} {\mathrm{K}}^{-1}$	48 (@ $T_m$)	[37]
Emissivity	$ϵ$		0.2	[38]
Electrical conductivity	${\sigma }_{\mathrm{el}}$	$\mathrm{S} {\mathrm{m}}^{-1}$	$1.4\times {10}^6$	[39]
Dynamic viscosity	$\mu$	$\mathrm{Pa}\mathrm{s}$	$6.8\times {10}^{-4}$ (@ $T_m$)	[40]
Prandtl number	Pr		0.05
Marangoni coefficient		$\mathrm{N} {\mathrm{m}}^{-1} {\mathrm{K}}^{-1}$	$7.32\times {10}^{-5}$	[41]
Latent heat	$\Delta H$	$\mathrm{J} \mathrm{k}{\mathrm{g}}^{-1}$	$4.65\times {10}^5$	[42]
Laplace constant	$a_L$	m	$5.3\times {10}^{-3}$

${}^{*}$ Intermediate value between measurements of [35,40].

. We used a direct simulation of the melt flow and applied the steady-state Navier–Stokes equation without a turbulent model. For this reason we used an increased viscosity ($\mu =3\times {10}^{-2}\mathrm{Pa}\mathrm{s}$), which ensures a convergence of the steady-state Navier–Stokes solution. The melting point is $T_m=1210 \mathrm{K}$.

For the crystal, we used temperature dependent parameters if possible, because in the upper part of the crystal the temperature becomes lower than 700 K at later growth stages. All parameters for the Ge crystal are given in

Table 2. Physical parameters of the Ge crystal.

Density	${\mathit{\rho }}_{\mathit{m}}$	$\mathrm{kg} {\mathrm{m}}^{-3}$	5370 (@$\mathit{T}=300 \mathrm{K}$)	${}^{†}$
Heat capacity	$c_p$	$\mathrm{J} \mathrm{k}{\mathrm{g}}^{-1} {\mathrm{K}}^{-1}$	418 (@ $T_m$)	[42]
Thermal expansion	${\alpha }_{\mathrm{L}}$	${\mathrm{K}}^{-1}$	$\left(7.6-4.0e^{-T/365\mathrm{K}}\right)\times {10}^{-6}$	${}^{*}$
Thermal conductivity	$\lambda$	$\mathrm{W} {\mathrm{m}}^{-1} {\mathrm{K}}^{-1}$
		$T\le 300\mathrm{K}$	$4.669\times {10}^4{\left(1{\mathrm{K}}^{-1}T\right)}^{-1.16}$	${}^{*}$
		$30\mathrm{K}<T\le 1000\mathrm{K}$	$2.614\times {10}^4{\left(1{\mathrm{K}}^{-1}T\right)}^{-0.729}$	${}^{*}$
		$1000\mathrm{K}<T$	17	${}^{*}$
Emissivity	$ϵ$		0.55	[43]
Electrical conductivity	${\sigma }_{\mathrm{el}}$	$\mathrm{S} {\mathrm{m}}^{-1}$	$e^{15.255-4317.3\mathrm{K}/T}$	[44]
Elastic constants	$C_{11}$	Pa	$(13.2\times {10}^{10}-1.7\times {10}^7{\mathrm{K}}^{-1}T)$	${}^{*}$
	$C_{12}$	Pa	$(4.3\times {10}^{10}-4.2\times {10}^6{\mathrm{K}}^{-1}T)$	${}^{*}$
	$C_{44}$	Pa	$(7.1\times {10}^{10}-1.07\times {10}^7{\mathrm{K}}^{-1}T)$	${}^{*}$
Peierls energy	Q	eV	1.62	[45,46]
Prefactor in Equation (2)	K	$\mathrm{m} {\mathrm{N}}^{-1}$	$4.0\times {10}^{-3}$	[12]
Prefactor in Equation (3)	$k_0$	$\mathrm{m} {\mathrm{s}}^{-1} {\mathrm{Pa}}^{-1}$	$6.9\times {10}^{-3}$	[12]
Exponent in Equation (2)	l		1	[12]
Exponent in Equation (3)	m		1	[12]

${}^{*}$ See Appendix A. ${}^{†}$ for historical reasons we used a slightly larger value than given by [47]: 5330 $\mathrm{k}\mathrm{g} {\mathrm{m}}^{-3}$. For transient computations, only ${\rho }_m{c}_p$ is relevant. For steady-state calculations the density does not matter at all.

. More details on the origin of the functions given in the table are presented in Appendix A.

The parameters for the Alexander–Haasen model (last section in Table 2) are derived from the discrete dislocation dynamics simulations by Gradwohl et al. [12], except the Peierls energy, which was taken from experiments [45,46].

3.2. Other Materials

The physical parameters for steel and quartz were taken as given in [9]. The quartz was treated as completely transparent for the thermal radiation. The data for the graphite susceptor are given in

Table 3. Physical parameters used for graphite.

Density	$\mathit{\rho }$	$\mathrm{kg} {\mathrm{m}}^{-3}$	1750
Heat capacity	$c_p$	$\mathrm{J} \mathrm{k}{\mathrm{g}}^{-1} {\mathrm{K}}^{-1}$	527
Emissivity	$ϵ$		0.8
Thermal conductivity	$\lambda$	$\mathrm{W} {\mathrm{m}}^{-1} {\mathrm{K}}^{-1}$	$90(1.29-1.09\times {10}^{-3}{\mathrm{K}}^{-1}T$
			$+4.016\times {10}^{-7}{{\mathrm{K}}^{-2}\mathrm{T}}^2-4.77\times {10}^{-11}{\mathrm{K}}^{-3}{T}^3)$	[48]
Elec. conductivity	${\sigma }_{\mathrm{el}}$	$\mathrm{S} {\mathrm{m}}^{-1}$
		$T\le 1136\mathrm{K}$	${10}^6/(8.4+{(T-1173\mathrm{K})}^2/1.9\times {10}^5{\mathrm{K}}^2)$	[48]
		$T>1136\mathrm{K}$	${10}^6/(8.95+0.65\tanh ((T-1673\mathrm{K}/350\mathrm{K})))$	[48]

.

For growing HPGe crystals, hydrogen is used as the gas atmosphere (typical pressure is $P=0.1\mathrm{MPa}$). The parameters were set as shown in

Table 4. Physical parameters of H${}_2$ gas.

Density	$\mathit{\rho }$	$\mathrm{kg} {\mathrm{m}}^{-3}$	0.0617 (@$\mathit{P}=0.1\mathrm{MPa}$, $\mathit{T}=400 \mathrm{K}$)
Heat capacity	$c_p$	$\mathrm{J} \mathrm{k}{\mathrm{g}}^{-1} {\mathrm{K}}^{-1}$	14,500 (@$P=0.1 \mathrm{MPa}$, $T=400 \mathrm{K}$)
Heat conductivity	$\lambda$	$\mathrm{W} {\mathrm{m}}^{-1} {\mathrm{K}}^{-1}$	$0.08+3.9\times {10}^{-4}{\mathrm{K}}^{-1}T$	[49]
			(@$P=0.25-0.4\mathrm{MPa}$)
Thermal expansion	${\alpha }_L$	${\mathrm{K}}^{-1}$	$\frac{1}{3}{T}^{-1}$	ideal gas
Dynamic viscosity	$\mu$	$\mathrm{Pa}\mathrm{s}$	$1\times {10}^{-3}$	[50]
			(@$P=0.1\mathrm{MPa}$, $T=350 \mathrm{K}$)

.

4. Experimental Setup

HPGe crystals are grown in equipment with inductive heating. The crucible holder is made of graphite and acts as a susceptor. All other parts in the inner region of the furnace consist of high-purity fused silica (referred to as quartz). A general description is given elsewhere [22] and a sketch is provided in Figure 1. The equipment is used to grow 3″ crystals in the [001] direction using a quartz crucible for the melt. Crystals are pulled with a constant pulling velocity of $v_{\mathrm{pull}}=0.5 \mathrm{mm} {\min }^{-1}$, except when the end cone is grown to empty the crucible. The position of the crucible is fixed during the entire process of crystal growth. This means that due to the decreasing melt level, the growth velocity becomes larger than the pulling velocity. The diameter of the crystal is controlled manually by adapting the heating power. In this paper, we consider an experiment with an ingot mass of 3091 g. During the growth experiment, the rotation rates of the crystal and crucible were 17 rpm and −3 rpm, respectively (counter rotation). The final crystal is shown in Figure 2. We use the shape of this crystal for our calculations. The shape was extracted from the projection on millimetre graph paper, as indicated by the black line at the upper part of the crystal in Figure 2. The two measured diameters (white lines in Figure 2) acted as a reference.

The dislocation density along the length of the crystal was investigated by EPD analysis according to the IEEE Standard Test Procedures for High-Purity Germanium Crystals for Radiation Detectors [51]. For this purpose, three oriented $\left(001\right)$ wafers were cut from the crystal, as indicated by the thick orange lines in Figure 2. The wafers were chemo-mechanically polished and etched. A $\left(001\right)$ HPGe dislocation etch was prepared, consisting of a parent solution of 25 g Cu${\left({\mathrm{NO}}_3\right)}_2$$·$3H${}_2$O in 500 mL deionized water ( 17 $\mathrm{M}$$\Omega$), 250 mL HNO${}_3$ (70% per weight), and 500 mL HF (49% per weight). The wafers were immersed in the etch for 90 s, so that distinct dislocation etch pits were formed. The etch pits were counted and analysed using optical microscopy and differential interference contrast.

5. Results and Discussion

5.1. Global Simulation (Elmer)

For HPGe, no lateral photovoltage scanning (LPS) measurements [52] are possible and, therefore, the interface deflection during growth is unknown. One way to compare the simulation with the experiment is the change in required power over time. In the computation, the power is measured as the power of the electrical current in the coils. The power is adapted to ensure the melting point temperature at the three-phase point. Please note that in the axisymmetric calculation, the coil is represented by rings—the number of rings corresponds to the number of windings in the real coil. In addition, all additional parts, such as current supplies, are not considered in the simulation. Furthermore, power losses in the generator are not considered in the model. Therefore, the absolute power is different in the experiment and simulation. Assuming that this difference is constant in the range of power change during the growth process, we can compare this change. Recently, investigations in model experiments on this problem have given a more detailed picture on power losses in induction heating systems [53], but there is not yet a general rule to be applied to our experiment.

In Figure 3, the change in power $\Delta \mathcal{P}$ is shown for three numerical runs and the experiment. The reference for $\Delta \mathcal{P}=0$ is the minimal power in the experiment (at a time of 346 min) and in the simulations (at a time of 368 min). The black curve (computation without gas convection) is only slightly different from the two others. Thus, gas convection has only a weak influence on the power consumption. Hydrogen, as the high-purity gas, has quite a high thermal conductivity and, therefore, there is already a significant thermal transport by diffusion. When argon is used, gas convection plays a much stronger role. We also checked a different time stepping, i.e., using $\Delta l=2\mathrm{mm}$ in one case (red curve) and $\Delta l=1\mathrm{mm}$ in the other (orange curve). Nearly no difference can be seen in Figure 3.

In Figure 4 (left), the drop in the melt height $h_{\mathrm{melt}}$ is shown, which causes a higher growth rate than the pulling velocity since the crucible position is fixed. On the right-hand side, the increase in the growth velocity $v_{\mathrm{growth}}$ with increasing crystal length can be seen. There are only small differences in $h_{\mathrm{melt}}$ and $v_{\mathrm{growth}}$ between the runs with $\Delta l=2\mathrm{mm}$ and $\Delta l=1\mathrm{mm}$.

At this point, we would like to mention the details of the volume computation for the melt and crystal, which is performed numerically. In both cases, we sum up the volumes of horizontal trapezoidal slices of height $\Delta h_X$ and $\Delta h_m$, for the crystal and melt, respectively. For $\Delta l=2\mathrm{mm}$ we had $\Delta h_X=1\mathrm{mm}$ and $\Delta h_m=1\times {10}^{-3}\mathrm{mm}$, and for $\Delta l=1\mathrm{mm}$ we chose $\Delta h_X=0.5\mathrm{mm}$ and $\Delta h_m=1\times {10}^{-6}\mathrm{mm}$. $\Delta h_m$ has to be quite small because the melt level is decreasing only very little during necking in the beginning of the growth process. A small value of $\Delta h_m$ ensures a smooth increase in the growth velocity. At later times, the increase in the growth velocity is quite significant (by up to 45%) and leads to a larger production of latent heat. In particular, we use $v_{\mathrm{growth}}$ computed at crystal length l for the latent heat calculation of the computation at $l+\Delta l$. For completeness, we also show the time step $\Delta t$, which is computed from the obtained growth velocity and the fixed increase in crystal length. We obtained the upper curve for $\Delta l=2\mathrm{mm}$ (blue curve) and the lower for $\Delta l=1\mathrm{mm}$ (cyan curve). The time step, even for $\Delta l=1\mathrm{mm}$, is much larger than required for the computation of the dislocation density in MACPLAS.

The solid/liquid interface shapes calculated by the global model are plotted at the bottom of Figure 5. The corresponding deflection is plotted at the top of Figure 5. Initially, the interface is almost flat. During the shoulder growth, the interface becomes more and more convex. This is at least partly caused by the melt convection, as can be seen from Figure 6 for the velocity field and Figure 7 for the temperature field. The inductive heating causes the hottest region to be in the right lower part of the melt. As a result a convection roll is established, transporting the heat from this right lower part to the three-phase point of crystal, melt, and gas. The rotation of the crystal prevents this convection roll entering the region beneath the crystal, leaving this region with cooler melt. With increasing crystal and decreasing melt level, the buoyancy-induced convection roll becomes weaker and the interface becomes flatter. With the increasing melt level, the counter vortex in the upper outer part of the melt disappears. At $l=50 \mathrm{mm}$ in Figure 6 this can be clearly seen, leading to a cold region in the top outer part of the melt (see Figure 7). Note that calculations with neglected melt flow were not performed, so that some effects of the changing global temperature field on the interface deflection cannot be excluded.

Looking in more detail at the deflection in the top part of Figure 5 one sees some oscillations in the deflection curve, less pronounced for $\Delta l=1\mathrm{mm}$ than for $\Delta l=2\mathrm{mm}$. In both cases, the mesh size in crystal and melt was $1\mathrm{mm}\times 1\mathrm{mm}$. In addition to questions of spatial resolution, it should be mentioned that we performed quasi-stationary calculations that do not take into consideration the change interface shape from one time step to the next for the computation of the latent heat release at the interface. We had only the latent heat production from the computed overall growth velocity $v_{\mathrm{growth}}$, which is the same all along the interface.

5.2. Local Simulation (MACPLAS)

In MACPLAS we used the crystal shape, including the computed solid/liquid interface shape, as the geometrical input, and use the temperature at the rim of the crystal as the thermal boundary condition. As a check, we compared the stresses (isotropic) for a crystal of length 228 mm. For this purpose a steady-state temperature computation was performed with MACPLAS. The elastic stress calculated by Elmer and MACPLAS was found to be in excellent agreement.

At this point, we have to discuss the stress–strain relation in the axisymmetric case in more detail. For this case Lambropoulos derived modified matrices of elastic constants for two growth directions, $\langle 001\rangle$ and $\langle 111\rangle$ [54]. The typical growth direction for HPGe crystals is $\langle 001\rangle$ and this was also used in our experiment. For this growth direction, the stress–strain relation is

$$\left(\begin{array}{c}{\sigma }_{rr}\\ {\sigma }_{\varphi \varphi }\\ {\sigma }_{zz}\\ {\sigma }_{rz}\end{array}\right)=\left(\begin{array}{cccc}{C}_{11}+\frac{1}{4}H& C_{12}-\frac{1}{4}H& C_{12}& 0\\ C_{12}-\frac{1}{4}H& C_{11}+\frac{1}{4}H& C_{12}& 0\\ C_{12}& C_{12}& C_{11}& 0\\ 0& 0& 0& C_{44}\end{array}\right)\left(\begin{array}{c}{\epsilon }_{rr}-{\epsilon }_T\\ {\epsilon }_{\varphi \varphi }-{\epsilon }_T\\ {\epsilon }_{zz}-{\epsilon }_T\\ {\epsilon }_{rz}\end{array}\right),$$

(7)

where H is given by $H=2C_{44}+C_{12}-C_{11}$, i.e., $H=(5.3\times {10}^{10}-8.6\times {10}^6{\mathrm{K}}^{-1}T)\mathrm{Pa}$ for Ge. ${\sigma }_{ij}$ and ${\epsilon }_{ij}$ are the stress and strain components, respectively. The thermal strain is ${\epsilon }_T=\overline{\alpha }(T-T_{\mathrm{ref}})$, where $T_{\mathrm{ref}}$ is the reference temperature and $\overline{\alpha }$ is the averaged thermal expansion coefficient (calculated from the differential expansion coefficient ${\alpha }_L$; see [34]).

Currently, only the isotropic model with $H=0$ (i.e., the influence of the growth direction is not considered) is implemented in MACPLAS:

$$\left(\begin{array}{c}{\sigma }_{rr}\\ {\sigma }_{\varphi \varphi }\\ {\sigma }_{zz}\\ {\sigma }_{rz}\end{array}\right)=\left(\begin{array}{cccc}{C}_{11}& C_{12}& C_{12}& 0\\ C_{12}& C_{11}& C_{12}& 0\\ C_{12}& C_{12}& C_{11}& 0\\ 0& 0& 0& C_{44}\end{array}\right)\left(\begin{array}{c}{\epsilon }_{rr}-{\epsilon }_T-{\epsilon }_{rr}^c\\ {\epsilon }_{\varphi \varphi }-{\epsilon }_T-{\epsilon }_{\varphi \varphi }^c\\ {\epsilon }_{zz}-{\epsilon }_T-{\epsilon }_{zz}^c\\ {\epsilon }_{rz}-{\epsilon }_{rz}^c\end{array}\right),$$

(8)

where the plastic strain components ${\epsilon }_{ij}^c$ are obtained by integrating Equation (6). Elastic simulations can be carried out by setting ${\epsilon }_{ij}^c=0$. This means that all the results presented later for the dislocation densities are based on the stress computation with $H=0$. In order to estimate the range of error with this approach, i.e., to evaluate the importance of the anisotropy, we checked the impact of the correction term H on the stress field in a steady-state Elmer computation for a crystal of length 228 mm. The stress distribution is similar in both cases, as can be seen in Figure 8. However, the stress is about 15% larger for the case $H\ne 0$.

Using the crystal shape and temperature field at the crystal side surface calculated by Elmer, fully transient temperature and dislocation density simulations were carried out using MACPLAS. Since the maximum time step was set to $0.2$ $\mathrm{s}$ (about three orders of magnitude smaller than the temporal resolution in the Elmer model; see Figure 4), the original Elmer data were interpolated in time, as described in Section 2.2.

The calculated dislocation density in the crystal is shown in Figure 9. The initial density in the entire seed was set to $N_0=1 \mathrm{c}{\mathrm{m}}^{-2}$. At the final time, the maximum dislocation density is three orders of magnitude higher than the initial one. The highest dislocation density is at the crystal side surface; a local $N_m$ maximum develops on the crystal axis, thus, the radial dislocation density distribution has a W-shaped profile (see also [25]). The W-shape can be even better seen in Figure 10, where the radial distribution of the dislocation density is shown for lines at three different vertical positions in the final crystal. These positions correspond to the horizontal cuts for the wafers from the crystal that was grown, and are indicated in Figure 2.

The computed dislocation densities are compared to the EPD data on the three wafers. EPD data were taken at the centre and near the edge of the wafer. As can be seen from Figure 10 and Figure 11, a good match is obtained for the first cut ($l=80 \mathrm{mm}$). For the other positions, the dislocation density in the experiment was higher. This discrepancy, even to a greater extent, was already observed in another investigation [25]. We also see in Figure 11 that for the dislocation densities at $l>180\mathrm{mm}$ there is a significant difference between the runs with $\Delta l=2 \mathrm{mm}$ (dashed line) and with $\Delta l=1 \mathrm{mm}$ (solid line). This indicates that small differences in the interface shapes might lead to a larger difference in dislocation densities at later times. From Figure 5, one can see that the differences are in the sub-millimetre range. Thus, we need a very fine mesh in the region of the interface to reach a high accuracy. Another issue is the computation of the melt convection. Higher accuracy might be achieved by coupling the global simulation of Elmer with the computation of convection by a finite volume code such as OpenFOAM (https://www.openfoam.com (accessed on 27 September 2023)), which is better suited for this purpose. A third point on the global simulation is the missing latent heat effect when the interface shape is changing from one time step to the next. This also requires a very fine mesh in the vicinity of the interface, as already mentioned above, and in addition, more inner iteration steps are necessary to converge the coupled system of geometrical change, latent heat production, and melt convection.

There might also be several reasons at the level of the local calculations for the discrepancy between the computed dislocation density and EPD:

• We use a quasi-transient approach for the global thermal calculation.
• For the stress calculation, $H=0$ was used in Equation (7) and not Lambropoulos’ approximation for growth in the $\langle 100\rangle$ direction.
• We used a simple AH model without considering glide planes [55].
• The AH model is a local model in the sense that dislocations do not propagate from one computational node to a neighbouring one.

The last point in the list is important from an experimental point of view. In the end of the growth process, the crystal is detached from the melt or the melt is consumed completely. In both cases, there is a strong tendency of dislocation multiplication at the end of the crystal running back into the bulk volume. Such a propagation of dislocation from the end of the crystal towards the central bulk region is not covered by the Alexander–Haasen model. This could significantly influence the final dislocation density at $l=205 \mathrm{mm}$. In order to study this in detail, a mapping of the etch pit densities over larger areas is required. This would allow a broader comparison between the experimental and numerical results. A better understanding of the dislocation dynamics in the middle of the growth process, both by further developing the numerical approach as well as a more detailed dislocation density distribution study along the crystal, will be the prerequisite for defining possible strategies for dislocation engineering.

6. Conclusions

We have developed a uni-directional coupling of global axisymmetric calculations using the open source package Elmer with local axisymmetric calculations of dislocation densities in the crystal using the software tool MACPLAS. The local simulations apply a series of crystal shapes (including melt/crystal interface) as the geometrical input and the series of temperature distributions at the rim of the crystal as thermal boundary conditions for the calculations. The local simulation is a transient calculation with a much smaller time step than the series of global computations. Therefore, geometrical and thermal boundary conditions are interpolated.

This numerical approach has been applied to the Czochralski growth of high-purity germanium, where inductive heating is used. We performed one experiment for comparison with the computations. The shape of the crystal grown in the experiment was used in the computations. In the global simulation of the temperature field, the heater power was changed by an algorithm to keep the temperature at the three-phase point at melting point temperature. This change in power was in good agreement with the change in the power in the experiment.

The etch pit density (EPD) was measured for three wafers from different parts of the crystal. The computed dislocation density at the top part (shoulder) agrees well with the EPD. In the further growth phases, the dislocation density in the simulation does not increase so much as in the experiment and, thus, we have the largest disagreement for the wafer at the end of the cylindrical part. There might be several reasons for this discrepancy. On the level of global simulations the calculation of the melt convection, i.e., crystallisation interface shape should be improved. In the computation of dislocation densities there are several options to extend the simple Alexander–Haasen model. A more systematic study, with more experiments exploiting the EPD at several locations, will be important and will be carried out as a next step.