Possible causes of electrical resistivity distribution inhomogeneity in Czochralski grown single crystal silicon

Electrical resistivity distribution maps have been constructed for single crystal silicon wafers cut out of different parts of Czochralski grown ingots. The general inhomogeneity of the wafers has proven to be relatively high, the resistivity scatter reaching 1–3 %. Two electrical resistivity distribution inhomogeneity types have been revealed: azimuthal and radial. Experiments have been carried out for crystal growth from transparent simulating fluids with hydrodynamic and thermophysical parameters close to those for Czochralski growth of silicon single crystals. We show that a possible cause of azimuthal electrical resistivity distribution inhomogeneity is the swirl-like structure of the melt under the crystallization front (CF), while a possible cause of radial electrical resistivity distribution inhomogeneity is the CF curvature. In a specific range of the Grashof, Marangoni and Reynolds numbers which depend on the ratio of melt height and growing crystal radius, a system of well-developed radially oriented swirls may emerge under the rotating CF. In the absence of such swirls the melt is displaced from under the crystallization front in a homogeneous manner to form thermal and concentration boundary layers which are homogeneous in azimuthal direction but have clear radial inhomogeneity. Once swirls emerge the melt is displaced from the center to the periphery, and simultaneous fluid motion in azimuthal direction occurs. The overall melt motion becomes helical as a result. The number of swirls (two to ten) agrees with the number of azimuthally directed electrical resistivity distribution inhomogeneities observed in the experiments. Comparison of numerical simulation results in a wide range of Prandtl numbers with the experimental data suggests that the phenomena observed in transparent fluids are universal and can be used for theoretical interpretation of imperfections in silicon single crystals.


Introduction
The growth of silicon single crystals with homogeneous impurity distributions using the Czochralski technique requires maintaining growth conditions providing for a flat crystallization front (CF), minimum radial temperature gradients at the crystallization front and absence of solutal undercooling of the melt under the crystallization front [1][2][3]. This problem becomes increasingly complex with an increase in the crystal diameter. The requirements to the homogeneity of the electrophysical properties of the single crystals become increasingly strict because the semiconductor industry transits to submicron technologies [4][5]. Optimization of single crystal growth process modes requires a combined effort including physical and mathematical simulation of the growth processes as well as analysis of the homogeneity of the as-grown crystals. Numerical simulation of single crystal growth processes for the Czochralski technique is currently the focus of intense research for various growth chamber design options [6][7][8][9][10][11][12], yet there are no works that would combine numerical simulation of the hydrodynamic growth conditions with electrical resistivity simulation and measurement in growing crystals. It has been believed that, once a flat crystallization front is achieved, the main cause of radial inhomogeneities is the crystallographic symmetry of planes in the growing crystal. Therefore the radial impurity distribution is expected to contain 4 or 8 characteristic surfaces [13]. The aim of this work is to study the homogeneity of electrical resistivity distribution in Czochralski grown silicon single crystals and to analyze the causes of various inhomogeneity types on the basis of crystal growth physical simulation data.

