Numerical Analysis of Phosphorus Concentration Distribution in a Silicon Crystal during Directional Solidiﬁcation Process

: For bulk doping, boron and phosphorus are usually used as p-type and n-type dopants, respectively. The distribution of these dopant concentrations in a silicon crystal along the vertical direction is governed by the segregation phenomena. As the segregation coefﬁcient of phosphorus is small, phosphorus concentration distribution in a silicon crystal becomes inhomogeneous; inhomogeneous phosphorus concentration distribution affects the distribution of resistivity in the crystal. Therefore, it is important to control the phosphorus concentration distribution in a silicon crystal and make it uniform. In this study, by numerical analysis, we investigated the effect of the evaporation ﬂux at the melt surface on phosphorus concentration distribution during the directional solidiﬁcation process. To obtain a homogeneous phosphorus concentration distribution in the silicon crystal, we had to control the phosphorous evaporation ﬂux at the melt surface and maintain approximately the same phosphorus concentration in the melt during the entire solidiﬁcation process even though the growth rate was always changing.


Introduction
The photovoltaic market has rapidly advanced over the last two decades. Compared with multicrystalline materials, single-crystalline silicon has proven advantages such as high quality and very low dislocation density, but the production cost of single-crystalline silicon is high. Therefore, multicrystalline materials are typically used in solar cells because they are more cost-effective than their single-crystalline silicon counterparts. The most widely used method for growing multicrystalline materials is directional solidification (DS). Moreover, technological progress regarding solar cells has been developed significantly within the last decade, leading to enhanced efficiency of solar cells and accounting for the increased use of phosphorus-doped n-type crystalline silicon materials [1,2].
To reduce production costs, the impurity distribution along the axis direction should be controlled and kept constant. However, ingots grown by most manufacturers have an inhomogeneous phosphorus concentration distribution along their height because of the small segregation coefficient of phosphorus. A small segregation coefficient leads to the low production yield and high production cost of silicon crystals because of the specific range of resistivity required in the crystals by industrial manufacturers. The resistivity of a p-doped CZ silicon crystal is widely distributed over the length of the crystal. The minority carrier lifetime of solar cells decreases as a function of the resistivity along the silicon crystal length [3]. Resistivity distribution is very effective regarding solar cell capability [4][5][6].
One method of improving the inhomogeneous resistivity distribution along the vertical direction could be simultaneous doping using two different impurities. Using numerical analysis, Lee conducted a study and reported that codoping Ga and Bi could improve the vertical distribution of Ga [7]. Using numerical simulation, Wang et al. conducted a study and presented that the codoping method is useful for improving silicon crystal growth [8].
Forster et al. proposed codoping with P, B, and Ga during ingot crystallization [9]. However, compensational doping may have negative effects such as the degradation of minority carrier lifetime and formation of recombination centers. Thus, it may not be suitable for high-efficiency solar cells.
Another method, the enhancement of dopant evaporation at the melt surface, is also effective in controlling the vertical resistivity distribution of grown crystals. Buchovska et al. have reported that traveling magnetic fields can enhance melt stirring and evaporation rate to control the phosphorus distribution along the vertical direction of the crystal [10]. In CZ silicon crystal growth, the growth rate, which is related to the incorporation flux of the dopant impurity at the melt-crystal interface, could be controlled accurately and kept constant during the growth process. However, the growth rate in the directional solidification (DS) method will always change during the entire solidification process; therefore, it is difficult to precisely control the growth rate.
In this study, we focus on the incorporation flux of phosphorus at the melt-crystal interface and evaporation flux of phosphorus at the melt surface. By numerical analysis, we investigate how to obtain a homogeneous phosphorus concentration distribution along the crystal's vertical direction via the DS method, even when the growth rate is not constant.

