The possibility of growing crystals with homogeneous impurity distribution over crystal height has been demonstrated in a study of segregation during silicon and germanium growth from thin melt layers using the submerged heater method. Numeric simulation of 200 mm diam. antimony-doped germanium crystallization has shown that, beginning from a 40 mm melt layer thickness, the exact problem solution with convection allowance is identical to the unidimensional heat exchange problem solution in the central ingot part. The conditions under which convection can be ignored in mass transport calculation are more rigorous: the melt layer height must be within 20 mm. In this case Tiller’s ratio can be used for calculating the longitudinal impurity distribution for predominantly diffusion-controlled mass transport pattern. Analysis of the existing attempts to describe the experimental crystal growth results using the simplified formulae shows that they only yield acceptable results if the actual growth rate or change in melt layer thickness during crystallization are taken into account, e.g. as in the formula suggested by Marchenko et al. One can therefore analytically describe the longitudinal impurity distribution in the ingot, e.g. B and P distribution in silicon, and recommend the degree of additional doping of the melt zone under the heater so that to provide a constant impurity concentration over the ingot height. Homogeneous material can be obtained after residual layer solidification in the end portion of the ingot if the growth rate is controlled through varying the cooling rate.

mathematical simulation, impurity diffusion, directional crystallization, submerged heater, germanium, multicrystalline silicon

The directional crystallization method used for the production of multicrystalline silicon does not provide a homogeneous doping impurity distribution along the growing crystal. Detailed analysis [1] showed that the concentration of impurities such as boron and phosphorus which is described adequately well by Shell’s ratio grows slightly at an early growth stage, accumulates during the crystal growth and increases abruptly in the end portion of the ingot. This latter zone is typically 15–20 % of the ingot length which makes at least 30–50 mm of the material at the top of the ingot low-quality. Impurity concentration change along the ingot affects the resistivity of silicon and is therefore detrimental but this problem cannot be solved with standard methods. However, the task of producing a height-homogeneous ingot can be solved by using the melt-submerged heater technique for directional crystallization growth which is implemented either as the submerged heater method (SHM, [2]) or the axial heat flow method (AHP, [3]). The heater is placed at a short distance from the phase boundary and partitions the melt zone in two (the crystallization zone under the heater and the feeding zone above the heater), thus fundamentally changing the stream pattern. Analysis of the prerequisites for the technical implementation of this approach to producing multicrystalline silicon using the AHP method was reported earlier [4]. The abovementioned paper [4] reported results of testing a protective coating [5] of submerged heater casing for silicon crystallization from a standard G1 crucible.

Crystal growth from a thin melt layer significantly reduces convection intensity [6] and provides for a segregation pattern that is close to a diffusion one [7]. It was shown with assumptions for experimentally grown Te-doped GaSb [8] that even for crystals grown on Earth the impurity concentration distribution C_S (x) can be described adequately well by an expression derived for the specific implementation conditions, i.e., the so-called diffusion controlled mode [9]:

(1)

where C₀ is the initial impurity concentration in the melt, V is the growth rate, k is the equilibrium segregation ratio, D is the diffusion coefficient and x is the coordinate along the crystal axis.

Segregation in a thin melt layer reaches an equilibrium condition considerably faster [10], and under specific conditions the effective segregation ratio k_eff becomes unity [1]. This process was analytically described earlier [12] with the following formula:

(2)

It is assumed herein that different melt zones (under the heater and above the heater) had the impurity concentrations C₁ and C₂, respectively, C_S is the initial impurity concentration in the melt, and h is the melt layer thickness. Obviously, if the feeding zone is doped to the concentration C₂ = C_S and the growth zone is doped so that the ratio C₁ = C₂/k is true, the crystal will have a constant impurity concentration equal to the impurity concentration in the feeding zone. Unfortunately, at the final stage of crystallization when there is no more fresh melt in the feeding zone and the submerged heater does not move relative to the crucible, Eq. (2) is no longer valid. However, nor true is the earlier suggested ratio [13] for the complete mixing case since for the SHM or the AHP methods the submerged heater is not removed from the melt even after crucible pulling is stopped. There is no free surface and hence no Marangoni convection, and the diffusion contribution to mass transport remains prevailing.