Materials and experimental
We studied electrical resistivity distribution maps for 17 wafers cut out at the distance X from the beginning of five boron-doped p-type 1 Ohm·cm (KDB 1) Grade single crystal silicon ingots grown using the Czochralski technique in REDMET type growth plants [14][15].
The wafer diameter was 154 ± 2 mm, the thickness being 2.4 ± 0.2 mm. The measurements were conducted on a VIK UES-A instrument using the four-probe method with a linear probe setup. Electrical resistivity calculations were carried out taking into account all the corrections as per the applicable standards [16][17] for accurate evaluation of electrical resistivity distribution across the wafers. The overall measurement error in the 1-10 Ohm × cm range as determined in the course of calibration by the Standard Authority was within 1.3%, the random measurement error being 0.5%. The numbers of measurement points on the surfaces were 49 and 121. Single-factor dispersion analysis data were used for verification and confirmed with a high probability the hypothesis of non-random electrical resistivity distribution for all the test wafers. Table 1 shows calculated average electrical resistivity (r av ), RMS deviation (s) and maximum scatter d = r max -r min as percentage of r av . The thermophysical conditions of the growth were analyzed using results of physical simulation of Czochralski crystal growth for a growth plant the parameters of which were reported earlier [18]. Table 1 shows electrical resistivity distribution map parameters for the 17 test wafers showing distances from ingot beginning (x) at which the wafers were cut off, and the respective inhomogeneity types. Except for one wafer cut from the ingot beginning at the edge of the growth cone, all the wafers were quite homogeneous (the standard deviation was within 2.8%) but two inhomogeneity types could be traced from their electrical resistivity distribution patterns: azimuthal and radial. The azimuthal electrical resistivity inhomogeneity was distinguished by the presence of several distribution peaks at the outer radius side along the wafer circumference ( Fig. 1). Table 1 shows the number of the peaks observed. The radial inhomogeneity was expressed by the asymmetrical hill in the electrical resistivity at one wafer side (Fig. 2). Most of the specimens having radial electrical resistivity inhomogeneity also exhibited traces of quite a symmetrical azimuthal electrical resistivity distribution. The inhomogeneity types are also shown in Table 1. The number of electrical resistivity hills that can be resolved in the 3D map is also shown for the radial distribution. Specimen No. 8 cut from the beginning of the ingot is differed from the other specimens by a dramatic drop in the electrical resistivity at 20 mm from the wafer edge ( Fig. 3), which caused a large scatter of the electrical resistivity (28.8%) and RMS deviation (9.2%) compared with the other wafers. A single ingot could exhibit both single-pattern electrical resistivity distribution maps and variable map patterns. The latter observation suggests variation of the growth parameters during crystal pulling from the melt.