Materials and Methods
The DS furnace used in this study consisted of a crystal, silicon melt, crucibles, pedestals, heat shields, and heaters, as shown in Figure 1a. The crystal's height was 80 mm and its diameter was 108 mm. Figure 1b shows the computational grid of the meltcrystal domain. The melt-crystal interface is calculated using a dynamic interface tracking method. The melt-crystal interface is moving upward as a function of time. The grid size of boundary areas such as the melt-crystal interface, melt surface and melt-crucible wall interface is 0.007-0.1 mm. We developed our own simulation code as follows. A transient global model was developed for the DS process to study the global heat transfer in the entire furnace, flow field in the melt, and shape of the melt-crystal interface as a function of time. The global heat transfer in the entire furnace included the convective heat transfer of the silicon melt, conductive heat transfer in all of the solid components, and radiative heat transfer in all of the diffusive surfaces of the furnace. A dynamic interface tracking method was also included to obtain the shape of the melt-crystal interface. The following major assumptions were made for the simulations: (i) The geometry of the furnace configuration is axisymmetric; (ii) the melt flow is incompressible and laminar; (iii) the radiative heat transfer can be modeled as a diffuse-gray surface; and (iv) the effect of the gas flow in the furnace is negligible and can be ignored. The governing equations for melt flow under these conditions are as follows: where V, ρ, p, µ, g , β T , c and k represent melt velocity, melt density, melt pressure, melt viscosity, gravitational acceleration, thermal expansion coefficient, heat capacity and thermal conductivity, respectively. Radiative heat exchange between diffuse surfaces is important in heat transport in a unidirectional solidification furnace [11]. The calculation details, including the impurity transport in the melt, have already been reported elsewhere [12][13][14]. details, including the impurity transport in the melt, have already been reported elsewhere [12][13][14].
(a) (b) The governing equation of phosphorus transport in the silicon melt is where and are the phosphorus concentration and diffusion coefficient [5] in the silicon melt, respectively.
We considered the phosphorus evaporation flux at the melt surface to be represented as follows [15]: where is the diffusion coefficient in the silicon melt and is the phosphorus concentration at the melt surface. We applied the phosphorus incorporation flux at the meltcrystal interface, Fluxinterface, to the phosphorus evaporation flux at the melt surface, Fluxmelt surface. The phosphorus incorporation flux at the melt-crystal interface was represented as follows: where is the growth rate, is the phosphorus concentration at the melt-crystal interface, and is the segregation coefficient of phosphorus. We take = 0.35 [16]. The initial value of the phosphorus concentration in the melt was 10 16 (atoms/cm 3 ).  Figure 4b was almost homogeneous. In the other two cases, inhomogeneous distribution was noticed. Figure 5a shows the phosphorus concentration distribution in the center of the crystal along the vertical direction with different evaporation fluxes at the melt sur- The governing equation of phosphorus transport in the silicon melt is

Results and Discussion
where C and D are the phosphorus concentration and diffusion coefficient [5] in the silicon melt, respectively. We considered the phosphorus evaporation flux at the melt surface to be represented as follows [15]: where D is the diffusion coefficient in the silicon melt and C is the phosphorus concentration at the melt surface. We applied the phosphorus incorporation flux at the melt-crystal interface, Flux interface , to the phosphorus evaporation flux at the melt surface, Flux melt surface . The phosphorus incorporation flux at the melt-crystal interface was represented as follows: where V g is the growth rate, C i is the phosphorus concentration at the melt-crystal interface, and k is the segregation coefficient of phosphorus. We take k = 0.35 [16]. The initial value of the phosphorus concentration in the melt was 10 16 (atoms/cm 3 ).  The distribution of phosphorus concentration in the crystal for the case shown in Figure 4b was almost homogeneous. In the other two cases, inhomogeneous distribution was noticed. Figure 5a shows the phosphorus concentration distribution in the center of the crystal along the vertical direction with different evaporation fluxes at the melt surface. The distribution in case (i) was close to the normal freezing model because the phosphorus evaporation flux at the melt surface was zero. The distribution in case (iii) decreased as a function of the height from the bottom because phosphorus evaporation flux at the melt surface was larger than phosphorus incorporation flux at the melt-crystal interface. However, the distribution of phosphorus concentration in the center of the crystal along the vertical direction in case (ii) was almost homogeneous because the phosphorus evaporation flux at the melt surface was the same as the phosphorus incorporation flux at the melt-crystal interface. This meant that the incoming and outgoing phosphorus concentrations in the melt were the same even though the growth rate always changes during the solidification process. We also investigated antimony concentration distribution with different evaporation fluxes to confirm our simulation result even when the segregation coefficient is different. Figure 5b shows the antimony concentration distribution in the center of the crystal along the vertical direction with different evaporation fluxes at the melt surface. The segregation coefficient of antimony is 0.023 [16]. Tendency is similar as that of phosphorus.