Thus there is a need for a numeric simulation to study longitudinal segregation under the conditions set forth above and for determining the optimum crystallization heating mode providing for the growth of a homogeneous crystal over the maximum possible ingot distance. This study will determine the applicability limits of the impurity distribution analytical expressions to be further used for engineering calculations of silicon growth. Though only the “residual” melt portion is in question, the length of this zone may be as large as 100 mm or greater if large-size multicrystalline silicon ingots are grown because the melt layer height for the submerged heater method is chosen based on the h/d ~ 0.15 ratio. The model material for this study was antimony doped germanium since the antimony segregation ratio is 0.003, i.e., a priori lower than any of the segregation ratios of commercially important impurities in silicon, and furthermore the problem in hand becomes the most complex for small k. Moreover, the properties of germanium have been well studied and there is the possibility to conduct, if necessary, experiments at lower temperatures than would be required for silicon.

For submerged heater crystallization, by analogy with the conventional Bridgman method, the crucible is lowered into the cold zone of the growth chamber (Fig. 1a) and the crystal is grown until the melt in the zone above the heater is over. The unidimensional heat transfer problem can be described with the following equation:

Figure 1.

Crystal growth setup at (a) crucible pulling stage relative to submerged AHP heater and (b) at the end after stopping, and (c) domain schematic.

(3)

where λ_L,S is the liquid and solid state heat conductivity, ρ_L is the melt density, J is the crystallization heat, T_m is the melting point, T_hot and T_cool are the temperatures at the hot and cool boundaries of the melt/crystal system at the AHP heater bottom and the crucible, respectively, and H is the height of the grown crystal. To ensure constant crystallization conditions at the phase boundary, including constant axial temperature gradient, one should maintain h and T_hot constant. Meanwhile the cool boundary temperature is gradually reduced; for nontransparent melts, e.g. germanium and silicon, the reduction obeys a linear law in accordance with Eq. (3).

The heater is typically removed from the melt at the end of the process to avoid its freezing into the crystal, so the final portion crystallizes relatively rapidly and is later cut off from the ingot. Below we will consider a problem (Fig. 1b) for which crystal pulling is already stopped but the submerged heater is still in contact with the melt. This configuration allows controlling crystallization until its completion thus maintaining the required temperature at the top T_top(t) and bottom T_bot(t) boundaries of the domain. It can be seen from Fig. 1 c that the domain for an axially symmetrical configuration is a half-disc with the diameter d and the height Δ which is cooled through simultaneous reduction of the temperatures at its ends, with the cooling rate being time-variable in the general case. The initial impurity distribution in the melt is accepted to be homogeneous over the entire domain and is set to unity.

We used the ANSYS Fluent software module for calculating crystallization and melting with a constant grid [14]. This module is designed for melts having a solidification zone between the solid and the liquid phases. To justify the use of this approach in this work, we considered a small temperature range (1 K) around the melting point, and therefore the solidus T_S and liquidus T_L temperatures were set as follows:

T _S = T_m – 0.5; T_L = T_m + 0.5.

The specific mixture enthalpy h_mix expression can be written as follows:

where c_S and c_L are the specific heat capacities of the solid and the liquid phases, J is the latent melting heat and f_L is the liquid phase fraction which is

The energy conservation equation will be written as follows:

(4)

Equation (4) is solved taking into account a linear temperature change at the bottom T_bot(t) = T_bot(0) – R_coolt and the top T_top(t) = T_top(0) – R_coolt surfaces, respectively. We assume that cooling has the constant rate R_cool. The convection heat exchange at the side boundary is

where α is heat convection coefficient.

The momentum conservation equation is solved using the porosity enthalpy method which considers the solidification zone as a porous medium. The momentum absorption term added to Eq. (4) as a flow for bringing the rates in the solid phase zones to zero becomes as follows:

where ɛ = 0.001 is the small number used for avoiding division by zero and A_mush is the solidification zone constant.

The properties of germanium (density, heat capacity and heat conductivity) used for the calculations are shown in Table 1 with allowance for their stepwise change at the phase boundary:

