- Email: [email protected]
Micro structures in the grain evolution during solidification of silicon: Phase field calculations
Acta Materialia 140 (2017) 1e9
Contents lists available at ScienceDirect
Acta Materialia journal homepage: www.elsevier.com/locate/actamat
Full length article
Micro structures in the grain evolution during solidification of silicon: Phase field calculations W. Miller a, *, A. Popescu b a b
Leibniz Institute for Crystal Growth (IKZ), Max-Born-Str.2, 12489 Berlin, Germany Faculty of Physics, West University of Timisoara, Bd. V. Parvan 4, 300223 Timisoara, Romania
a r t i c l e i n f o
a b s t r a c t
Article history: Received 3 March 2017 Received in revised form 8 August 2017 Accepted 10 August 2017 Available online 12 August 2017
Two dimensional phase field simulations have been performed to study the influence of the growth kinetics and the surface energy on the growth behaviour of grains during solidification of Si. In particular, we studied the groove between two grains as recently observed by in situ observation [1]. Furthermore, we performed computations for different S boundaries and discuss the interplay between equilibrium of interfacial energies and growth kinetics. © 2017 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved.
Keywords: A1. Computer simulation A1. Directional solidification A1. Growth models B2 Semiconducting silicon
1. Introduction Despite the great importance of direct solidification of silicon and the large number of investigations there are still many open questions concerning the grain evolution. Post-mortem analysis of the grain structure is one way for a better understanding of the growth kinetics. Most often these analyses are performed on wafers giving some statistics at a certain height but not the continuous evolution in time. Analysis in vertical direction has been performed e.g. by Prakash et al. [2] and more recently by Lin et al. [3]. However, only the situation after solidification and cooling can be analysed. A better path towards understanding the details is in-situ investigation which became recently available by improved synchrotron X-ray imaging to visualize the melt/solid interface of silicon [4]. Thus, it was possible to follow the creation and disappearance of facetted/facetted grooves [1]. For an understanding of the basic mechanisms a simulation on the microscopic length scale is required, which includes the effect of interface curvature. Here, phase field methods are the preferential ones. Chen et al. studied the grain competition in an undercooled melt with a simple anisotropy function for the solid/ melt interfacial energy. Cantù et al. and Miller et al. studied a more
* Corresponding author. E-mail address: [email protected] (W. Miller). http://dx.doi.org/10.1016/j.actamat.2017.08.025 1359-6454/© 2017 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved.
advanced anisotropy function but simplified growth kinetics [5,6]. All aforentioned calculations have been two-dimensional approaches. Lin & Lan used a three-dimensional model, still with simplified anisotropy function and growth kinetics [7]. So far, all phase field calculations were based on the model of Warren et al. for computing the grain growth [8,9]. Also in this paper we use such a model and discuss in the end the limitations of this approach concerning the grain boundary orientation. Compared to our previous investigations [5,6] we update the solid/melt interfacial energy according to recent numerical and experimental results [10,11]. We also set the growth kinetics based on the results of molecular dynamics calculations [12]. The much more realistic input parameters enabled us to compare our calculations with the results of in situ investigations by Tandjaoui et al. [1]. 2. Interfacial energies In this paper we consider the (110) plane of silicon for 2D calculations (see Fig. 1). We need two kinds of interfacial energies for our computations: that of the melt/crystal interface and those of the grain/grain boundaries. In particular, we need the melt/crystal interface energy as a function of the orientation. Recently, the results of a molecular dynamics study [10] and of an experiment with in-situ observation [11] have been published. In both studies the equilibrium shape in the (110) plane was the result (see Fig. 2). The equilibrium shape represents the relative change of the interface
2
W. Miller, A. Popescu / Acta Materialia 140 (2017) 1e9
Fig. 1. Plane for calculations.
Fig. 2. Equilibrium shape for interface energy of Eq. (1) (blue line) in comparison with result of MD simulation (left [10]) and experiment (right [11]). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
energy with respect to orientation. Absolute value of interface energies were computed earlier for three orientations by means of molecular dynamics method [13]. Apte employed the cleaving wall method to extract the interfacial energies as follows: 0:42±0:02 J=m2 , 0:34±0:02 J=m2 , 0:35±0:03 J=m2 for the {100}, {111}, and {110} planes, respectively. In the computations we need to represent the interface energy by a function, which is differentiable twice with respect to orientation. In order to recover an equilibrium shape close to the observed ones we made the following. appoximation:
.
n
g ¼ 0:335 J m2 1: þ 0:16ð1 tanhð20ðw 50 ÞÞ þ 0:03ð1
for {100} and {110} planes and two deep for {111} and {449} planes. Grain boundary energies for silicon were computed by Kohyama and co-workers using a bond counting method [15]. More recently, molecular dynamics have been employed to obtain the energy at certain S boundaries [16]. In Fig. 4 we show the data points of the two kinds of calculations in the ð111Þ plane. The original values of Kohyama et al. have been scaled by 0.65 to adapt them to the values of the molecular dynamics calculations. Thus, they are consistent with the absolute values of the melt/solid interfacial energy, which are also based on the results of molecular dynamics calculations. For small misorientations one can apply the Read-Shockley theory, Dq 1 ln Dq where gGB ¼ g0GB Dq , where Dq and Dqm are the Dq m
tanhð18ð60 wÞÞ 0:02 sin6 w 1 0:05 cos8 ðwÞ o 0:04 cos128 ðw 34 Þ : (1) The equilibrium shape given by this equation is shown in Fig. 2 as a blue line. Besides the equilibrium shape the three absolute values for the interface energy as given above were a criterium for setting up Eq. (1). We almost match the values given by molecular dynamics calculations as Eq. (1) yields 0:425 J=m2 ({100}), 0:337 J=m2 ({111}), and 0:348 J=m2 ({110}). The interface energy as a function of the orientation according to Eq. (1) is shown in Fig. 3. The curve is similar to the curve derived from bond density calculations [14]. We have four minima: two flat
m
missorientation and the saturation, respectively. In Fig. 4 this energy is shown for g0GB ¼ 0:68 J=m2 and Dqm ¼ 28o as a light grey curve. For the computations we use the values of the blue curve. This is discussed later in the context of the phase field model we use for computing the phase transition.
3. Growth kinetics The growth velocity of a melt/solid interface is given by
vn ¼ bðDTÞDT;
(2)
where DT is the undercooling at the considered point of the interface. The kinetic coefficient b depends on the growth
W. Miller, A. Popescu / Acta Materialia 140 (2017) 1e9
3
Fig. 5. Kinetic coefficient as a function of the missorientation from the 〈111〉 direction. Red diamonds are the data obtained by molecular dynamics calculations [12]. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 3. Melt crystal interfacial energy as given by Eq. (1).
Below 1:2o either Eq. (3) or Eq. (4) for dislocation or 2D nucleation is applied. 4. Numerical methods 4.1. Phase transition The phase transition is treated by the phase-field method as explained in Ref. [5]. It is based on the model of Warren et al. for computing the phase transition and evolution of grains [8,9]. Here, we shortly repeat the main equations. The evolution equations for the order parameter of solid (f ¼ 1) and liquid phase (f ¼ 0) and for the grain orientation (p=2 < j p=2) are:
tf vt f ¼ x2f V2 f þ 4Wfð1 fÞðf 1=2Þ DHV
þ Fig. 4. Grain boundary energy as a function of the misorientation. The red routes are the values obtained by Kohyama with the bond-counting method [17]. The original energies were scaled by 0.65. The two triangles are from molecular dynamics calculations and are taken as a reference for the data of Kohyama [16]. The blue curve represents the curve given by Eq. (13). The light grey curve is the grain boundary energy according to the Read-Shockley theory. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
mechanism. For growth on atomically rough surfaces b does not depend on DT. For Si(100) it was computed by means of molecular dynamics as bð100Þ ¼ 0:12 m=ðs KÞ [12]. At the 〈111〉 facet either a dislocation induced growth can occur, where the kinetic coefficient obeys:
bðDTÞ ¼ bdisl DT;
(3)
or a 2D nucleation, where b is given by
bðDTÞ ¼ aA expfaB =DTgDT 1=3 m=ðs KÞ:
2 Tm
For the {111} plane experimental and numerical attempts have been made to obtain the parameters aA and aB [18e20]. Once a nucleus exists a step growth starts. This type of growth has been investigated by molecular dynamics calculations [12]. The data can be interpolated by lines in two regimes as shown in Fig. 5.
(5)
2sfjVjj x2j fjVjj2 s þ x2j Vj ; f2 tj Vt j ¼ V f2 jVjj
(6)
where xf and xj are the phase field parameters related to the liquid/ solid and grain/grain energy, respectively. DHV and Tm are the latent heat per volume and the melting-point temperature, respectively. W is the barrier for the solid/liquid phase transition at T ¼ Tm and s is a coupling constant. tf is defined by the kinetic coefficient, which was discussed in the previous section. tj is the relaxation time for the grain-grain phase transition and is chosen ten times smaller than tf . Eqs. (5) and (6) are solved by finite difference and a kinetic method, respectively [5]. The solid-liquid interface energy is given by
pffiffiffiffiffiffiffiffi
1
(4)
30f2 ð1 fÞ2 ðT Tm Þ
gsl ¼ pffiffiffiTm xf 2W ; 3 2
(7)
and the grain-grain boundary energy is related to this via
gGB ¼ 2ssl 1 f3min :
(8)
4
W. Miller, A. Popescu / Acta Materialia 140 (2017) 1e9
fmin is the minimum value of f between two grains. From a 1D analysis for T≡Tm one obtains [8]:
fmin ¼ 1
Dq ; Dqc
(9)
where the critical difference in crystallographic orientations is
Dqc ¼
pffiffiffiffiffiffiffiffi 2W xf 1 ¼ : facs facx s
(10)
facx ¼ xj =xf is a model parameter. We chose facx ¼ 8=15 as it was pffiffiffiffiffiffiffiffi done by Warren et al. [8]. facs ¼ s=ð 2W xj Þ defines the stiffness of the grain-grain boundary. For a fixed facs the grain boundary energy as given by Eq. (8) obeys a similar behaviour as those of the Read-Shockley model:
gGB ¼ E0
Dq Dq 1 ln Dqm Dqm
:
(11)
In the Read-Shockley model the grain boundary energy saturates at E0 for large misfits (larger than the material parameter Dqm ). However, gGB exhibits an oscillating behaviour as discussed in section 1 (see Fig. 4). We adapt to this behaviour by varying facs as a function of the misfit angle. In particular, we divide the misfit into two regimes: up to a misfit of Dq ¼ 20o the grain boundary energy can be described by Eq. (11). Beyond we use the values given in Ref. [17] to define the grain boundary energy. In particular the grain boundary energy is given by:
! . Dq 3 2 gGB ðDqÞ ¼ 0:67 J m 1 1 ;
Dqc
(12)
where Dqc is given by:
Dqc ¼
8 > < > : 20o þ
Dq
2
4:47o
30o þ 0:5Dq sinð28Dq 280o Þ
if Dq 20o if Dq > 20o (13)
We have performed computations for a 2D domain with a uniform grid (see Fig. 6). At the bottom we have two grains with a solid/melt interface at x ¼ xint , where x ¼ 0 is at the top and x ¼ xend at the bottom. The width of the diffuse interface (0:05 f 0:95) is set ot be 10 grid points. 4.2. Computational domain and heat transport
R cool Gmelt T
;
transport equation
vT vf DHV ¼ VðDT VÞT þ vt vt rcp
(15)
is solved by a kinetic scheme [26]. In the diffuse interface the ð1 fÞ þ Dsolid f. thermal diffusivity is given as DT ¼ Dmelt T T As sketched in Fig. 6 we start from an artificially initial configuration of grain boundary and solid/melt interface shape. In the end we like to obtain a steady state situation, where the interface is moving without changing its shape and the grain boundary keeps its orientation. In order to reach this state a large number of time steps is required and the solid/melt interface is moving a significantly distance in the meantime. In order to keep the configuration as close to the original one and allow a large number of time steps we perform a rezoning once the average interface moved by one grid spacing [5]. In the solid part the last line of grid points is leaving the computational domain while in the melt part a new line enters. All variables are set to the same value at the adjacent line of grid points. When rezoning the temperature at the boundary x ¼ 0 is reset to its original value Tbc ðt ¼ 0Þ. Only in the steady-state situation the interface is moving approximately with the velocity vn based on the given cooling rate. In the mean time the interface velocity might be quite different because of effects of growth kinetics or changes in the grain boundary. Thus, solving everything in a moving reference frame is not useful. The entire code was parallelized by using the PETSC library version 3.1 [27].
5. Results and discussions
We try to mimic the typical process in a solidification experiment as known from both large scale production experiments (see e.g. Ref. [25]) or small scale research experiments (see e.g. Ref. [4]). The procedure of setting the boundary condition is described in detail in Ref. [5]. The computational domain with boundary conditions is sketched in Fig. 6. The heat flux at the bottom boundary is given by:
j ¼ Dmelt rcp Gmelt þ DHV T T
Fig. 6. Domain for computations with thermal boundary conditions.
(14)
The region around the three phase junction of grains and melt is a complex system, where the interfacial energy effects are coupled with the local temperature field and the growth kinetics. Before discussing the numerical results we discuss the situation at the three phase junction with respect to the interfacial energies without considering the growth kinetics. At the three phase junction between two grains and the melt the interfacial forces between the grains, grain 1 and melt, and grain 2 and melt should be balanced in equilibrium. For the anisotropic system this balance is given by the following equation [28]:
where Dmelt , r, and cp are the thermal diffusivity, the density, and T the specific heat. The cooling rate R cool and the temperature
3 X
are the two process parameters. In exgradient in the melt Gmelt T periments, the latter is typically known only indirectly by
j¼1
¼ R cool =vn . The heat measuring the growth velocity vn : Gmelt T
!
gj b j þ
vgj ! nj ¼ 0 vw
(16)
For a groove the situation is schematically drawn in Fig. 7. In 2D we end up with:
W. Miller, A. Popescu / Acta Materialia 140 (2017) 1e9
Fig. 7. Angles and orientations in Young-Herring's force balance.
vgGB vgsl þ gsl ðq1 Þsinðq1 Þ þ cosðq1 Þ vw vw q1 vgsl cosðq2 Þ ¼ f1 gsl ðq2 Þsinðq2 Þ vw q2
(17)
vgsl sinðq1 Þ vw q1 vgsl sinðq2 Þ ¼ f2 : þgsl ðq2 Þcosðq2 Þ vw q2
(18)
gGB þ gsl ðq1 Þcosðq1 Þ
In Eqs. (17) and (18) we use as a more general version with f1 and f2 on the right hand side. In the equilibrium where all forces cancel f1 ≡f2 ¼ 0. In all phase field computations we have gGB independent vgGB vw ≡0.
The equations of the orientation of the grain boundary, i.e. can be used in different ways. On one hand we can compute the angles q1 and q2 for given grain boundary energy. This enables us to countercheck the results of the numerical calculations and see if they obey the equilibrium or if there is an influence of growth kinetics. The equations can be used also in the opposite way. Knowing q1 and q2 e.g. from experimental results the corresponding grain boundary energy can be computed (see e.g. Ref. [29]). The latter we did for the computation of a facetted/facetted groove as observed during in-situ observation of solidification. 5.1. Computation for facetted/facetted interface Firstly, we consider the case of facetted interfaces at both sides
5
of the groove. We refer to a case obtained by Tandjoui et al. in an experiment with in-situ observation [1]. They used synchrotron Xray radiation to visualize the solid/melt interface exploiting the differences in density of molten and solid silicon [4]. The solidification took place in a pyrolithic boron nitride crucible with a diameter of 6 mm applying a cooling rate of R cool ¼ 0:5 K=min. The growth velocity could be directly obtained from the in-situ measurement yielding vg ¼ 9:2 mm=s. Assuming that the cooling rate at the interface is the same, the temperature gradient can be computed as GT ¼ R cool vg , which yields GT ¼ 900 K=m. In particular, we consider the groove as shown on the right hand side of Fig. 8 and apply our numerical approach to compute the interface shape and grain boundary orientation. As no EBSD images are available for the crystal of this experiment we do not know the orientation of the grains. However, we know that the facets are {111} planes. We assume that the imaging plane in the experiment was a {110} plane to be consistent with our 2D model. According to Fig. 1 there are always two possibilities for the nearest {100} plane. For the facet of left grain (grain 1) the 〈100〉 orientation is either 50o or 20o off the main growth direction. For the other facet the possibilities are 10o or 80o. Thus, in total four combinations are possible. In Fig. 8 the configuration for 20o and 80o is shown. In case that the 〈100〉 direction of grain 1 is 50 off with respect to the main growth direction there is a {111} facet just 4.7 off the flat interface as shown in Fig. 9. This facet will be formed, which results in an undercooling at this facet and also at the left side of the groove. Thus, there is no chance to establish a long facet at the groove and a deep groove as observed in the experiment. Therefore, we computed only cases with 20o off for the left grain. In the following we denote the case with 20o and 10o (Dq ¼ 10o ) as case B and the case with 20o and 80o (Dq ¼ 60o ) as case A. In the experiment a depth of dgroove ¼ 240 mm was observed (see Fig. 8). Using this value the undercooling at the bottom can be ¼ 0:22 K. This suggest a dislocation estimated as DT ¼ dgroove Gmelt T growth mechanism rather than 2D nucleation. Therefore, we use
bdisl ¼ 1 105 m s1 K2 . Beyond an angle of 1:2o off the 〈111〉 direction we apply the step growth rules as shown in Fig. 5. This holds up to 25o , when it reaches the value for rough surfaces (b ¼ 0:12 m s1 K1 ). We performed runs for case A and case B for a domain of 250 250 mm2 with a lattice spacing of Dx ¼ 0:5 mm (parameter see Table 1). Thus, the width of the diffuse interface was wf ¼ 5:0 mm. As we know the geometry of the groove from experiment we can compute the grain boundary energy by Eq. (18). In fact, we used Eqs. (17) and (18) to compute q1 and q2 as a function of the grain boundary energy gGB. We start with case A: the target values for both q1 and q2 are 30o . For a grain boundary energy of gGB ¼ 0:59 J=m2 the computation yields q1 ¼ 28:9o and q2 ¼ 30:6o . According to Fig. 5 for these orientations we have a dislocation growth mechanism. The scaled value of Kohyama for the misfit of
Fig. 8. Left: Experimental observation by Tandjaoui et al. [1] (part of Fig. 2). Right: Orientations in the facetted/facetted groove.
6
W. Miller, A. Popescu / Acta Materialia 140 (2017) 1e9
Fig. 9. Result of computation for the case that the 〈100〉 direction of grain 1 is 50 off with respect to the main growth direction. A facet is formed at the interface, causing an undercooling in this region, which prevents the formation of a deep groove.
Table 1 Material parameters of Si used in simulations. thermal diffusivity in melt
Dmelt T
2:6 105 m2 =s
thermal diffusivity in solid
Dsolid T DHV rcp
9 106 m2 =s
[22]
4:22 109 J=m3
[23]
2:2 108 J=ðK m3 Þ 1683 K
[24]
latent heat specific heat melting-point temperature Gibbs-Thomson coefficient
Tm
G
[21]
1:8 107 K m
Dq ¼ 60o is gGB ¼ 0:65 J=m2 and thus it is slightly larger than the
computed one. In the phase field computations we used a fixed value of Dqcrit ¼ 116o representing gGB ¼ 0:59 J=m2 . As the initial configuration for the run we use the geometry known from the experiment, i.e. a groove with 〈111〉 facets on both sides a grain boundary with a 15 tilt with respect to the main growth direction. After consolidation of the temperature field and establishing the curved interface at the upper boundaries of the groove a steady state growth of the crystal was observed. The depth of the groove in this steady motion was dgroove ¼ 28 mm. The interface shape at two different time steps is shown in Fig. 10 demonstrating that there is no further change in the interface shape. The temperature profile along a line crossing the bottom in the groove is shown in the right part of Fig. 11 giving an undercooling of DT ¼ 0:112 K at the bottom of the groove. This is much larger than the estimation using Gmelt T and the depth of the groove which yields DT ¼ 0:025 K. We will come later back to this point. Here we like to note that the depth in the calculation is much smaller than that observed in the experiment (dgroove ¼ 240 mm). The reason for the discrepancy is the limitation of undercooling at the bottom of the groove in the numerical calculation due to the spatial resolution. At the bottom we have curved interface of curvature k. At that point the melting point temperature Tm is reduced by DTcurv ¼ 2Gk, where G ¼ gDHV =Tm is the Gibbs-Thomson coefficient. For a facetted/facetted groove k should be very large. However, it is limited by the thickness of the diffuse interface: kmax z2=wf . In the case of Dx ¼ 0:5 mm we have wf ¼ 5 mm and thus DTz0:15 K. This corresponds approximately with the undercooling in the computational run. We increased the resolution to Dx ¼ 0:2 mm taking a domain with 2000 2000 grid points. Now we obtained a depth of dgroove ¼ 103 mm with an undercooling of DT ¼ 0:285 K. Still the depth is smaller than the
Fig. 10. Interface lines in Case A for two time steps. The system is in a steady-state mode. The angles of the facets and the grain boundary are in the range of the experimentally observed ones (see top right sketch). The colors of the lines indicate the kinetic coefficient: red and blue means large and small value, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
one observed in the experiment. For a more realistic computation a local grid refinement is required. In order to estimate the undercooling at the bottom of the groove in the experiment we extrapolated the computed undercoolings for two different depth of the groove to the on measured in the experiment. Here, we took into account the limit for a vanishing groove: dgroove ¼ 00DT ¼ 0. We obtained DTz0:6 K. This undercooling is much smaller than required for 2D nucleation. It is much more likely that the growth on the facets was induced by a dislocation. Using bdislocation ¼ 2:8 104 m s1 K1 given by Voronkov [30] the necessary undercooling for vn ¼ 9:2 mm=s is 0.03 K. This is appoximately the undercooling at top of the facets. Thus, we expect that the growth on the facets was induced by at least one dislocation on each facet. The undercooling has also an impact on the probability of twin nucleation on the {111} facets. In the experiments it was observed that the facetted/facetted grooves most often disappear by the nucleation of a twin. According to the estimations by Duffar and Nadri [31] a reasonable probability for twinning is reached for undercoolings larger than 1.5 K, which is obviously much larger than observed in experiments. Recently, Lan and co-workers developed an advanced model for a twin seed, yielding a probability of about 1 109 already for an undercooling of DT ¼ 0:5 K. At DT ¼ 0:6 K it is already about 1 107 [32]. Therefore, the occurrence of twinning in the experiments is explainable by this model in conjunction with the temperature computation at the groove. In case B the misorientation is small ðDq ¼ 10o Þ and in the ReadShockley regime. From Fig. 4 one can see that gGB < 0:45 J=m2 . Performing the computations with Eqs. (17) and (18) we obtain gGB ¼ 0:58 J=m2 for angles q1 ≡q2 ¼ 28:8o . This means that the computed value using the geometry and our approach for the melt/ solid interface energy is much larger than the one from our model for the grain boundary energy. Either our models for the interface energies are not appropriate or this configuration was not realized in the experiment. Nevertheless, we performed runs for this case using a fixed value of Dqcrit ¼ 20o , representing gGB ¼ 0:58 J=m2 for Dq ¼ 10o . At the end we like to make some comments on the growth kinetics. Looking in more detail on the orientation of the two facets one observes that the angle with respect to the 〈111〉 direction is
W. Miller, A. Popescu / Acta Materialia 140 (2017) 1e9
7
Fig. 11. Temperature profile along line through deepest point of groove. Left: Line in the configuration of run. The light blue curve indicates the isoline of the melting point temperature. Right: temperature profile in red and phase field profile in blue. Regions a,b, and c mark the regions of groove, interface, and grain boundary, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
about 1:2o , i.e. near the boundary where the dislocation growth changes to step growth. There is step flow growth all along the facets except at the top end and at the bottom. 5.2. S boundaries As we mentioned already in all our calculations we use a unique value for the grain boundary energy for a given misfit. Therefore, the grain boundary energy does not depend on the orientation of the grain boundary. Of course, this is a simplification of the reality. Especially for coincidence lattice boundaries (S boundaries) one orientation is energetically favoured with respect to all others. This energy minimizing is valid for the bulk. The question is what happens at the three phase junction, where the orientation of the grain boundary to fulfill the force balance might be completely different from the lattice coincidence orientation (see Fig. 12). The most prominent S boundary is the S 3ð111Þ twin boundary. All computations agree that its boundary energy is negligible. Thus, it is most likely that two grains in this twin relation always form a (111) plane at the interface regardless of the local situation at the interface. It is also expected that the interface is almost flat with
only a small groove at the three-phase junction. However, the question is if this holds for all orientations of the growth direction with respect to the S 3ð111Þ twin boundary. We used Eq. (18) to obtain for a given q1 the corresponding angle q2 . In Fig. 13 the sketch of angles and the angle of the groove as a function of the angle between grain boundary and growth direction is shown. It is 4i ¼ 90o qi (i ¼ 1; 2) and 42 is negative in the figure. The growth direction is assumed to be symmetric with respect to the groove and in such a way wgrain is defined. A value wgroove > 180o means not a groove but a hillock, which will be unstable during growth. In Fig. 13 only the results are shown where the resulting interfaces at grain 1 and grain 2 are both non-missing interfaces. If one of the interfaces is a plane of a missing direction according to Eq. (1) (see Fig. 12) the situation is rather unstable and the values for wgroove are fluctuating. In Fig. 13 we show all data points used and it can be clearly seen that there are certain ranges of wgrain where there are no points, i.e. at least one of the direction of the planes is missing one. Next, we consider the S19 boundary, which has a grain boundary energy of 0:52 J=m2 (at the coincidence plane (133)). We performed runs for two different growth directions as shown in Fig. 14. In Case 1 there is only a small deviation between the computed grain boundary (white line) and the (133) coincidence plane (dotted orange line). Both sides of the groove are 〈111〉 facets. In Case 2 the deviation is larger. The groove is not as deep as in Case 1, because on both sides the interface is slightly off the (111) plane exhibiting a step growth mode with a larger kinetic coefficient b. Please note that in our phase field model there is no variation of the grain boundary energy with orientation and thus, there is no energy minimization with respect to this energy. In the future we will adapt the multi phase model for our purposes, which allows the definition of an anisotropy in the grain boundary energy. 6. Conclusions
Fig. 12. Stiffness of melt/solid interface. There are three regions of missing directions. The region around w ¼ 20o is enlarged for better visibility.
We applied a phase field model to compute the grain evolution during the solidification of silicon on the microscopic scale. We could recover the experimental observation by Tandjaoui et al. [1] of a groove with two 〈111〉 facets. The main issue is a very fine resolution at the bottom of the groove to reflect a large curvature at that point sustaining the undercooling caused by growth kinetics. The undercooling at the bottom is about 0:6 K. According to the recent work of Lin et al. [32] this undercooling is enough for a
8
W. Miller, A. Popescu / Acta Materialia 140 (2017) 1e9
Fig. 13. Computation of the angle of the groove wgroove as a function of the grain boundary with respect to the growth direction. For a given f1 we compute f2 according to Eq. (18). gGB Using Eq. (17) one can compute the required vvw .
Fig. 14. S19 grain boundary for two orientations. The dotted orange line indicates the S19 boundary, i.e. the (133) coincidence plane. The yellow line represents the isoline of Tm . In Case 1 the computed grain boundary is very close to the (133) coincidence plane. In Case 2 the deviation is much larger. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
reasonable probability of twin nucleation. This corresponds with the in-situ observations in the experiments [1]. We also investigated cases with S boundaries. Here, one would expect that the grain boundary is the coincidence plane of the two grains. However, the situation at the three phase junction with respect to interfacial energies and in addition growth kinetics might prefer a different orientation of the grain boundary. Part of the story is the weakness of the currently used phase field model, which does not allow to include an orientation dependent grain boundary energy. Therefore, we will change to a multi phase field model in the future. Acknowledgments
[5] [6] [7]
[8]
[9] [10]
[11]
W.M. likes to acknowledge fruitful discussions with Nathalie €el and her research group. Mangelinck-No
[12]
References
[13]
[1] A. Tandjaoui, N. Mangelinck-Noel, G. Reinhart, B. Billia, T. Lafford, J. Baruchel, Investigation of grain boundary grooves at the solideliquid interface during directional solidification of multi-crystalline silicon: in situ characterization by X-ray imaging, J. Cryst. Growth 377 (2013) 203e211. [2] R.R. Prakash, K. Jiptner, J. Chen, Y. Miyamura, H. Harada, T. Sekiguchi, Grain boundary interactions in multicrystalline silicon grown from small randomly oriented seeds, Appl. Phys. Express 8 (3) (2015) 035502. [3] H.K. Lin, M.C. Wu, C.C. Chen, C.W. Lan, Evolution of grain structures during directional solidification of silicon wafers, J. Cryst. Growth 439 (2016) 40e46. €l, G. Reinhart, J.-J. Furter, B. Billia, T. Lafford, [4] A. Tandjaoui, N. Mangelinck-Noe J. Baruchel, X. Guichard, Real time observation of the directional solidification
[14] [15]
[16]
[17] [18]
of multicrystalline silicon: X-ray imaging characterization, Energy Procedia 27 (2012) 82e87. G. Cantù, A. Popescu, W. Miller, Grain growth of silicon, Acta Mater 60 (2012) 6755e6761, http://dx.doi.org/10.1016/j.actamat.2012.08.048. W. Miller, A. Popescu, G. Cantù, Solidification of multicrystalline silicon simulation of micro structures, J. Cryst. Growth 385 (2014) 127e133. H.K. Lin, C.W. Lan, Three-dimensional phase field modeling of silicon thin-film growth during directional solidification: facet formation and grain competition, J. Cryst. Growth 401 (2014) 740e747. J.A. Warren, R. Kobayashi, A.E. Lobkovsky, W.C. Carter, Extending phase field models of solidification to polycrystalline materials, Acta Mater 51 (2003) 6035. J.A. Warren, R. Kobayashi, W.C. Carter, Modeling grain boundaries using a phase-field technique, J. Cryst. Growth 211 (2000) 18e20. G. Barinovs, A. Sabanskis, A. Muiznieks, Study of silicon crystal surface formation based on molecular dynamics simulation results, J. Cryst. Growth 391 (2014) 13e17. X. Yang, K. Fujiwara, K. Maeda, J. Nozawa, H. Koizumi, S. Uda, Crystal growth and equilibrium crystal shapes of silicon in the melt, Prog. Photovolt. Res. Appl. 22 (2014) 574e580. D. Buta, M. Asta, J.J. Hoyt, Kinetic coefficient of steps at the Si(111) crystalmelt interface from molecular dynamics simulations, J. Chem. Phys. 127 (07) (2007) 074703. P.A. Apte, X.C. Zeng, Anisotropy of crystal-melt interfacial free energy of silicon by simulation, Appl. Phys. Lett. 92 (22) (2008) 221903. P.J. Hesketh, C. Ju, S. Gowda, E. Zanoria, S. Danyuk, Surface free energy model of silicon anisotropic etching, J. Electrochem. Soc. 140 (4) (1993) 1080e1085. M. Kohyama, R. Yamamoto, M. Doyama, Reconstructed structures of symmetrical 〈011〉 tilt grain boundaries in silicon, Phys. Stat. Sol. B 138 (1986) 387e397. X. Chen, L. Xiong, A. Chernatynskiy, Y. Chen, A molecular dynamics study of tilt grain boundary resistance to slip and heat transfer in nanocrystalline silicon, J. Appl. Phys. 116 (2014) 244309. M. Kohyama, Computational studies of grain boundaries in covalent materials, Modelling Simul, Mater. Sci. Eng. 10 (2002) R31eR59. K.A. Jackson, Response to: some remarks on the undercooling of the Si(111)
W. Miller, A. Popescu / Acta Materialia 140 (2017) 1e9
[19]
[20] [21]
[22] [23] [24] [25]
facet and the “Monte Carlo modeling of silicon crystal growth” by Kirk M. Beatty & Kenneth A. Jackson, J. Cryst. Growth 211 (2000), 13 by W. Miller, J. Crystal Growth 325 (1) (2011) 104. W. Miller, Some remarks on the undercooling of the Si(111) facet and the “Monte Carlo modeling of silicon crystal growth” by Kirk M. Beatty & Kenneth A. Jackson, J. Cryst. Growth 211 (2000) 13. J. Crystal Growth 325 (1) (2011) 101e103. K.M. Beatty, K.A. Jackson, Monte Carlo modeling of silicon crystal growth, J. Cryst. Growth 211 (2000) 13e17. Y. Inatomi, F. Onishi, K. Nagashio, K. Kuribayashi, Density and thermal conductivity measurements for silicon melt by electromagnetic levitation under a static magnetic field, Int. J. Thermophys. 28 (1) (2007) 44e59. A. Virzi, Computer modelling of heat transfer in Czochralski silicon crystal growth, J. Cryst. Growth 112 (1991) 699e722. € rnstein, vol. 17, SpringerK.-H. Hellwege, O. Madelung (Eds.), Landolt-bo Verlag, Berlin, 1984. M. Watanabe, Talk at ICCG-16, 2010. Beijing. C. Kudla, A.T. Blumenau, F. Büllesfeld, N. Dropka, C. Frank-Rotsch, F. Kiessling,
[26] [27] [28]
[29] [30] [31] [32]
9
O. Klein, P. Lange, W. Miller, U. Rehse, U. Sahr, M. Schellhorn, G. Weidemann, M. Ziem, G. Bethin, R. Fornari, M. Müller, J. Sprekels, V. Trautmann, P. Rudolph, Crystallization of 640 kg mc-silicon ingots under traveling magnetic field by using a heater-magnet module, J. Cryst. Growth 365 (2013) 54e58. I. Rasin, S. Succi, W. Miller, A multi-relaxation lattice kinetic method for passive scalar diffusion, J. Comput. Phys. 206 (2) (2005) 453e462. http://www.mcs.anl.gov/petsc/. C. Herring, The use of classical macroscopic concepts in surface-energy problems, in: R. Gomer, C.S. Smith (Eds.), Structure and Properties of Solid Surfaces, The University of Chicago Press, 1953, pp. 5e41. G.S. Rohrer, Grain boundary energy anisotropy: a review, J. Mater. Sci. 46 (2011) 5881e5895. V.V. Voronkov, Kristallografiya 17 (5) (1972) 909e917 (in Russian). T. Duffar, A. Nadri, On the twinning occurrence in bulk semiconductor crystal growth, Scr. Mater 62 (12) (2010) 955e960. H.K. Lin, C.W. Lan, Revisiting the twinning mechanism in directional solidification of multi-crystalline silicon sheet, Acta Mater 131 (2017) 1e10.
Contents lists available at ScienceDirect
Acta Materialia journal homepage: www.elsevier.com/locate/actamat
Full length article
Micro structures in the grain evolution during solidification of silicon: Phase field calculations W. Miller a, *, A. Popescu b a b
Leibniz Institute for Crystal Growth (IKZ), Max-Born-Str.2, 12489 Berlin, Germany Faculty of Physics, West University of Timisoara, Bd. V. Parvan 4, 300223 Timisoara, Romania
a r t i c l e i n f o
a b s t r a c t
Article history: Received 3 March 2017 Received in revised form 8 August 2017 Accepted 10 August 2017 Available online 12 August 2017
Two dimensional phase field simulations have been performed to study the influence of the growth kinetics and the surface energy on the growth behaviour of grains during solidification of Si. In particular, we studied the groove between two grains as recently observed by in situ observation [1]. Furthermore, we performed computations for different S boundaries and discuss the interplay between equilibrium of interfacial energies and growth kinetics. © 2017 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved.
Keywords: A1. Computer simulation A1. Directional solidification A1. Growth models B2 Semiconducting silicon
1. Introduction Despite the great importance of direct solidification of silicon and the large number of investigations there are still many open questions concerning the grain evolution. Post-mortem analysis of the grain structure is one way for a better understanding of the growth kinetics. Most often these analyses are performed on wafers giving some statistics at a certain height but not the continuous evolution in time. Analysis in vertical direction has been performed e.g. by Prakash et al. [2] and more recently by Lin et al. [3]. However, only the situation after solidification and cooling can be analysed. A better path towards understanding the details is in-situ investigation which became recently available by improved synchrotron X-ray imaging to visualize the melt/solid interface of silicon [4]. Thus, it was possible to follow the creation and disappearance of facetted/facetted grooves [1]. For an understanding of the basic mechanisms a simulation on the microscopic length scale is required, which includes the effect of interface curvature. Here, phase field methods are the preferential ones. Chen et al. studied the grain competition in an undercooled melt with a simple anisotropy function for the solid/ melt interfacial energy. Cantù et al. and Miller et al. studied a more
* Corresponding author. E-mail address: [email protected] (W. Miller). http://dx.doi.org/10.1016/j.actamat.2017.08.025 1359-6454/© 2017 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved.
advanced anisotropy function but simplified growth kinetics [5,6]. All aforentioned calculations have been two-dimensional approaches. Lin & Lan used a three-dimensional model, still with simplified anisotropy function and growth kinetics [7]. So far, all phase field calculations were based on the model of Warren et al. for computing the grain growth [8,9]. Also in this paper we use such a model and discuss in the end the limitations of this approach concerning the grain boundary orientation. Compared to our previous investigations [5,6] we update the solid/melt interfacial energy according to recent numerical and experimental results [10,11]. We also set the growth kinetics based on the results of molecular dynamics calculations [12]. The much more realistic input parameters enabled us to compare our calculations with the results of in situ investigations by Tandjaoui et al. [1]. 2. Interfacial energies In this paper we consider the (110) plane of silicon for 2D calculations (see Fig. 1). We need two kinds of interfacial energies for our computations: that of the melt/crystal interface and those of the grain/grain boundaries. In particular, we need the melt/crystal interface energy as a function of the orientation. Recently, the results of a molecular dynamics study [10] and of an experiment with in-situ observation [11] have been published. In both studies the equilibrium shape in the (110) plane was the result (see Fig. 2). The equilibrium shape represents the relative change of the interface
2
W. Miller, A. Popescu / Acta Materialia 140 (2017) 1e9
Fig. 1. Plane for calculations.
Fig. 2. Equilibrium shape for interface energy of Eq. (1) (blue line) in comparison with result of MD simulation (left [10]) and experiment (right [11]). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
energy with respect to orientation. Absolute value of interface energies were computed earlier for three orientations by means of molecular dynamics method [13]. Apte employed the cleaving wall method to extract the interfacial energies as follows: 0:42±0:02 J=m2 , 0:34±0:02 J=m2 , 0:35±0:03 J=m2 for the {100}, {111}, and {110} planes, respectively. In the computations we need to represent the interface energy by a function, which is differentiable twice with respect to orientation. In order to recover an equilibrium shape close to the observed ones we made the following. appoximation:
.
n
g ¼ 0:335 J m2 1: þ 0:16ð1 tanhð20ðw 50 ÞÞ þ 0:03ð1
for {100} and {110} planes and two deep for {111} and {449} planes. Grain boundary energies for silicon were computed by Kohyama and co-workers using a bond counting method [15]. More recently, molecular dynamics have been employed to obtain the energy at certain S boundaries [16]. In Fig. 4 we show the data points of the two kinds of calculations in the ð111Þ plane. The original values of Kohyama et al. have been scaled by 0.65 to adapt them to the values of the molecular dynamics calculations. Thus, they are consistent with the absolute values of the melt/solid interfacial energy, which are also based on the results of molecular dynamics calculations. For small misorientations one can apply the Read-Shockley theory, Dq 1 ln Dq where gGB ¼ g0GB Dq , where Dq and Dqm are the Dq m
tanhð18ð60 wÞÞ 0:02 sin6 w 1 0:05 cos8 ðwÞ o 0:04 cos128 ðw 34 Þ : (1) The equilibrium shape given by this equation is shown in Fig. 2 as a blue line. Besides the equilibrium shape the three absolute values for the interface energy as given above were a criterium for setting up Eq. (1). We almost match the values given by molecular dynamics calculations as Eq. (1) yields 0:425 J=m2 ({100}), 0:337 J=m2 ({111}), and 0:348 J=m2 ({110}). The interface energy as a function of the orientation according to Eq. (1) is shown in Fig. 3. The curve is similar to the curve derived from bond density calculations [14]. We have four minima: two flat
m
missorientation and the saturation, respectively. In Fig. 4 this energy is shown for g0GB ¼ 0:68 J=m2 and Dqm ¼ 28o as a light grey curve. For the computations we use the values of the blue curve. This is discussed later in the context of the phase field model we use for computing the phase transition.
3. Growth kinetics The growth velocity of a melt/solid interface is given by
vn ¼ bðDTÞDT;
(2)
where DT is the undercooling at the considered point of the interface. The kinetic coefficient b depends on the growth
W. Miller, A. Popescu / Acta Materialia 140 (2017) 1e9
3
Fig. 5. Kinetic coefficient as a function of the missorientation from the 〈111〉 direction. Red diamonds are the data obtained by molecular dynamics calculations [12]. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 3. Melt crystal interfacial energy as given by Eq. (1).
Below 1:2o either Eq. (3) or Eq. (4) for dislocation or 2D nucleation is applied. 4. Numerical methods 4.1. Phase transition The phase transition is treated by the phase-field method as explained in Ref. [5]. It is based on the model of Warren et al. for computing the phase transition and evolution of grains [8,9]. Here, we shortly repeat the main equations. The evolution equations for the order parameter of solid (f ¼ 1) and liquid phase (f ¼ 0) and for the grain orientation (p=2 < j p=2) are:
tf vt f ¼ x2f V2 f þ 4Wfð1 fÞðf 1=2Þ DHV
þ Fig. 4. Grain boundary energy as a function of the misorientation. The red routes are the values obtained by Kohyama with the bond-counting method [17]. The original energies were scaled by 0.65. The two triangles are from molecular dynamics calculations and are taken as a reference for the data of Kohyama [16]. The blue curve represents the curve given by Eq. (13). The light grey curve is the grain boundary energy according to the Read-Shockley theory. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
mechanism. For growth on atomically rough surfaces b does not depend on DT. For Si(100) it was computed by means of molecular dynamics as bð100Þ ¼ 0:12 m=ðs KÞ [12]. At the 〈111〉 facet either a dislocation induced growth can occur, where the kinetic coefficient obeys:
bðDTÞ ¼ bdisl DT;
(3)
or a 2D nucleation, where b is given by
bðDTÞ ¼ aA expfaB =DTgDT 1=3 m=ðs KÞ:
2 Tm
For the {111} plane experimental and numerical attempts have been made to obtain the parameters aA and aB [18e20]. Once a nucleus exists a step growth starts. This type of growth has been investigated by molecular dynamics calculations [12]. The data can be interpolated by lines in two regimes as shown in Fig. 5.
(5)
2sfjVjj x2j fjVjj2 s þ x2j Vj ; f2 tj Vt j ¼ V f2 jVjj
(6)
where xf and xj are the phase field parameters related to the liquid/ solid and grain/grain energy, respectively. DHV and Tm are the latent heat per volume and the melting-point temperature, respectively. W is the barrier for the solid/liquid phase transition at T ¼ Tm and s is a coupling constant. tf is defined by the kinetic coefficient, which was discussed in the previous section. tj is the relaxation time for the grain-grain phase transition and is chosen ten times smaller than tf . Eqs. (5) and (6) are solved by finite difference and a kinetic method, respectively [5]. The solid-liquid interface energy is given by
pffiffiffiffiffiffiffiffi
1
(4)
30f2 ð1 fÞ2 ðT Tm Þ
gsl ¼ pffiffiffiTm xf 2W ; 3 2
(7)
and the grain-grain boundary energy is related to this via
gGB ¼ 2ssl 1 f3min :
(8)
4
W. Miller, A. Popescu / Acta Materialia 140 (2017) 1e9
fmin is the minimum value of f between two grains. From a 1D analysis for T≡Tm one obtains [8]:
fmin ¼ 1
Dq ; Dqc
(9)
where the critical difference in crystallographic orientations is
Dqc ¼
pffiffiffiffiffiffiffiffi 2W xf 1 ¼ : facs facx s
(10)
facx ¼ xj =xf is a model parameter. We chose facx ¼ 8=15 as it was pffiffiffiffiffiffiffiffi done by Warren et al. [8]. facs ¼ s=ð 2W xj Þ defines the stiffness of the grain-grain boundary. For a fixed facs the grain boundary energy as given by Eq. (8) obeys a similar behaviour as those of the Read-Shockley model:
gGB ¼ E0
Dq Dq 1 ln Dqm Dqm
:
(11)
In the Read-Shockley model the grain boundary energy saturates at E0 for large misfits (larger than the material parameter Dqm ). However, gGB exhibits an oscillating behaviour as discussed in section 1 (see Fig. 4). We adapt to this behaviour by varying facs as a function of the misfit angle. In particular, we divide the misfit into two regimes: up to a misfit of Dq ¼ 20o the grain boundary energy can be described by Eq. (11). Beyond we use the values given in Ref. [17] to define the grain boundary energy. In particular the grain boundary energy is given by:
! . Dq 3 2 gGB ðDqÞ ¼ 0:67 J m 1 1 ;
Dqc
(12)
where Dqc is given by:
Dqc ¼
8 > < > : 20o þ
Dq
2
4:47o
30o þ 0:5Dq sinð28Dq 280o Þ
if Dq 20o if Dq > 20o (13)
We have performed computations for a 2D domain with a uniform grid (see Fig. 6). At the bottom we have two grains with a solid/melt interface at x ¼ xint , where x ¼ 0 is at the top and x ¼ xend at the bottom. The width of the diffuse interface (0:05 f 0:95) is set ot be 10 grid points. 4.2. Computational domain and heat transport
R cool Gmelt T
;
transport equation
vT vf DHV ¼ VðDT VÞT þ vt vt rcp
(15)
is solved by a kinetic scheme [26]. In the diffuse interface the ð1 fÞ þ Dsolid f. thermal diffusivity is given as DT ¼ Dmelt T T As sketched in Fig. 6 we start from an artificially initial configuration of grain boundary and solid/melt interface shape. In the end we like to obtain a steady state situation, where the interface is moving without changing its shape and the grain boundary keeps its orientation. In order to reach this state a large number of time steps is required and the solid/melt interface is moving a significantly distance in the meantime. In order to keep the configuration as close to the original one and allow a large number of time steps we perform a rezoning once the average interface moved by one grid spacing [5]. In the solid part the last line of grid points is leaving the computational domain while in the melt part a new line enters. All variables are set to the same value at the adjacent line of grid points. When rezoning the temperature at the boundary x ¼ 0 is reset to its original value Tbc ðt ¼ 0Þ. Only in the steady-state situation the interface is moving approximately with the velocity vn based on the given cooling rate. In the mean time the interface velocity might be quite different because of effects of growth kinetics or changes in the grain boundary. Thus, solving everything in a moving reference frame is not useful. The entire code was parallelized by using the PETSC library version 3.1 [27].
5. Results and discussions
We try to mimic the typical process in a solidification experiment as known from both large scale production experiments (see e.g. Ref. [25]) or small scale research experiments (see e.g. Ref. [4]). The procedure of setting the boundary condition is described in detail in Ref. [5]. The computational domain with boundary conditions is sketched in Fig. 6. The heat flux at the bottom boundary is given by:
j ¼ Dmelt rcp Gmelt þ DHV T T
Fig. 6. Domain for computations with thermal boundary conditions.
(14)
The region around the three phase junction of grains and melt is a complex system, where the interfacial energy effects are coupled with the local temperature field and the growth kinetics. Before discussing the numerical results we discuss the situation at the three phase junction with respect to the interfacial energies without considering the growth kinetics. At the three phase junction between two grains and the melt the interfacial forces between the grains, grain 1 and melt, and grain 2 and melt should be balanced in equilibrium. For the anisotropic system this balance is given by the following equation [28]:
where Dmelt , r, and cp are the thermal diffusivity, the density, and T the specific heat. The cooling rate R cool and the temperature
3 X
are the two process parameters. In exgradient in the melt Gmelt T periments, the latter is typically known only indirectly by
j¼1
¼ R cool =vn . The heat measuring the growth velocity vn : Gmelt T
!
gj b j þ
vgj ! nj ¼ 0 vw
(16)
For a groove the situation is schematically drawn in Fig. 7. In 2D we end up with:
W. Miller, A. Popescu / Acta Materialia 140 (2017) 1e9
Fig. 7. Angles and orientations in Young-Herring's force balance.
vgGB vgsl þ gsl ðq1 Þsinðq1 Þ þ cosðq1 Þ vw vw q1 vgsl cosðq2 Þ ¼ f1 gsl ðq2 Þsinðq2 Þ vw q2
(17)
vgsl sinðq1 Þ vw q1 vgsl sinðq2 Þ ¼ f2 : þgsl ðq2 Þcosðq2 Þ vw q2
(18)
gGB þ gsl ðq1 Þcosðq1 Þ
In Eqs. (17) and (18) we use as a more general version with f1 and f2 on the right hand side. In the equilibrium where all forces cancel f1 ≡f2 ¼ 0. In all phase field computations we have gGB independent vgGB vw ≡0.
The equations of the orientation of the grain boundary, i.e. can be used in different ways. On one hand we can compute the angles q1 and q2 for given grain boundary energy. This enables us to countercheck the results of the numerical calculations and see if they obey the equilibrium or if there is an influence of growth kinetics. The equations can be used also in the opposite way. Knowing q1 and q2 e.g. from experimental results the corresponding grain boundary energy can be computed (see e.g. Ref. [29]). The latter we did for the computation of a facetted/facetted groove as observed during in-situ observation of solidification. 5.1. Computation for facetted/facetted interface Firstly, we consider the case of facetted interfaces at both sides
5
of the groove. We refer to a case obtained by Tandjoui et al. in an experiment with in-situ observation [1]. They used synchrotron Xray radiation to visualize the solid/melt interface exploiting the differences in density of molten and solid silicon [4]. The solidification took place in a pyrolithic boron nitride crucible with a diameter of 6 mm applying a cooling rate of R cool ¼ 0:5 K=min. The growth velocity could be directly obtained from the in-situ measurement yielding vg ¼ 9:2 mm=s. Assuming that the cooling rate at the interface is the same, the temperature gradient can be computed as GT ¼ R cool vg , which yields GT ¼ 900 K=m. In particular, we consider the groove as shown on the right hand side of Fig. 8 and apply our numerical approach to compute the interface shape and grain boundary orientation. As no EBSD images are available for the crystal of this experiment we do not know the orientation of the grains. However, we know that the facets are {111} planes. We assume that the imaging plane in the experiment was a {110} plane to be consistent with our 2D model. According to Fig. 1 there are always two possibilities for the nearest {100} plane. For the facet of left grain (grain 1) the 〈100〉 orientation is either 50o or 20o off the main growth direction. For the other facet the possibilities are 10o or 80o. Thus, in total four combinations are possible. In Fig. 8 the configuration for 20o and 80o is shown. In case that the 〈100〉 direction of grain 1 is 50 off with respect to the main growth direction there is a {111} facet just 4.7 off the flat interface as shown in Fig. 9. This facet will be formed, which results in an undercooling at this facet and also at the left side of the groove. Thus, there is no chance to establish a long facet at the groove and a deep groove as observed in the experiment. Therefore, we computed only cases with 20o off for the left grain. In the following we denote the case with 20o and 10o (Dq ¼ 10o ) as case B and the case with 20o and 80o (Dq ¼ 60o ) as case A. In the experiment a depth of dgroove ¼ 240 mm was observed (see Fig. 8). Using this value the undercooling at the bottom can be ¼ 0:22 K. This suggest a dislocation estimated as DT ¼ dgroove Gmelt T growth mechanism rather than 2D nucleation. Therefore, we use
bdisl ¼ 1 105 m s1 K2 . Beyond an angle of 1:2o off the 〈111〉 direction we apply the step growth rules as shown in Fig. 5. This holds up to 25o , when it reaches the value for rough surfaces (b ¼ 0:12 m s1 K1 ). We performed runs for case A and case B for a domain of 250 250 mm2 with a lattice spacing of Dx ¼ 0:5 mm (parameter see Table 1). Thus, the width of the diffuse interface was wf ¼ 5:0 mm. As we know the geometry of the groove from experiment we can compute the grain boundary energy by Eq. (18). In fact, we used Eqs. (17) and (18) to compute q1 and q2 as a function of the grain boundary energy gGB. We start with case A: the target values for both q1 and q2 are 30o . For a grain boundary energy of gGB ¼ 0:59 J=m2 the computation yields q1 ¼ 28:9o and q2 ¼ 30:6o . According to Fig. 5 for these orientations we have a dislocation growth mechanism. The scaled value of Kohyama for the misfit of
Fig. 8. Left: Experimental observation by Tandjaoui et al. [1] (part of Fig. 2). Right: Orientations in the facetted/facetted groove.
6
W. Miller, A. Popescu / Acta Materialia 140 (2017) 1e9
Fig. 9. Result of computation for the case that the 〈100〉 direction of grain 1 is 50 off with respect to the main growth direction. A facet is formed at the interface, causing an undercooling in this region, which prevents the formation of a deep groove.
Table 1 Material parameters of Si used in simulations. thermal diffusivity in melt
Dmelt T
2:6 105 m2 =s
thermal diffusivity in solid
Dsolid T DHV rcp
9 106 m2 =s
[22]
4:22 109 J=m3
[23]
2:2 108 J=ðK m3 Þ 1683 K
[24]
latent heat specific heat melting-point temperature Gibbs-Thomson coefficient
Tm
G
[21]
1:8 107 K m
Dq ¼ 60o is gGB ¼ 0:65 J=m2 and thus it is slightly larger than the
computed one. In the phase field computations we used a fixed value of Dqcrit ¼ 116o representing gGB ¼ 0:59 J=m2 . As the initial configuration for the run we use the geometry known from the experiment, i.e. a groove with 〈111〉 facets on both sides a grain boundary with a 15 tilt with respect to the main growth direction. After consolidation of the temperature field and establishing the curved interface at the upper boundaries of the groove a steady state growth of the crystal was observed. The depth of the groove in this steady motion was dgroove ¼ 28 mm. The interface shape at two different time steps is shown in Fig. 10 demonstrating that there is no further change in the interface shape. The temperature profile along a line crossing the bottom in the groove is shown in the right part of Fig. 11 giving an undercooling of DT ¼ 0:112 K at the bottom of the groove. This is much larger than the estimation using Gmelt T and the depth of the groove which yields DT ¼ 0:025 K. We will come later back to this point. Here we like to note that the depth in the calculation is much smaller than that observed in the experiment (dgroove ¼ 240 mm). The reason for the discrepancy is the limitation of undercooling at the bottom of the groove in the numerical calculation due to the spatial resolution. At the bottom we have curved interface of curvature k. At that point the melting point temperature Tm is reduced by DTcurv ¼ 2Gk, where G ¼ gDHV =Tm is the Gibbs-Thomson coefficient. For a facetted/facetted groove k should be very large. However, it is limited by the thickness of the diffuse interface: kmax z2=wf . In the case of Dx ¼ 0:5 mm we have wf ¼ 5 mm and thus DTz0:15 K. This corresponds approximately with the undercooling in the computational run. We increased the resolution to Dx ¼ 0:2 mm taking a domain with 2000 2000 grid points. Now we obtained a depth of dgroove ¼ 103 mm with an undercooling of DT ¼ 0:285 K. Still the depth is smaller than the
Fig. 10. Interface lines in Case A for two time steps. The system is in a steady-state mode. The angles of the facets and the grain boundary are in the range of the experimentally observed ones (see top right sketch). The colors of the lines indicate the kinetic coefficient: red and blue means large and small value, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
one observed in the experiment. For a more realistic computation a local grid refinement is required. In order to estimate the undercooling at the bottom of the groove in the experiment we extrapolated the computed undercoolings for two different depth of the groove to the on measured in the experiment. Here, we took into account the limit for a vanishing groove: dgroove ¼ 00DT ¼ 0. We obtained DTz0:6 K. This undercooling is much smaller than required for 2D nucleation. It is much more likely that the growth on the facets was induced by a dislocation. Using bdislocation ¼ 2:8 104 m s1 K1 given by Voronkov [30] the necessary undercooling for vn ¼ 9:2 mm=s is 0.03 K. This is appoximately the undercooling at top of the facets. Thus, we expect that the growth on the facets was induced by at least one dislocation on each facet. The undercooling has also an impact on the probability of twin nucleation on the {111} facets. In the experiments it was observed that the facetted/facetted grooves most often disappear by the nucleation of a twin. According to the estimations by Duffar and Nadri [31] a reasonable probability for twinning is reached for undercoolings larger than 1.5 K, which is obviously much larger than observed in experiments. Recently, Lan and co-workers developed an advanced model for a twin seed, yielding a probability of about 1 109 already for an undercooling of DT ¼ 0:5 K. At DT ¼ 0:6 K it is already about 1 107 [32]. Therefore, the occurrence of twinning in the experiments is explainable by this model in conjunction with the temperature computation at the groove. In case B the misorientation is small ðDq ¼ 10o Þ and in the ReadShockley regime. From Fig. 4 one can see that gGB < 0:45 J=m2 . Performing the computations with Eqs. (17) and (18) we obtain gGB ¼ 0:58 J=m2 for angles q1 ≡q2 ¼ 28:8o . This means that the computed value using the geometry and our approach for the melt/ solid interface energy is much larger than the one from our model for the grain boundary energy. Either our models for the interface energies are not appropriate or this configuration was not realized in the experiment. Nevertheless, we performed runs for this case using a fixed value of Dqcrit ¼ 20o , representing gGB ¼ 0:58 J=m2 for Dq ¼ 10o . At the end we like to make some comments on the growth kinetics. Looking in more detail on the orientation of the two facets one observes that the angle with respect to the 〈111〉 direction is
W. Miller, A. Popescu / Acta Materialia 140 (2017) 1e9
7
Fig. 11. Temperature profile along line through deepest point of groove. Left: Line in the configuration of run. The light blue curve indicates the isoline of the melting point temperature. Right: temperature profile in red and phase field profile in blue. Regions a,b, and c mark the regions of groove, interface, and grain boundary, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
about 1:2o , i.e. near the boundary where the dislocation growth changes to step growth. There is step flow growth all along the facets except at the top end and at the bottom. 5.2. S boundaries As we mentioned already in all our calculations we use a unique value for the grain boundary energy for a given misfit. Therefore, the grain boundary energy does not depend on the orientation of the grain boundary. Of course, this is a simplification of the reality. Especially for coincidence lattice boundaries (S boundaries) one orientation is energetically favoured with respect to all others. This energy minimizing is valid for the bulk. The question is what happens at the three phase junction, where the orientation of the grain boundary to fulfill the force balance might be completely different from the lattice coincidence orientation (see Fig. 12). The most prominent S boundary is the S 3ð111Þ twin boundary. All computations agree that its boundary energy is negligible. Thus, it is most likely that two grains in this twin relation always form a (111) plane at the interface regardless of the local situation at the interface. It is also expected that the interface is almost flat with
only a small groove at the three-phase junction. However, the question is if this holds for all orientations of the growth direction with respect to the S 3ð111Þ twin boundary. We used Eq. (18) to obtain for a given q1 the corresponding angle q2 . In Fig. 13 the sketch of angles and the angle of the groove as a function of the angle between grain boundary and growth direction is shown. It is 4i ¼ 90o qi (i ¼ 1; 2) and 42 is negative in the figure. The growth direction is assumed to be symmetric with respect to the groove and in such a way wgrain is defined. A value wgroove > 180o means not a groove but a hillock, which will be unstable during growth. In Fig. 13 only the results are shown where the resulting interfaces at grain 1 and grain 2 are both non-missing interfaces. If one of the interfaces is a plane of a missing direction according to Eq. (1) (see Fig. 12) the situation is rather unstable and the values for wgroove are fluctuating. In Fig. 13 we show all data points used and it can be clearly seen that there are certain ranges of wgrain where there are no points, i.e. at least one of the direction of the planes is missing one. Next, we consider the S19 boundary, which has a grain boundary energy of 0:52 J=m2 (at the coincidence plane (133)). We performed runs for two different growth directions as shown in Fig. 14. In Case 1 there is only a small deviation between the computed grain boundary (white line) and the (133) coincidence plane (dotted orange line). Both sides of the groove are 〈111〉 facets. In Case 2 the deviation is larger. The groove is not as deep as in Case 1, because on both sides the interface is slightly off the (111) plane exhibiting a step growth mode with a larger kinetic coefficient b. Please note that in our phase field model there is no variation of the grain boundary energy with orientation and thus, there is no energy minimization with respect to this energy. In the future we will adapt the multi phase model for our purposes, which allows the definition of an anisotropy in the grain boundary energy. 6. Conclusions
Fig. 12. Stiffness of melt/solid interface. There are three regions of missing directions. The region around w ¼ 20o is enlarged for better visibility.
We applied a phase field model to compute the grain evolution during the solidification of silicon on the microscopic scale. We could recover the experimental observation by Tandjaoui et al. [1] of a groove with two 〈111〉 facets. The main issue is a very fine resolution at the bottom of the groove to reflect a large curvature at that point sustaining the undercooling caused by growth kinetics. The undercooling at the bottom is about 0:6 K. According to the recent work of Lin et al. [32] this undercooling is enough for a
8
W. Miller, A. Popescu / Acta Materialia 140 (2017) 1e9
Fig. 13. Computation of the angle of the groove wgroove as a function of the grain boundary with respect to the growth direction. For a given f1 we compute f2 according to Eq. (18). gGB Using Eq. (17) one can compute the required vvw .
Fig. 14. S19 grain boundary for two orientations. The dotted orange line indicates the S19 boundary, i.e. the (133) coincidence plane. The yellow line represents the isoline of Tm . In Case 1 the computed grain boundary is very close to the (133) coincidence plane. In Case 2 the deviation is much larger. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
reasonable probability of twin nucleation. This corresponds with the in-situ observations in the experiments [1]. We also investigated cases with S boundaries. Here, one would expect that the grain boundary is the coincidence plane of the two grains. However, the situation at the three phase junction with respect to interfacial energies and in addition growth kinetics might prefer a different orientation of the grain boundary. Part of the story is the weakness of the currently used phase field model, which does not allow to include an orientation dependent grain boundary energy. Therefore, we will change to a multi phase field model in the future. Acknowledgments
[5] [6] [7]
[8]
[9] [10]
[11]
W.M. likes to acknowledge fruitful discussions with Nathalie €el and her research group. Mangelinck-No
[12]
References
[13]
[1] A. Tandjaoui, N. Mangelinck-Noel, G. Reinhart, B. Billia, T. Lafford, J. Baruchel, Investigation of grain boundary grooves at the solideliquid interface during directional solidification of multi-crystalline silicon: in situ characterization by X-ray imaging, J. Cryst. Growth 377 (2013) 203e211. [2] R.R. Prakash, K. Jiptner, J. Chen, Y. Miyamura, H. Harada, T. Sekiguchi, Grain boundary interactions in multicrystalline silicon grown from small randomly oriented seeds, Appl. Phys. Express 8 (3) (2015) 035502. [3] H.K. Lin, M.C. Wu, C.C. Chen, C.W. Lan, Evolution of grain structures during directional solidification of silicon wafers, J. Cryst. Growth 439 (2016) 40e46. €l, G. Reinhart, J.-J. Furter, B. Billia, T. Lafford, [4] A. Tandjaoui, N. Mangelinck-Noe J. Baruchel, X. Guichard, Real time observation of the directional solidification
[14] [15]
[16]
[17] [18]
of multicrystalline silicon: X-ray imaging characterization, Energy Procedia 27 (2012) 82e87. G. Cantù, A. Popescu, W. Miller, Grain growth of silicon, Acta Mater 60 (2012) 6755e6761, http://dx.doi.org/10.1016/j.actamat.2012.08.048. W. Miller, A. Popescu, G. Cantù, Solidification of multicrystalline silicon simulation of micro structures, J. Cryst. Growth 385 (2014) 127e133. H.K. Lin, C.W. Lan, Three-dimensional phase field modeling of silicon thin-film growth during directional solidification: facet formation and grain competition, J. Cryst. Growth 401 (2014) 740e747. J.A. Warren, R. Kobayashi, A.E. Lobkovsky, W.C. Carter, Extending phase field models of solidification to polycrystalline materials, Acta Mater 51 (2003) 6035. J.A. Warren, R. Kobayashi, W.C. Carter, Modeling grain boundaries using a phase-field technique, J. Cryst. Growth 211 (2000) 18e20. G. Barinovs, A. Sabanskis, A. Muiznieks, Study of silicon crystal surface formation based on molecular dynamics simulation results, J. Cryst. Growth 391 (2014) 13e17. X. Yang, K. Fujiwara, K. Maeda, J. Nozawa, H. Koizumi, S. Uda, Crystal growth and equilibrium crystal shapes of silicon in the melt, Prog. Photovolt. Res. Appl. 22 (2014) 574e580. D. Buta, M. Asta, J.J. Hoyt, Kinetic coefficient of steps at the Si(111) crystalmelt interface from molecular dynamics simulations, J. Chem. Phys. 127 (07) (2007) 074703. P.A. Apte, X.C. Zeng, Anisotropy of crystal-melt interfacial free energy of silicon by simulation, Appl. Phys. Lett. 92 (22) (2008) 221903. P.J. Hesketh, C. Ju, S. Gowda, E. Zanoria, S. Danyuk, Surface free energy model of silicon anisotropic etching, J. Electrochem. Soc. 140 (4) (1993) 1080e1085. M. Kohyama, R. Yamamoto, M. Doyama, Reconstructed structures of symmetrical 〈011〉 tilt grain boundaries in silicon, Phys. Stat. Sol. B 138 (1986) 387e397. X. Chen, L. Xiong, A. Chernatynskiy, Y. Chen, A molecular dynamics study of tilt grain boundary resistance to slip and heat transfer in nanocrystalline silicon, J. Appl. Phys. 116 (2014) 244309. M. Kohyama, Computational studies of grain boundaries in covalent materials, Modelling Simul, Mater. Sci. Eng. 10 (2002) R31eR59. K.A. Jackson, Response to: some remarks on the undercooling of the Si(111)
W. Miller, A. Popescu / Acta Materialia 140 (2017) 1e9
[19]
[20] [21]
[22] [23] [24] [25]
facet and the “Monte Carlo modeling of silicon crystal growth” by Kirk M. Beatty & Kenneth A. Jackson, J. Cryst. Growth 211 (2000), 13 by W. Miller, J. Crystal Growth 325 (1) (2011) 104. W. Miller, Some remarks on the undercooling of the Si(111) facet and the “Monte Carlo modeling of silicon crystal growth” by Kirk M. Beatty & Kenneth A. Jackson, J. Cryst. Growth 211 (2000) 13. J. Crystal Growth 325 (1) (2011) 101e103. K.M. Beatty, K.A. Jackson, Monte Carlo modeling of silicon crystal growth, J. Cryst. Growth 211 (2000) 13e17. Y. Inatomi, F. Onishi, K. Nagashio, K. Kuribayashi, Density and thermal conductivity measurements for silicon melt by electromagnetic levitation under a static magnetic field, Int. J. Thermophys. 28 (1) (2007) 44e59. A. Virzi, Computer modelling of heat transfer in Czochralski silicon crystal growth, J. Cryst. Growth 112 (1991) 699e722. € rnstein, vol. 17, SpringerK.-H. Hellwege, O. Madelung (Eds.), Landolt-bo Verlag, Berlin, 1984. M. Watanabe, Talk at ICCG-16, 2010. Beijing. C. Kudla, A.T. Blumenau, F. Büllesfeld, N. Dropka, C. Frank-Rotsch, F. Kiessling,
[26] [27] [28]
[29] [30] [31] [32]
9
O. Klein, P. Lange, W. Miller, U. Rehse, U. Sahr, M. Schellhorn, G. Weidemann, M. Ziem, G. Bethin, R. Fornari, M. Müller, J. Sprekels, V. Trautmann, P. Rudolph, Crystallization of 640 kg mc-silicon ingots under traveling magnetic field by using a heater-magnet module, J. Cryst. Growth 365 (2013) 54e58. I. Rasin, S. Succi, W. Miller, A multi-relaxation lattice kinetic method for passive scalar diffusion, J. Comput. Phys. 206 (2) (2005) 453e462. http://www.mcs.anl.gov/petsc/. C. Herring, The use of classical macroscopic concepts in surface-energy problems, in: R. Gomer, C.S. Smith (Eds.), Structure and Properties of Solid Surfaces, The University of Chicago Press, 1953, pp. 5e41. G.S. Rohrer, Grain boundary energy anisotropy: a review, J. Mater. Sci. 46 (2011) 5881e5895. V.V. Voronkov, Kristallografiya 17 (5) (1972) 909e917 (in Russian). T. Duffar, A. Nadri, On the twinning occurrence in bulk semiconductor crystal growth, Scr. Mater 62 (12) (2010) 955e960. H.K. Lin, C.W. Lan, Revisiting the twinning mechanism in directional solidification of multi-crystalline silicon sheet, Acta Mater 131 (2017) 1e10.