Results and Discussion
face. The distribution in case (i) was close to the normal freezing model because the phosphorus evaporation flux at the melt surface was zero. The distribution in case (iii) decreased as a function of the height from the bottom because phosphorus evaporation flux at the melt surface was larger than phosphorus incorporation flux at the melt-crystal interface. However, the distribution of phosphorus concentration in the center of the crystal along the vertical direction in case (ii) was almost homogeneous because the phosphorus evaporation flux at the melt surface was the same as the phosphorus incorporation flux at the melt-crystal interface. This meant that the incoming and outgoing phosphorus concentrations in the melt were the same even though the growth rate always changes during the solidification process. We also investigated antimony concentration distribution with different evaporation fluxes to confirm our simulation result even when the segregation coefficient is different. Figure 5b shows the antimony concentration distribution in the center of the crystal along the vertical direction with different evaporation fluxes at the melt surface. The segregation coefficient of antimony is 0.023 [16]. Tendency is similar as that of phosphorus.  face. The distribution in case (i) was close to the normal freezing model because the phosphorus evaporation flux at the melt surface was zero. The distribution in case (iii) decreased as a function of the height from the bottom because phosphorus evaporation flux at the melt surface was larger than phosphorus incorporation flux at the melt-crystal interface. However, the distribution of phosphorus concentration in the center of the crystal along the vertical direction in case (ii) was almost homogeneous because the phosphorus evaporation flux at the melt surface was the same as the phosphorus incorporation flux at the melt-crystal interface. This meant that the incoming and outgoing phosphorus concentrations in the melt were the same even though the growth rate always changes during the solidification process. We also investigated antimony concentration distribution with different evaporation fluxes to confirm our simulation result even when the segregation coefficient is different. Figure 5b shows the antimony concentration distribution in the center of the crystal along the vertical direction with different evaporation fluxes at the melt surface. The segregation coefficient of antimony is 0.023 [16]. Tendency is similar as that of phosphorus.      Figure 6 shows the difference between the maximum and minimum phosphorus concentrations in the melt as a function of time; Figure 7 displays the history of the phosphorus evaporation flux at the melt surface as a function of time. At the beginning of the solidification time, the difference between the maximum and minimum phosphorus concentrations in the melt in case (i) was small because the phosphorus evaporation flux at the melt surface was zero. However, at the end of the solidification time, the difference between the maximum and minimum phosphorus concentrations in the melt in case (i)  Figure 6 shows the difference between the maximum and minimum phosphorus concentrations in the melt as a function of time; Figure 7 displays the history of the phosphorus evaporation flux at the melt surface as a function of time. At the beginning of the solidification time, the difference between the maximum and minimum phosphorus concentrations in the melt in case (i) was small because the phosphorus evaporation flux at the melt surface was zero. However, at the end of the solidification time, the difference between the maximum and minimum phosphorus concentrations in the melt in case (i) rapidly increased because of the effect of the small segregation coefficient of phosphorus. Therefore, the phosphorus concentration in the melt in case (i) substantially changes from the beginning to the end of the solidification time. Moreover, at the beginning of the solidification time, the difference between the maximum and minimum phosphorus concentrations in the melt in case (iii) was larger than that in the other two cases because of the large evaporation flux at the melt surface. However, at the end of the solidification time, the difference between the maximum and minimum phosphorus concentrations in the melt in case (iii) was smaller than that in the other two cases. The phosphorus evaporation flux at the melt surface was related to the growth rate and phosphorus concentration at the melt-crystal interface, as shown in Equation (6). In this case, the phosphorus concentration in the melt had a large effect on the behavior of the phosphorus evaporation flux at the melt surface in case (iii) because the growth rates in cases (ii) and (iii) were the same. The phosphorus concentration in the melt in case (iii) was lower at the end of the solidification time because the phosphorus had already evaporated considerably at the beginning. Therefore, the phosphorus concentration in the melt in case (iii) also changed substantially from the beginning to the end of the solidification time. In contrast, the difference between the maximum and minimum phosphorus concentrations in the melt in case (ii) remained almost constant throughout the solidification time. This caused the enhancement of the phosphorus evaporation at the beginning of the solidification time and almost kept the same phosphorus evaporation flux until the end of the solidification time, as shown in Figure 7. Therefore, we needed to maintain nearly the same phosphorus concentration in the melt during the whole solidification process to obtain a crystal with a homogeneous phosphorus concentration distribution along the vertical direction.     Figure  9 shows the distribution of the phosphorus concentration in the center of the melt along the vertical direction with different evaporation fluxes at the melt surface at one hour after the start of solidification. We could identify that there was a large difference in the phosphorus concentration at the melt surface because the phosphorous evaporation flux was very effective. However, the phosphorus concentration at the melt-crystal interface was   Figures 6 and 8a-c. The effect of the phosphorus evaporation flux at the melt surface on the phosphorus concentration at the melt-crystal interface was small because the melt volume was large, velocity of the melt convection was small, and solidification time was short. Figure 9 shows the distribution of the phosphorus concentration in the center of the melt along the vertical direction with different evaporation fluxes at the melt surface at one hour after the start of solidification. We could identify that there was a large difference in the phosphorus concentration at the melt surface because the phosphorous evaporation flux was very effective. However, the phosphorus concentration at the melt-crystal interface was close even though the evaporation flux at the melt surface was different because the phosphorous evaporation flux at the melt surface had little effect on the phosphorus concentration at the melt-crystal interface. Figure 8d-f show the phosphorus concentration distribution in the melt with different evaporation fluxes at the melt surface three hours after the start of solidification. In these three hours, the melt volume decreased and solidification time increased. Thus, the phosphorus concentration in the melt could have been more homogeneous, and the difference between the maximum and minimum phosphorus concentrations in the melt was small because the melt mixing was enhanced and the effect of the evaporation flux of phosphorus at the melt surface on the phosphorus concentration at the melt-crystal interface subsequent to three hours after the start of solidification was larger than that subsequent to one hour after the start of solidification. Figure 10 shows the distribution of the phosphorus concentration in the center of the melt along the vertical direction with different evaporation fluxes at the melt surface at three hours after the start of solidification. We can say that there was a considerable difference in the phosphorous concentration not only at the melt surface but also at the melt-crystal interface, and the phosphorous concentration in the melt subsequent to three hours after the start of solidification was more homogeneous than that subsequent to one hour after the start of solidification. Thus, the effect of phosphorus evaporation flux at the melt surface on the phosphorus concentration at the melt-crystal interface at the beginning of the solidification process is minor. of solidification. We can say that there was a considerable difference in the phosphorous concentration not only at the melt surface but also at the melt-crystal interface, and the phosphorous concentration in the melt subsequent to three hours after the start of solidification was more homogeneous than that subsequent to one hour after the start of solidification. Thus, the effect of phosphorus evaporation flux at the melt surface on the phosphorus concentration at the melt-crystal interface at the beginning of the solidification process is minor.  (c) (f) Figure 8. Distribution of the phosphorus concentration in the melt with different evaporation fluxes at the melt surface after one hour and three hours of solidification. The phosphorous evaporation flux at the melt surface was (i) 0 (a,d), (ii) the same as the phosphorus incorporation flux at the melt-crystal interface (b,e), and (iii) twice as large as that of case (ii) (c,f).