Table 1.Download as CSV 

Physical properties of germanium and silicon used in calculations.

Properties	Notation	Germanium [15–17]	Silicon [18,19]
Melting point, K	Tm	1210	1683
Density, kg/m3:
Crystal	ρS	5260	2330
Melt	pL	5550	2530
Latent melting heat, kJ/kg	J	460	164
Heat capacity, J/(kg × K):
Crystal	C S	400	700
Melt	C L	1140	1100
Heat conductivity, W/(m × K):
Mass specific heat capacity, J/(kg × K)
Crystal	λS	17	22
Melt	λL	39	67
Kinematic viscosity, m2/s	ν	1.35 × 10-5	1.5 × 10-5
Heat expansion coefficient, K-1	β	0.97 × 10-4	1.32 × 10-4
Segregation ratio	k
Sb		0.003	–
P		–	0.35
B		–	0.8
Diffusion coefficient, cm2/s	D
Sb		1 × 10-4	–
P		–	3.3 × 10-4
B		–	2.7 × 10-4

(5)

The impurity distribution is set in a similar manner using mixture concentration description:

where the generalized scalar transfer equation

(6)

allows calculating the impurity distribution by introducing the UDS user defined function for the concentration. The transfer equation for the solidification zone imitating the crystallization front allows for a segregation-related impurity source which was calculated as follows:

(7)

The parameters for which the numeric calculations were carried out are summarized in Table 2.

Table 2.Download as CSV 

Numeric calculation parameters.

Parameter	Notation	Germanium	Silicon
Heat convection coefficient, W/(m2 × K)	α	25	25
Ambient temperature, K	T∞	800	1300
Cooling rate, K/h	R cool	10–60	15–40
Initial impurity concentration in eq. (1), cm-3	С 0
Sb		2.7 × 1017	–
P		–	6.0 × 1017
B		–	2.62 × 1017
Initial impurity concentration in Eq. (2) above AHP Heater, cm-3	С 2
P			2.1 × 1017
B			2.1 × 1017
Initial impurity concentration in Eq. (2) under AHP heater, cm-3	С 1
P		–	2.62 × 1017
B		–	6.0 × 1017

Figure 2 shows typical established temperature field patterns in the domain, stream function in the melt and velocity vectors certain time after crystallization onset. Since the top and bottom boundaries of a thin layer (Δ/d < 0.15) are maintained under isothermal conditions, no convection occurs in the middle of the domain. Convection driving force only emerges at the edge, due to the difference between the side boundary temperature and the ambient temperature, and a vertex with a clockwise stream forms at the wall. The bulk temperature field is unidimensional and only undergoes distortion at the edge. The swirl only exists at the side wall and can be considered to be localized in the domain having Δ = 10 mm. With an increase in Δ (Fig. 3) the stream penetrates deeper and can even affect the middle part. This however does not influence the temperature profile in most of the domain. Therefore convection can be ignored for heat transfer calculations with such Δ/d.

Figure 2.

(a) Temperature distribution. (b) stream function and (c) velocity vectors at periphery of domain in 500 s after Ge solidification onset.

Figure 3.

Stream velocity (m/s) in the domain at an early crystallization stage for different melt thicknesses: (a) 10 mm, (b) 20 mm and (c) 40 mm.

Interestingly, the stream which forms in quite a short time even for Δ = 40 mm as can be seen from Fig. 4 changes but slightly thereafter, until certain portion of the domain crystallizes and the melt layer starts depleting noticeably. This is well illustrated by Fig. 5 where it can be seen how the vertex affected zone decreases in size, the convection intensity drops and the stream structure degenerates with the growth of the crystal. A completely different pattern is observed if one considers the effect of convection on mass transport. It can be seen from Fig. 4b that the initially homogeneous impurity distribution undergoes distortion in a short time due to convection and continues changing, including in the middle of the domain. It is therefore difficult to conclude about a predominantly diffusion pattern of mass transport, at least for Δ = 40 mm. The situation is completely different if the domain height is 10 mm (Fig. 6). The convection emerging at the side wall causes transfer impurity such that its excess quantity is accumulated immediately behind the vortex where an antimony-enriched crystal grows. Meanwhile closer to the middle of the domain the inhomogeneity is completely eliminated and the concentration isocontours are parallel to the flat phase boundary in the middle part. Thus a diffusion crystallization mode is established in almost the entire zone and the axial antimony distribution in germanium calculated using Eq. (1) will hold in almost the whole ingot. This is confirmed by the axial impurity profile pattern shown in Fig. 7. It can well be seen that the antimony distribution profiles over the height of the domain (in the crystal and in the melt) is almost the same for different crystallization stages until a 75 mm radius is achieved. Only the extreme edge portions exhibit differences in the absolute values and the pattern of the profile.