Results and discussion
The general tendency of the electrical resistivity map evolution can be accounted for based on the results of physical and numerical simulations of hydrodynamics and heat exchange and their effect on the shape of the crystallization front as follows. CF curvature and local impurity distributions in the crystal depend largely on the melt hydrodynamics and the hydrodynamics-controlled local thickness of the heat and concentration boundary layers. The azimuthal melt inhomogeneity is predominant in the free convection mode. Starting from uniform crystal rotation and with the transition to the mixed convection allows controlling melt hydrodynamics and local heat and mass exchange at the CF by varying the relative contributions of free and forced convection [18][19][20]. The dynamic parameters of the melt are described by the Grashof (G r ), Marangoni (M a ) and Reynolds (Re c ) numbers, and the thermophysical melt properties, by the Prandtl number (P r ). G r describes the intensity of thermal gravity convection in the melt, M a describes the contribution of the thermocapillary effect and Re c depends on the intensity of the forced convection caused by the rotation of the growing crystal. Both numerical and physical simulations for axial symmetry problem statement suggest that there are some specific ratios of the thermodynamic parameters in the free convection mode, where the shape of the sub-crystal isotherm is the most planar and the radial distribution of the local heat flow in the boundary layer is the most homogeneous [21][22][23][24]. Thus rule holds for P r = 16.0, 45.6 and 0.05. The latter figure is typical for liquid metals and close to P r for silicon melt.
Real crystallization experiments with the model fluid being heptadecane having P r = 45.6 showed that at dynamic parameters of the melt flow close to those for the real single crystal growth conditions in the Czochralski technique, the CF shape under single phase convection conditions is almost perfectly flat. The experiment also showed that the axially symmetrical laminar 3D shape of the flow only exists until the threshold G r , M a and Re c numbers are reached. Once the G r , M a and Re c numbers exceeded the threshold level at the same ratio between parameters, the boundary layers under the rotating CF may lose stability, causing a system of azimuthal swirls to develop under the crystallization front. Figure 4 shows examples of 3D flow shapes after the development of such instability.
The photos and diagrams of flow patterns shown in Fig. 4 are as follows: Fig. 4a shows the side view of the flow in the central part of the crucible (in the height/radius plane, r-z). The transparent melt-simulating fluid was visualized with highly light reflecting tracer particles. The  visualized fluid layer is sectioned in the drawing plane by a flat light beam of about 1 mm in thickness. As a result, the flow only occurs in that thin layer. The flow pattern in Fig. 4b is shown by projections of fluid motion trajectories on the r-z plane. The trajectories are 3D and can be described as coils wound on irregular torus rings. The diagram shows the inner torus with an upward flow in the center, formed by centrifugal forces (forced convection region). The melt is attracted toward the rotating CF in this region. The outer torus is the free convection region. These two flows collide at the free melt surface or under the crystal: the free convective flow from the crucible walls and the centrifugal flow from under the crystal. Arrows show the melt flow direction in the meridian plane. The respective thin flow layer in the sub-crystal melt region is illuminated in a similar manner (Fig. 4c-e) but for the plane parallel to the transparent crucible bottom. Depending on the G r , M a and Re c ratios the flow collision boundary shape varies from elliptic (Fig. 4c, Curve 1) to polygonal (Fig. 4d).
With an increase in Re c (due to an increase in the angular velocity of rotation of the crystal ), the pentagon transforms to a decagon, the flow collision boundary shifting toward the crucible walls. The melt flows not only from the region under the crystal in radial directions; it also has an azimuthal velocity component. Secondary helical flows emerge against the background of radial melt spreading. In Fig. 4c these flows are directed along the major ellipse axis (Curve 2), and in Fig. 4d, toward the pentagon corners. Thus, heated melt flows within the boundary layer towards the CF with respective angular periods, and cooled melt flows from the crystallization front. Furthermore, at the inner torus periphery (Fig. 4a, b) the downward flow is no longer homogeneous in the azimuthal direction: more intense local streams emerge which are co-current with the main flow (they flow in the same direction but at a different velocity [25]). The downward portions of these streams are denoted with the + signs in the flow diagrams in Fig. 5c, d, the outer and inner circles being the crucible and crystal boundaries, respectively. Thus, in a specific G r , M a and Re c range depending on the relative melt height H/R c (R c being the crucible radius) and relative radius R c /R x (R x being the growing crystal radius), a system of well-developed radially oriented swirls may emerge under the rotating CF. In the absence of such swirls the melt is displaced from under the crystallization front in a homogeneous manner to form thermal and concentration boundary layers which are homogeneous in azimuthal direction but have clear radial inhomogeneity. Once swirls emerge the melt is displaced from the center to the periphery, and simultaneous fluid motion in azimuthal direction occurs. The resulting overall melt motion is helical. Heated melt impinges onto the crystallization front in upward flows and may locally melt the CF, while in downward flows the melt is cooled, so the CF may become locally convex toward the melt. The structures of the boundary concentration layers and the respective concentration flows are similar. This seems to be the case for convection modes which produce single crystal portions with azimuthal electrical resistivity inhomogeneity. The number of swirls for the preset Prandtl number depends on the Grashof, Marangoni and Reynolds numbers ratios and on the H/R c and R c /R x ratios under the rotating CF. Thus, considerable inhomogeneities with azimuthal symmetry components emerge in the temperature fields, thermal stresses and concentration fields.
Simulation of crystallization front shape dependence on local hydrodynamics also showed [21] that even for  a highly symmetrical temperature field a large-scale azimuthal CF instability may emerge, i.e., the front can be partially molten at one side (Fig. 5). This portion of the crystallization front remains partially molten as the crystal grows, this phenomenon being typical both for almost flat CF surfaces and for convexo-concave ones, causing the experimentally observed radial-azimuthal electrical resistivity distribution asymmetry. Figure 5 shows photos of the flow in the central part of the crucible for convection modes in which upward melt flows form under the uniformly rotating CF. In these modes free convection of thermal gravity capillary origin is almost completely suppressed. The flat crystallization front shape in this specific case exists at G r = 3575, M a = 3454 and Re c = 40.4, the 3D flow pattern for this mode being in a general agreement with the flow shape shown in Fig. 4a.

Conclusion
Analysis of the electrical resistivity distribution maps for Czochralski grown 150 mm silicon single crystals revealed two inhomogeneity types in the wafers: radial and azimuthal ones. The experimentally observed distribution features agree with the experimental results for transparent melt-simulating fluids. This suggests that the cause of the azimuthal inhomogeneities is the emergence of a swirl flow structure in the melt under the crystallization front, and the cause of the radial inhomogeneities is the curvature of the crystallization front arising due to local features of heat transfer. Comparison of numerical simulation results in a wide range of Prandtl numbers with the physical experimental data suggests that the phenomena observed in transparent fluids are universal and can be used for theoretical interpretation of imperfections in silicon single crystals.