Conclusions
We investigated the relationship between the phosphorus concentration distribution in a silicon crystal and phosphorus evaporation flux at the melt surface during the DS process by numerical analysis. We focused on the incorporation flux of phosphorus at the melt-crystal interface and phosphorus evaporation flux at the melt surface. To realize a homogeneous phosphorus concentration distribution in the silicon crystal, we had to control the evaporation flux of phosphorous at the melt surface and maintain approximately the same phosphorus concentration in the melt during the entire solidification process and we had to consider the incoming and outgoing phosphorus concentrations in the melt, such as the phosphorus incorporation flux at the melt-crystal interface and evaporation flux of phosphorous at the melt surface during the solidification process even though the growth rate was always changing.

Conclusions
We investigated the relationship between the phosphorus concentration distribution in a silicon crystal and phosphorus evaporation flux at the melt surface during the DS process by numerical analysis. We focused on the incorporation flux of phosphorus at the melt-crystal interface and phosphorus evaporation flux at the melt surface. To realize a homogeneous phosphorus concentration distribution in the silicon crystal, we had to control the evaporation flux of phosphorous at the melt surface and maintain approximately the same phosphorus concentration in the melt during the entire solidification process and we had to consider the incoming and outgoing phosphorus concentrations in the melt, such as the phosphorus incorporation flux at the melt-crystal interface and evaporation flux of phosphorous at the melt surface during the solidification process even though the growth rate was always changing.