Figure 4.

Formation dynamics for (a) velocity field and (b) antimony concentration in germanium during crystal growth from 40 mm thick melt layer in 500, 750 and 1000 s after solidification onset. Initial concentration in the domain is unity.

Figure 5.

Stream pattern in germanium melt for (a) ∆ = 10 mm and (b) ∆ = 20 mm and its degeneration with decreasing crystal growth melt layer thickness.

Figure 6.

Change in antimony relative concentration in melt and growing germanium crystal in (a) 500, (b) 1250 and (c) 1750 s after solidification onset for domain with ∆ = 10 mm.

Figure 7.

Axial change in antimony relative concentration in (a, c, e, g) 1000 and (b, d, f, h) 2000 s after solidification onset for different crystal radii: (a, b) 25 mm, (c, d) 60 mm, (e, f) 75 mm and (g, h) 95 mm).

In practice only A.G. Ostrogorsky et al. used Tiller’s ratio to describe the SHM growth experimental results. The best fit with theory was achieved by analyzing the Te distribution in seed grown small-diameter GaSb single crystals [7] and later [8] where a disc-shaped seed with a diameter equal to the crucible diameter was used by analogy with the AHP method. However to achieve a closer matching the authors had to use the diffusion coefficient D = 1 × 10^-5 cm²/s for the calculations as per Eq. (1) which is several times smaller than the conventional value for tellurium [20, 21]. The authors justified their choice by the fact that the diffusion coefficient in the vicinity of the growing crystal within the boundary diffusion layer is far lower compared with that in the liquid bulk.

We however would like to offer a simpler explanation related with uncontrolled melt layer thickness change and hence growth rate in the experiments. This explanation can well be accepted taking into account that for the submerged heater design in question A.G. Ostrogorsky et al. never used a heating element or reduced the furnace temperature but only pulled the crucible to the cool zone within the preset temperature gradient. The latter technical solution was dictated by the fact that there is not thermocouple at the crucible bottom in the SHM method which would provide T_cool temperature readings and thus allow temperature control during crystal growth in accordance with Eq. (3). Therefore just like in many other methods there is no ground to assume that the growth rate (phase boundary velocity) V is equal to the pulling rate P. Although the crucible bottom temperature decreases, crystallization delays by several minutes or a longer time for the following reasons. First, there is a transient process of establishing during which the melt/crystal system has to be sufficiently cooled to a certain new condition. To start crystal growth, it is insufficient to reduce the temperature solely by lowering the crucible within the existing temperature gradient. The temperature should be somewhat lower in order to allow crystallization heat removal. Otherwise, as follows from Eq. (3), the T_cool temperature will decrease while the crystal height H will remain the same. Furthermore one can expect that the temperature gradient established before pulling will be closer to the temperature gradient in the melt since in the initial position [8] the melt zone was ~60 mm while the crystal was only 6 mm thick. In this case, as follows from Eq. (3), the melt layer h will grow with increasing distance between the heater bottom and the crucible bottom. This is illustrated by Table 3 demonstrating for the example of germanium crystallization that even if the furnace temperature gradient is the same as the one in the crystal before growth onset (20 K/cm) the growth will be accompanied by an increase in h to 10.87 mm after pulling for 1 h. Only setting a slightly higher temperature gradient (23.5 K/cm) will keep the layer thickness the same (h = 10.0, see Table 3) which should be ensured by choosing correct crystallization mode. If the furnace temperature gradient is close to the melt temperature gradient, then even with allowance for the crystallization heat the layer thickness will be 11.81 mm after pulling for 1 h (Table 3). An additional factor that is specific for submerged heater growth is the presence of a small but much overheated (by at least 30–50 K) melt portion moving from the zone above the heater to the zone under the heater immediately after the start of pulling. This also contributes to the increase in h. Therefore a general technical recommendation is to achieve the target puling rate gradually in order for the crystal to have sufficient time for cooling and thus to grow in accordance with the preset cooling rate.

Table 3.Download as CSV 

Change in germanium melt layer thickness h as a function of furnace temperature gradient for crucible pulling to cool zone*.

No.	Crystallization Stage	Temperature Reduction under Condition:	gradTc, K/cm	Cooling Zone Gradient, K/cm	T cool, °C	h, mm
1	Before pulling	Crystallization heat is constant	20.0	20.0	927.0	10.00
2	After lowering by 5 mm		20.0	20.0	917.0	10.00
3			20.0	10.0	922.0	11.43
4		Crystallization heat is released but ignored in cooling rate calculations	23.5	20.0	917.0	10.87
5			23.5	10.0	922.0	12.38
6		With account of crystallization heat release	23.5	23.5	913.5	10.00
7			23.5	10.0	920.2	11.81

* The hot boundary temperature T_hot = 947 °C is constant during the growth, the initial melt layer thickness is h₀ = 10 mm, the initial crystal length is H₀ = 5 mm; the temperature gradient is gradT_m = 10 K/cm for the melt and gradT_c = 20 K/cm for the crystal, the cooling rate is R = 5 mm/h; the solid and liquid phase germanium heat conductivities are rounded to 20 and 40 W/(m × K), respectively.

Once the melt is finally cooled to the temperature required for the onset of growth the crystallization front velocity V starts increasing but reaches the preset pulling rate R nonmonotonically. Calculations of V based on preset T_cool and T_hot temperatures in a CsI (Tl) single crystal growth experiment [22] showed that during the time of establishing (until the establishment of a new growth rate) the instantaneous growth rate was several times greater than the growth rates at the target growth rate change from 2 to 4 mm/h due to melt/crystal system cooling. This is also corroborated by our calculation results which will be analyzed below. Thus one can quite reasonably assume that the melt layer in the cited work [8] was not constant but quite probably started decreasing after a certain increase so the crystallization front almost approached the submerged heater bottom and followed it thereafter. Thus the growth rate finally became equal to the pulling rate after passing its peak value. Otherwise the experiment would be stopped due to heater freezing into the crystal. This seems to be the case in earlier work [23] which included unfinished InSb crystal growth experiments under microgravity conditions on board the ISS where a noticeable jump (by approx. 3–4 times) of the heat conductivity coefficient was observed during the solid to liquid state transition [24]. This leads to the formation of significantly different axial temperature gradients in the melt and in the crystal. Crystallization from a very thin layer [24] is suggested by the short time in which a steady-state Te distribution was achieved over the entire 15–20 mm crystal length which evidently testifies to a diffusion mass transport pattern. On the contrary, during growth of similar crystals on Earth [25] the axial temperature gradient in a similar furnace was insufficient for the crystallization front to follow the submerged heater. In all the experiments the melt seemed to occupy the entire ampoule in some time, reaching at least a 30–40 mm height. This melt portion most likely solidified rapidly during furnace cooling (after the end of the experiment) and was polycrystalline as the presented photos suggest. The single crystal regions were within 2–12 mm in length which is comparable with the size of the seed.

With the above assumptions about the variation patterns of the actual crystallization front velocity and the thickness of the melt layer from which the GaSb crystal grew, Curves 3 and 4 (Fig. 8) calculated in accordance with Eqs. (1) and (2) using the standard data on the tellurium diffusion coefficient D = 3 × 10^-5 cm²/s are quite close to Curve 2 calculated for D = 1 × 10^-5 cm²/s in the assumption that V = R. Thus Curves 3 and 4 (Fig. 8) describe earlier experimental data [8] adequately well. In comparison with Curve 2 they better describe the experimental data for the initial crystal part because Eq. (1) yields a convex curve unless the change of V in time is taken into account whereas the experimental curves have a concave bend, just like Curves 3 and 4. Thus, approximations of Eqs. (1) and (2) can be used in impurity distribution calculations for specific shapes of the crystallizing layer, provided however that the actual crystal growth rate and melt layer thickness patterns are taken into account.

Figure 8.

Calculation of Te distribution along GaSb axis using (1–3) Tiller’s ratio (Eq. (1)) and (4) Eq. (2): (1, 2) V = 5 mm/h, D = 3 × 10^-5 and 1 × 10^-5 cm²/s, respectively; (3) D = 3 × 10^-5 cm²/s for V growth from 0 to 15 mm/h and further drop to 5 mm/h; (4) h first reaches 12 mm, then drops to 2 mm when growth rate exceeds pulling rate.

Even if the melt zone under the submerged heater in the AHP method is not intentionally doped, the impurity distribution over the silicon ingot height as shown in Fig. 9 (Curves 1) will differ from that for standard directional melt crystallization because the thin ~ 10–20 mm melt layer achieves a steady impurity concentration relatively rapidly, this occurring much faster for boron with a close to unity segregation ratio than for phosphorus. However if the melt zones under and above the submerged AHP heater are doped as described above, the boron and phosphorus concentrations will be constant over the entire growth process in accordance with Eq. (2) regardless of the melt layer thickness which can be chosen as large as 50–100 mm depending on the crucible size. Certainly these element concentrations will be at their respective levels depending on the initial B and P concentrations in silicon.

Figure 9.

(a, b) boron and (c, d) phosphorus distributions over silicon ingot length: (b, d) initial establishing stage: (1) without special melt zone doping (C₁ = C₂ = C_S/k); (2) with additional doping of zone under AHP heater (C₂ = C_S, C₁ = C₂/k).

It can be seen from Fig. 10 that it is no problem to ensure constant boron concentration in the final portion of crystallizing silicon. Whichever crystallization rate in the 1 to 60 mm/h range is chosen the B concentration will change by 30 at.% worst. The situation is completely different for phosphorus distribution. To achieve at least the same result as for boron one should set a very low crystallization rate, ~1 mm/h. The process will take 50–100 which is almost one order of magnitude longer than the crystallization time of the main ingot part. Increasing the growth rate to 60 mm/h will lead to at least a threefold increase in the impurity concentration, to say nothing about possible structural quality problems of the material. Obviously to satisfy both the conditions, i.e., maintain the P concentration within the desired range and complete crystallization within a reasonable time, one should change the growth rate at a certain process stage. Figure 1 illustrates this possibility: the total crystallization time will be ~15 h if the growth rate is abruptly reduced almost to zero after the onset of crystallization. This can be achieved by changing the cooling rate in time, though this is not such a simple task and cannot be solved simply by replicating the same time function R_cool(t) as that for V. This is illustrated by Fig. 11a which suggests that the growth rate is a nonlinear function of the cooling rate.

Figure 10.

(a) boron and (b) phosphorus distributions in the final silicon portion for (a) 1, (b) 10 and (c) 60 mm/h crystallization rate.

Figure 11.

Typical phase boundary rate V change pattern for (a) germanium melt layer solidification and (b) phosphorus distribution for silicon crystallization at variable rate V = V₀(1 – 0,5x^0.3)

The ANSYS Fluent software package along with the in-house mathematical tools were used for the study of 2D segregation in thin melt layers. We chose antimony doped germanium as simulation object due to the low antimony segregation ratio of 0.003. Natural convection during AHP semiconductor growth can be efficiently suppressed and its effect on diffusion can be made negligible for some flat layer thicknesses, i.e. the so-called diffusion-controlled crystal growth mode can be established. This mode allows describing segregation using simple expressions which yield good results if the changes in the actual growth rate and melt layer thickness during crystallization are taken into account. Additional melt doping under the AHP heater provides for a constant impurity concentration in silicon during crucible pulling, and upon the solidification of the residual layer this can be achieved by varying the growth rate (from high to low) through changing the cooling rate in time. By and large the required cooling rate vs time pattern can be determined by solving the inverse task; this can be the objective of the next work stage.