Anisotropic solid–liquid interface kinetics in silicon: an atomistically informed phase-field model

Molecular dynamic calculations are performed using the SW interatomic potential for silicon [42], which describes the structure of the molten phase realistically, and reproduces the experimental melting point [25, 28, 29, 38]. We calculate thermal properties and interfacial velocities with the widely used LAMMPS code [35] and obtain free energy densities via using the MD++ code [37]. In table 1 we present a summary of the thermal properties obtained using the SW potential along with results form previous simulations and experiments.

Table 1.  Heat capacity, latent heat and melting point of silicon. T_m^MD is the melting point obtained from moving interface simulations, T_m^V is the melting point obtained at constant volume from the adiabatic switching method (see section 3.2).

	${c}_{p,{\rm{solid}}}$	${c}_{p,{\rm{liquid}}}$	${\rm{\Delta }}{H}_{{\rm{m}}}$(kJ mol−1)	TmMD(K)	TmV(K)
	(J mol−1 K−1)	(J mol−1 K−1)
SW-pot.	26.9	28.7	31.7	1683	1697.12
other sim.			30.9 [7]	1682 [17]
Exp.	26–29 [12]	27.2 ± 1.5 [12]	50.25 ± 0.6[12]	1683[42]
Data	(900–1687 K)

We initialize a simulation box containing 4096 atoms in the diamond structure and heat it up. For doing so, we apply a Nosé–Hoover thermostat with a rate of 10¹³ K s⁻¹. After melting occurs, we cool it down with the same rate and calculated the specific heat and latent heat from the average total energy. Clearly, the specific heat and melting point are in good agreement with the experimental finding. The latent heat, in contrast, is low compared to experimental measurements. The reason is that we describe the solid phase and the liquid phase by a single empirical model despite their different bonding mechanisms [37]. The melting point T_m^MD calculated from a simulation cell with solid–liquid phases co-existence is almost exact compared to the experimental value of 1683 K and in good agreement with earlier simulations [14, 17]. If the intersect of free energies calculated by the adiabatic switching method at constant volume is used (see section 3.2) the calculated melting point ${T}_{{\rm{m}}}^{{\rm{V}}}=1697.12\,{\rm{K}}$ is slightly higher. In order to be consistent with the free energy data, we use T_m^V as melting temperature for the phase-field model.

For calculating the free energy density, we use a supercell with 512 atoms and periodic boundary conditions in all three directions. The actual calculation of the free energy is performed using both adiabatic switching and reverse scaling [11, 45] as implemented in the MD++ software package by Ryu and Cai [37]. In order to be consistent with our phase-field model, we modify the code such that an NVT ensemble is used, the initial volume for both phases—solid and liquid—is equal to the equilibrium volume of the crystalline silicon supercell at 0 K. The free energy density is calculated then by dividing the obtained free energies by that same volume. Furthermore, in our approach we include the free energies of the amorphous state calculated by Broughton and Li [7]. The results are presented in figure 2. The melting point used in the phase-field model is the one obtained from this free energy calculation, ${T}_{{\rm{m}}}^{{\rm{V}}}=1697.12\,{\rm{K}}$, which is the temperature at which the amorphous/liquid and solid free energies intersect. We point out that this value does not correspond to the experimental value or to the values obtained from direct interface calculations (see table 1). The reasons for this discrepancy are the NVT ensemble we used for the free energy calculation and a numerical error in the adiabatic switching and reverse scaling methods, used for a large temperature intervals like the one we analyze.

Zoom In Zoom Out Reset image size

Interfacial free energies calculated with the SW potential are available in literature from Apte and Zeng [2], who used molecular dynamics to determine ${\gamma }_{\{100\}},{\gamma }_{\{110\}}$ and ${\gamma }_{\{111\}}$ at the melting point. Their mean values are given in table 2. The most densily packed $(111)$ orientation has the lowest interface energy, while the (110) direction exhibits nearly the same excess energy. This is obviously at odds with the equilibrium shape of Si grains embedded in a melt, which show a Wulff shape with (111) and (100) facets only. Thus, ${\gamma }_{\{110\}}$ is obviously underestimated by the SW-potential. Since the purpose of this study is to devise model parameter, which allows us to directly combine phase-field and molecular dynamics simulations with consistent model parameters corresponding to the SW potential, we adopt these values for our parametrization.

Table 2.  The crystal-melt interface free energies γ at the melting point T_m^V for the orientations {100}, {110} and {111} calculated by Apte and Zeng in [2].

Orientation	{100}	{110}	{111}
$\gamma ({T}_{{\rm{m}}}^{{\rm{V}}})$ in eV Å−2	0.026	0.0218	0.0212

We calculate interface velocities from molecular dynamics simulations of moving planar liquid–solid interfaces with different crystallographic orientations at constant temperature $T\in [800,2000\,{\rm{K}}]$. For this, we initialize a simulation box of about 43 × 43 × 130 Å, like the one shown in the left-hand side of figure 3. Such box contains 12 000–13 000 atoms, depending on the crystal orientation and has periodic boundary conditions in all three directions. The box sizes in x- and y-direction are adjusted in order to obtain a single crystal without lattice defects near the boundaries of the box. We start our simulation with an equilibration phase using a Nosé–Hoover thermo- and barostat for 10 ps at the desired temperature in order to consider thermal expansion of the box. One timestep corresponds to 1 fs. The box dimensions are left free to vary independently of each other. For the case of temperatures below the melting point, as shown in the left-hand side of figure 3, the crystalline part of the box (about $1/12$) is equilibrated at the desired temperature. The remaining atoms are melted at 1000 K above the melting point and then cooled to a temperature near T_m^MD. While melting, the box dimensions in x- and y-direction are fixed, but in z-direction the box is allowed to shrink or expand. Finally, we run the crystallization for some nanoseconds with a global thermostat at the desired temperature. The x- and y-dimensions are fixed again, but not in the growth direction. For the case of temperatures above the melting point the procedure is analogous. In this case, the lower part of the box (about $1/12$) is heated up to the desired temperature and, therefore, is melted, while the upper part of the box is kept crystalline with a temperature near T_m^MD. Then, for some nanoseconds, the complete simulation box is connected to global thermostat at the desired temperature above T_m^MD.

Zoom In Zoom Out Reset image size

To extract the velocity of the interface, one first has to determine its position at each timestep. There are numerous ways to do so (see [28]), for instance, by monitoring the particle density or the atomic potential energy. The observed parameter only has to fulfill the condition that it is sufficiently different in the solid and liquid phase. In this study, we choose the centrosymmetry parameter [24], which can be calculated for each atom within LAMMPS. It is zero for an atom in a perfect lattice, and gives a positive value for disturbed atomic environments. The average of the centrosymmetry parameter over one atomic layer perpendicular to the growth direction (corresponding to some hundreds atoms) results in a centrosymmetry of 8.5–11 for a liquid layer and 12–13.5 for a crystalline layer. We find that the centrosymmetry of melt and crystal is dependent on temperature, as shown in the right-hand of figure 3.

Therefore, we take the mean value of the centrosymmetry of crystal and liquid as the critical value CS_crit in order to distinguish liquid and solid atoms. Using the centrosymmetry method, the isothermal interface velocities are finally determined for certain temperatures in the range 800–2000 K for the SW potential.

The latent heat and the heat capacity determine (together with heat conductivity) how much heat is generated in a crystallizing sample at the moving interface and how fast it is conducted away. It was shown by Monk et al [33] that due to the release of latent heat, the actual interface temperature can differ from the one which is set by the thermostat. Thus, Monk et al proposed to use multiple thermostats, from which each one only sets the temperature for a volume element smaller than 20 Å in thickness. They simulated the scenario for pure $\mathrm{Ni}$. Our temperature calculations during crystallization show a flat temperature profile over the whole simulation box. This indicates that heat is taken away fast enough by one global thermostat and did not influence the crystallization velocity.

Another important feature of silicon is the presence of an amorphous phase, if there is a significant undercooling. This is captured by our isothermal conditions for the moving interface and results in a Vogel–Fulcher type dependence of the interface velocity on temperature [18, 41, 44].

However, in order to feed the PFM with these information, we need an analytical expression for the growth velocity. The growth velocity is described by the product of driving force P and mobility M, which is formulated for an atomically flat solid–liquid interface by Jackson [23] as

$$\begin{eqnarray}&&v=M\cdot P=f\cdot A\cdot \nu \cdot P,\end{eqnarray} \tag{ 4 }$$

where f represents the percentage of favorable growth sites (i. e. steps) on the crystal surface, $A=\sqrt[3]{{\rm{\Omega }}}$ the cube root of the atomic volume, ν the attempt frequency of atom jumps over the interface and P the driving force for crystallization. Within transition state theory [23] the interface velocity is given by the difference between the velocities of crystallization and melting:

$$\begin{eqnarray}&&v={v}_{0}\cdot {\rm{\exp }}\left(-\displaystyle \frac{Q}{{kT}}\right)\cdot \left[1-{\rm{\exp }}\left(-\displaystyle \frac{{\rm{\Delta }}G}{{kT}}\right)\right],\end{eqnarray} \tag{ 5 }$$

where Q is activation energy associated with atomic mobility and ${\rm{\Delta }}G$ is the Gibbs free energy difference between the two phases. The last term in brackets is the thermodynamic driving force for crystallization F, and can be approximated by a series expansion, which we develop upto second order to obtain

$$\begin{eqnarray}&&F=\left(1-\exp \left(-\displaystyle \frac{{\rm{\Delta }}G}{{kT}}\right)\right)\approx \displaystyle \frac{{\rm{\Delta }}G}{{kT}}.\end{eqnarray} \tag{ 6 }$$

In [47], Wilson derives the term $d\cdot \nu$ in equation (4) as $6D/d$ with the diffusion coefficient D. Frenkel [16] refines this expression further by replacing the diffusion coefficient with the Stokes–Einstein relation

$$\begin{eqnarray}&&D=\displaystyle \frac{{kT}}{3\pi \eta r},\end{eqnarray} \tag{ 7 }$$

which describes the diffusion of a spherical particle with radius r in a liquid with viscosity η. As a first approximation, the mobility of a particle in the liquid follows an Arrhenius function. However, from the enthalpy as a function of temperature from the section 3.1, we noted the occurrence of a glass transition. Specifically, when approaching and crossing T_g the mobility of the atoms in the melt is reduced and diffusion is slowed down drastically. The Arrhenius description does not describe this behavior of glass-forming melts. To overcome this problem, Vogel [44] and Fulcher [18] introduced an empirical relation allowing the increase in viscosity when approaching the glass transition

$$\begin{eqnarray}&&\eta ={\eta }_{0}\cdot \exp \left(\displaystyle \frac{A}{T-{T}_{{\rm{VF}}}}\right),\end{eqnarray} \tag{ 8 }$$

where ${\eta }_{0}$ and A are constants and the Vogel temperature, T_VF, lies about 50 K below the glass transition temperature. By replacing the Arrhenius- with the Vogel–Fulcher-expression, we finally obtain

$$\begin{eqnarray}&&v=f\cdot \displaystyle \frac{{\rm{\Delta }}G}{{kT}}\cdot \displaystyle \frac{6D}{d}=\displaystyle \frac{2f\cdot {\rm{\Delta }}{H}_{{\rm{m}}}}{3\pi {\eta }_{0}{rd}\cdot {T}_{{\rm{m}}}^{{\rm{V}}}}\cdot ({T}_{{\rm{m}}}^{{\rm{V}}}-T)\cdot \exp \left(-\displaystyle \frac{A}{T-{T}_{{\rm{VF}}}}\right).\end{eqnarray} \tag{ 9 }$$

The resulting velocity–temperature relationships are depicted in figure 4, where we fit our measurements with the Vogel–Fulcher expression (9).

Zoom In Zoom Out Reset image size

The small velocity of the {111} interface is related to its dense packed structure and low energy, which does not provide favorable sites for the attachment of atoms. Therefore, a nucleation step has to take place before a new {111} layer can grow.

On the contrary, the {110} and {100} interfaces are rough because of the formation of {111} facets. Since the growth in these directions is not nucleation limited, it is faster, which is in agreement with results reported in previous studies [3, 43]. In figure 5 the atomic configuration of an {110} interface under growth conditions is shown. The formation of crystalline {111} can be clearly seen as well as the formation of faulted planes on these facets. As a rule of thumb, the relation between the maximum velocities ${v}_{\{100\}}:{v}_{\{110\}}:{v}_{\{111\}}$ of the fitted curves is about 1:0.8:0.4 for this potential. The difference in growth velocity of {100} and {110} is related to the factor f in equation (4) since all other terms are bulk properties and not orientation dependent. To approach the factor f theoretically, we calculated the density of favorable sites by the density of {111} planes ending at a {100} surface. This is given by

$$\begin{eqnarray*}&&{\rho }_{\{100\}}=\displaystyle \frac{\sin {\alpha }_{\{100\}}}{d},\end{eqnarray*}$$

where d is the distance of {111} planes and ${\alpha }_{\{100\}}$ the angle between the {111} and {100} plane. If one compares ${\rho }_{\{100\}}$ and ${\rho }_{\{110\}}$ a relation of $1:1/\sqrt{2}$ or 1:0.7 is found, which is in rough agreement with the relation derived from the simulation. Alternatively one can count the number of broken bonds per area at the surface, which gives an identical result.

Zoom In Zoom Out Reset image size

Compared to literature values for MD simulations, we find values of 18–20 m s⁻¹ for Stillinger {100} and 9–14 m s⁻¹ for {111} [25, 28, 29], which is in good agreement with the above measured values. From experiments, velocities of 1.6 m s⁻¹ are reported by Kuo [27] and 14 m s⁻¹ by Ohdaira [34], so that we conclude that the results for the SW potential are a good representation of the anisotropic growth velocity of silicon crystals.

A free energy density $F(p,T)$, which has the form of a double-well potential in p, can be established as a polynomial of fourth degree, which is one of the common forms for F. Here, the coefficients may depend on the temperature and we assume the expression

$$\begin{eqnarray*}&&F(p,T)={a}_{0}(T)+{a}_{1}(T)p+{a}_{2}(T){p}^{2}+{a}_{3}(T){p}^{3}+{a}_{4}(T){p}^{4}.\end{eqnarray*}$$

Since the equilibrium states for the bulk free energy density are represented by the two minima of the double-well polynomial, we choose the coefficients ${a}_{0},\,\ldots ,\,{a}_{4}$, such that $(0,F(0,T))$ and $(1,F(1,T))$ are the minima points of F. Then, we equip the free energy density F of the phase-field model with the equilibrium values of the atomistic free energy and via polynomial interpolation, such that we obtain two polynomials, ${f}_{1}(T)$ for the crystalline values in figure 2 at p = 1 and ${f}_{0}(T)$ for the amorphous/liquid values at p = 0. With respect to the condition, that f₁ and f₀ represent the minima of F, the free energy density has the form

$$\begin{eqnarray}&&F(p,T)={f}_{0}(T)+a(T){p}^{2}-2(a(T)+2H(T)){p}^{3}+(a(T)+3H(T)){p}^{4},\end{eqnarray} \tag{ 11 }$$

where $H(T)={f}_{0}(T)-{f}_{1}(T)$. Note here, that $H(T\lt {T}_{{\rm{m}}}^{{\rm{V}}})\gt 0$. For the remaining degree of freedom a(T) in (11) yields:

$$\begin{eqnarray}a(T)\gt \left\{\begin{array}{ll}0 & T\in [0,{T}_{{\rm{m}}}^{{\rm{V}}}]\\ -6H(T) & T\gt {T}_{{\rm{m}}}^{{\rm{V}}}.\end{array}\right.\end{eqnarray} \tag{ 12 }$$

This expression for F fulfills $F(0,T)={f}_{0}(T)$ and $F(1,T)={f}_{1}(T)$ for the minima. Details of the derivation of (11) and (12) are given in the appendix A.1.

Furthermore, the maximum point of F is $(\mu (T),F(\mu (T),T))$ with

$$\begin{eqnarray}&&\mu (T)=\displaystyle \frac{a(T)}{2a(T)+6H(T)}\in (0,1).\end{eqnarray} \tag{ 13 }$$

The calculation of (13) is described in appendix A.2. Hence, as we indicated in section 2, for the line g(p) which is tangent to both of the minima of F, the energy barrier B_kin, that has to be exceeded to get from one equilibrium phase to the other, is the difference of the function values of the maximum $F(\mu )$ and $g(\mu )$. At the melting point T_m^V, the two minima of F have the same function value and thus ${B}_{{\rm{kin}}}({T}_{{\rm{m}}}^{{\rm{V}}})$ is the difference of the function value of the maximum of F and the function value of the minima. Since the function values of the minima are then ${f}_{0}({T}_{{\rm{m}}}^{{\rm{V}}})={f}_{1}({T}_{{\rm{m}}}^{{\rm{V}}})$, with equation (13) the maximum point at T_m^V has the simple expression

$$\begin{eqnarray*}&&(\mu ({T}_{{\rm{m}}}^{{\rm{V}}}),F(\mu ({T}_{{\rm{m}}}^{{\rm{V}}})))=(1/2,{f}_{0}({T}_{{\rm{m}}}^{{\rm{V}}})+a({T}_{{\rm{m}}}^{{\rm{V}}})/16).\end{eqnarray*}$$

Hence, the kinetic barrier at the melting point has the form

$$\begin{eqnarray}&&{B}_{{\rm{kin}}}({T}_{{\rm{m}}}^{{\rm{V}}})=\displaystyle \frac{a({T}_{{\rm{m}}}^{{\rm{V}}})}{16}.\end{eqnarray} \tag{ 14 }$$

The kinetic barrier closely relates to the interface energy γ, which we discuss in section 4.2, where we also determine B_kin and hence with (14) the remaining degree of freedom for the bulk energy.

We derive the gradient energy coefficient σ and the degree of freedom a in the free energy density (11) consistent with the interfacial energies obtained in Apte and Zeng [2], who used molecular dynamic simulations. We note that for the calculation of σ and a, we presently use only the values for the equilibrium state and hence make a independent on the temperature, since there is no literature with values of temperature dependent interface free energy calculated by means of molecular dynamics with SW potential for $\mathrm{Si}$. But after all, the mobility parameter ${M}^{{\rm{PF}}}$ will compensate the missing temperature dependence at this point.

Since Allen and Cahn [1], we know the relations between the interfacial energy coefficient σ in the PFM, the modeled interface thickness ε, the interfacial free energy γ and the modeled kinetic barrier B_kin.

At the melting point $T={T}_{{\rm{m}}}^{{\rm{V}}}$ the relations are

$$\begin{eqnarray*}&&{\gamma }_{\{{hkl}\}}({T}_{{\rm{m}}}^{{\rm{V}}})\approx {\sigma }_{\{{hkl}\}}\sqrt{\displaystyle \frac{1}{2}{B}_{{\rm{kin}}}({T}_{{\rm{m}}}^{{\rm{V}}})},\quad {\varepsilon }_{\{{hkl}\}}\approx \displaystyle \frac{{\sigma }_{\{{hkl}\}}}{\sqrt{2{B}_{{\rm{kin}}}({T}_{{\rm{m}}}^{{\rm{V}}})}}.\end{eqnarray*}$$

For convenience we set

$$\begin{eqnarray}&&\displaystyle \frac{a}{16}={B}_{{\rm{kin}}}({T}_{{\rm{m}}}^{{\rm{V}}})=\displaystyle \frac{{\gamma }_{\{111\}}({T}_{{\rm{m}}}^{{\rm{V}}})}{{\varepsilon }_{\{111\}}},\quad {\sigma }_{\{{hkl}\}}=\sqrt{2{\gamma }_{\{{hkl}\}}({T}_{{\rm{m}}}^{{\rm{V}}})\,{\varepsilon }_{\{{hkl}\}}}.\end{eqnarray} \tag{ 15 }$$

For the crystal orientations {100}, {110} and {111}, Apte and Zeng [2] obtained the crystal-melt interfacial free energy γ, see table 3. We assume for orientation {111} an interface thickness of ${\varepsilon }_{\{111\}}=1\,\mathrm{nm}$ and find a kinetic barrier of ${B}_{{\rm{kin}}}\approx 0.002\,\mathrm{eV}\,{\mathring{\rm A} }^{-3}$. Hence, the interface thickness of orientation {100} results in ${\varepsilon }_{\{100\}}\approx 12.35\,\mathring{\rm A}$ and of {110} in ${\varepsilon }_{\{110\}}\approx 10.29\,\mathring{\rm A}$. In table 3, the values $\gamma ,\varepsilon$, and the resulting model parameters σ and a are given for all three considered orientations. Figure 6 shows the resulting double-well potential at different temperatures.

Zoom In Zoom Out Reset image size

Table 3.  The crystal-melt interface free energies γ at the melting point T_m^V for the orientations {100}, {110} and {111} calculated by Apte and Zeng in [2]. With the help of γ and (15), we calculate the model constants a and σ for the three orientations. For that we choose an interface thickness of $\varepsilon =1\,\mathrm{nm}$ for orientation {111}.

Orientation	{100}	{110}	{111}
${\gamma }_{\{{hkl}\}}({T}_{{\rm{m}}}^{{\rm{V}}})$ in eV Å−2	0.0262	0.0218	0.0212
${\varepsilon }_{\{{hkl}\}}$ in Å	12.35	10.29	10
${\sigma }_{\{{hkl}\}}^{2}$ in eV Å−1	0.8042	0.6702	0.6512
a in eV Å−3	0.0339

In the previous subsections, we derived the bulk energy and the gradient energy of the phase-field model (10) by explicit use of the results from molecular dynamics. In this section, we calculate the mobility parameter ${M}_{\{{hkl}\}}^{{\rm{PF}}}(T)$ of the phase-field model (10), such that the model reproduces the atomistic velocities ${v}_{\{{hkl}\}}^{{\rm{MD}}}(T)$ from section 3.4. But in the phase field, the interface velocity is not a parameter which can be directly incorporated as the interface energy ${\gamma }_{\{{hkl}\}}$ in ${\sigma }_{\{{hkl}\}}$ or the free energy minima in F. Hence, we implement a shooting method, where we vary ${M}_{\{{hkl}\}}^{{\rm{PF}}}(T)$ until all required conditions are fulfilled for a fixed temperature and orientation. This procedure is repeated for all three considered directions (100), (110) and (111), where we make measurements for the temperatures $T=800,850,900,\ldots ,1950,2000\,{\rm{K}}$.

For the shooting method, we first define boundary values for (10), such that we have crystal material at the left boundary, and liquid material at the right:

$$\begin{eqnarray}&&0={M}_{\{{hkl}\}}^{{\rm{PF}}}(T)\,{\sigma }_{\{{hkl}\}}^{2}\displaystyle \frac{{\partial }^{2}p}{\partial {x}^{2}}-{M}_{\{{hkl}\}}^{{\rm{PF}}}(T)\,\displaystyle \frac{\partial F}{\partial p}-\displaystyle \frac{\partial p}{\partial t},\,p(-\infty ,t)=1,p(+\infty ,t)=0.\end{eqnarray} \tag{ 16 }$$

We now fix the temperature and the orientation arbitrarily. Hence, ${M}^{{\rm{PF}}},\sigma$ and also the velocity ${\nu }^{{\rm{M}}{\rm{D}}}$ are constants. The phase field is a traveling wave solution and can be expressed as

$$\begin{eqnarray*}&&p(x,t)=\phi (x-{{\nu }}^{{\rm{M}}{\rm{D}}}{t})=: \phi (\xi ).\end{eqnarray*}$$

Substitution of ϕ in (16) leads to the boundary value problem:

$$\begin{eqnarray}&&0={M}^{{\rm{PF}}}{\sigma }^{2}\phi ^{\prime\prime} -{M}^{{\rm{PF}}}F^{\prime} +{\nu }^{{\rm{M}}{\rm{D}}}\phi ^{\prime} ,\phi (-\infty )=1,\phi (\infty )=0,\end{eqnarray} \tag{ 17 }$$

where $\phi ^{\prime} =\tfrac{{\rm{d}}\phi }{{\rm{d}}\xi }$ and $F^{\prime} =\tfrac{\partial F}{\partial \phi }$.

To receive an initial-value problem, we integrate the ordinary differential equation in (17) with the left boundary $\phi (-\infty )=1$, which holds as our first initial condition, and the right boundary s, which represents the new space-dependence variable for the shooting method. Furthermore, we define $\tilde{F}^{\prime} (\xi ):= F^{\prime} (\phi (\xi ))$, thereby we are able to integrate the derivative of the potential F in ξ. For the calculation of the integral we use the second initial conditions $\phi ^{\prime} (-\infty )=0$. The both initial conditions guarantee, that we have crystal material at the left boundary. Applying all described conditions, we receive

$$\begin{eqnarray}&&{c}_{1}={M}^{{\rm{PF}}}{\sigma }^{2}{\int }_{-\infty }^{s}\phi ^{\prime\prime} \,{\rm{d}}\xi -{M}^{{\rm{PF}}}{\int }_{-\infty }^{s}\tilde{F}^{\prime} (\xi )\,{\rm{d}}\xi +{\nu }^{{\rm{M}}{\rm{D}}}{\int }_{-\infty }^{s}\phi ^{\prime} \,{\rm{d}}\xi ,\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{c}_{2}={M}^{{\rm{PF}}}{\sigma }^{2}(\phi ^{\prime} (s)-0)-{M}^{{\rm{PF}}}(\tilde{F}(s)-\tilde{F}(-\infty ))+{{\nu }^{{\rm{M}}{\rm{D}}}}_{}(\phi (s)-1).\end{eqnarray} \tag{ 19 }$$

By choosing ${c}_{2}={M}^{{\rm{PF}}}\tilde{F}(-\infty )$, our initial value problem has the form:

$$\begin{eqnarray}&&\left(\begin{array}{c}\tilde{F}^{\prime} (s)\\ \phi ^{\prime} (s)\end{array}\right)=\left(\begin{array}{c}F^{\prime} (\phi (s))\\ \displaystyle \frac{{\nu }^{{\rm{M}}{\rm{D}}}}{{M}^{{\rm{PF}}}{\sigma }^{2}}(1-\phi (s))+\displaystyle \frac{1}{{\sigma }^{2}}\tilde{F}(s)\end{array}\right),\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&\phi (-\infty )=1,\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&\phi ^{\prime} (-\infty )=0.\end{eqnarray} \tag{ 22 }$$

We calculate the mobility ${M}^{{\rm{PF}}}$ via a bisection method. Therefore we vary ${M}^{{\rm{PF}}}$ until the right boundary condition $\phi (\infty )=0$ is fulfilled for the numerical solution of (20)–(22). Thereby we solve the initial value problem with Runge-Kutta 4/5. We apply this method for each of the three orientations {100}, {110} and {111} and use the respective value of ${\sigma }_{\{{hkl}\}}$ listed in table 3. Thereby we calculate ${\nu }^{{\rm{M}}{\rm{D}}}$ for the respective temperature and orientation with the fits shown in figure 4. Our results are shown in figure 7 and listed in appendix A.3.

Zoom In Zoom Out Reset image size

We solve model (10) numerically at fixed temperatures $T=1050,1100,\,\ldots ,\,2000\,{\rm{K}}$ for the three crystallographic orientations. During the simulation, we measure the velocity of the interface region by interpolating the position of p = 0.5. In fact, the model reproduces the velocities ${v}_{\{{hkl}\}}^{{\rm{MD}}}(T)$ from molecular dynamics, see figure 8.

Zoom In Zoom Out Reset image size

As for the molecular dynamical simulations, we also have periodic boundary conditions for the phase-field model for $p(x,t)$, which we solve numerically using a Fourier spectral method. Our equidistant grid guarantees that enough grid points are located on the interface to secure an accurate solution. For velocities close to zero, we lower the time step for a better result.

As initial condition we simply define a jump function J(x): for $T\lt {T}_{{\rm{m}}}^{{\rm{V}}}$ we set J(x) equal to one close to the boundaries and zero in between. For $T\gt {T}_{{\rm{m}}}^{{\rm{V}}}$ we define J(x) the other way around. In both cases, one has to take care that the intervals where $J(x)=1$ or $J(x)=0$ are wider than the interface thickness, else the system evolves to an equilibrium state before the traveling wave is established.

Our results of the numerical velocity measurement match accurately with the results from molecular dynamics. On the left-hand side of figure 8, the velocities from the phase-field simulation of model (10) are located directly on the line of the Vogel–Fulcher fit of the molecular dynamical data. For a better comparison, we calculate the relative error, which is

$$\begin{eqnarray}&&{R}_{\{{hkl}\}}(T)=\left|\displaystyle \frac{{v}_{\{{hkl}\}}^{{\rm{MD}}}(T)-{v}_{\{{hkl}\}}^{{\rm{PF}}}(T)}{{v}_{\{{hkl}\}}^{{\rm{MD}}}(T)}\right|.\end{eqnarray} \tag{ 23 }$$

The relative error is shown on the right-hand side of figure 8. We observe, that the maximal relative error is 0.0013, which is at $T=1650\,{\rm{K}}$ for orientation {110}. 95% of the errors are even lower than 10⁻³.

Our approach for the free energy density is a polynomial of fourth degree in p

$$\begin{eqnarray}&&F(p,T)={a}_{0}(T)+{a}_{1}(T)p+{a}_{2}(T){p}^{2}+{a}_{3}(T){p}^{3}+{a}_{4}(T){p}^{4},\end{eqnarray} \tag{ A.1 }$$

which attains the values in figure 2 at the equilibrium states p = 0 (liquid phase) and p = 1 (crystalline phase). So we need to consider the derivatives of F with respect to p:

$$\begin{eqnarray*}\displaystyle \frac{\partial F}{\partial p}(p,T) & = & {a}_{1}(T)+2{a}_{2}(T)p+3{a}_{3}(T){p}^{2}+4{a}_{4}(T){p}^{3},\\ \displaystyle \frac{{\partial }^{2}F}{\partial {p}^{2}}(p,T) & = & 2{a}_{2}(T)+6{a}_{3}(T)p+12{a}_{4}(T){p}^{2}.\end{eqnarray*}$$

From the conditions for the existence of minima at p = 0 and p = 1 follows

$$\begin{eqnarray}&&{a}_{1}(T)\equiv 0,\end{eqnarray} \tag{ A.2 }$$

$$\begin{eqnarray}&&{a}_{2}(T)\gt 0,\end{eqnarray} \tag{ A.3 }$$

$$\begin{eqnarray}&&{a}_{3}(T)=-\displaystyle \frac{2}{3}{a}_{2}(T)-\displaystyle \frac{4}{3}{a}_{4}(T),\end{eqnarray} \tag{ A.4 }$$

$$\begin{eqnarray}&&{a}_{2}(T)\lt 2{a}_{4}(T).\end{eqnarray} \tag{ A.5 }$$

The polynomials ${f}_{0}(T)$ and ${f}_{1}(T)$ pass the equilibrium values in figure 2 for liquid/amorphous and crystalline $\mathrm{Si}$, respectively. For the minima, we need for (A.1) the equalities $F(0,T)={f}_{0}(T)$ and $F(1,T)={f}_{1}(T)$ and hence

$$\begin{eqnarray}&&{f}_{0}(T)=F(0,T)={a}_{0}(T),\end{eqnarray} \tag{ A.6 }$$

$$\begin{eqnarray}&&{f}_{1}(T)=F(1,T)\mathop{=}\limits^{({\rm{A}}.2,{\rm{A}}.4,{\rm{A}}.6)}{f}_{0}(T)+\displaystyle \frac{1}{3}{a}_{2}(T)-\displaystyle \frac{1}{3}{a}_{4}(T).\end{eqnarray} \tag{ A.7 }$$

With (A.4) and (A.7) the two coefficients a₃ and a₄ have the form

$$\begin{eqnarray}&&{a}_{3}(T)=-2{a}_{2}(T)-4({f}_{0}(T)-{f}_{1}(T))=-2{a}_{2}(T)-4H(T),\end{eqnarray} \tag{ A.8 }$$

$$\begin{eqnarray}&&{a}_{4}(T)={a}_{2}(T)+3({f}_{0}(T)-{f}_{1}(T))={a}_{2}(T)+3H(T).\end{eqnarray} \tag{ A.9 }$$

Together with (A.9) we can verify the inequality (A.5):

$$\begin{eqnarray}&&{a}_{2}(T)\gt -6H(T),\end{eqnarray} \tag{ A.10 }$$

where the right-hand side of (A.10) is negative if and only if $T\lt {T}_{{\rm{m}}}^{{\rm{V}}}$, as one can observe in figure 2. Hence, for temperatures below the melting point, a₂ has to fulfill (A.3). Finally, together with (A.8), (A.9), the double-well potential (A.1) has the form

$$\begin{eqnarray}&&F(p,T)={f}_{0}(T)+{a}_{2}(T){p}^{2}-2({a}_{2}(T)+2H(T)){p}^{3}+({a}_{2}(T)+3H(T)){p}^{4},\end{eqnarray} \tag{ A.11 }$$

where

$$\begin{eqnarray}{a}_{2}(T)\gt \left\{\begin{array}{ll}0 & T\leqslant {T}_{{\rm{m}}}^{{\rm{V}}}\\ -6{H}(T) & T\gt {T}_{{\rm{m}}}^{{\rm{V}}}.\end{array}\right.\end{eqnarray} \tag{ A.12 }$$

The formulations (A.11) and (A.12) are equivalent to (11) and (12), respectively, by renaming ${a}_{2}(T)$ to a(T).

In the following, the calculation of μ is described, where the notation of (11) and (12) is used. We first note that with (12) we have

$$\begin{eqnarray}&&a(T)+3H(T)\gt 0,\end{eqnarray} \tag{ A.13 }$$

which can be seen easily by using the restrictions on a in (12):

First of all we note that for $T\gt {T}_{{\rm{m}}}^{{\rm{V}}}$

$$\begin{eqnarray*}&&a(T)+3H(T)\mathop{\gt }\limits^{(12)}-6H(T)+3H(T)=-3H(T)\mathop{\gt }\limits^{T\gt {T}_{{\rm{m}}}^{{\rm{V}}}}0\end{eqnarray*}$$

and for $T\leqslant {T}_{{\rm{m}}}^{{\rm{V}}}$,

$$\begin{eqnarray*}&&a(T)+3H(T)\gt 0+3H(T)\geqslant 0.\end{eqnarray*}$$

Besides that p = 0 and p = 1, the first derivative of F

$$\begin{eqnarray*}&&\displaystyle \frac{\partial F}{\partial p}(p,T)=2p(1-p)[a(T)-(2a(T)+6H(T))p]\end{eqnarray*}$$

has a third root $\mu (T)$

$$\begin{eqnarray*}&&\mu (T)=\displaystyle \frac{a(T)}{2a(T)+6H(T)}.\end{eqnarray*}$$

Since (A.13) holds, the denominator of $\mu (T)$ can not become zero. Hence, with (12) we have $\mu (T)\gt 0$. Further, the condition

$$\begin{eqnarray}&&\mu (T)=\displaystyle \frac{a(T)}{2a(T)+6H(T)}\lt 1\end{eqnarray} \tag{ A.14 }$$

is equivalent to

$$\begin{eqnarray}&&a(T)\gt -6H(T)\end{eqnarray} \tag{ A.15 }$$

and together with (12) we have $\mu (T)\in (0,1)$ for all $T\geqslant 0$. In addition, for a maximum in $\mu (T)$, the second derivative of F in $\mu (T)$ has to be less than zero, i.e.

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}F}{\partial {p}^{2}}(\mu (T),T)=\displaystyle \frac{a(T)(a(T)+6H(T))}{-a(T)-3H(T)}\mathop{\lt }\limits^{!}0.\end{eqnarray} \tag{ A.16 }$$

Since (A.13) holds, the denominator of (A.16) is negative and we get, since $a(T)\gt 0$

$$\begin{eqnarray}&&a(T)\gt -6H(T).\end{eqnarray} \tag{ A.17 }$$

With (12), the last condition is a always true. Hence, the second derivative of F in $\mu (T)$ is less than zero, so $(\mu (T),F(\mu (T),T))$ is a maximum point of F.

In section 4.3 we describe the calculation of the mobility parameters for the phase-field model (10). The mobility values are shown in figure 7 and listed in table A1.

Table A1.  The calculated mobilities via a shooting method for orientations {100} {110} and {111}. They are shown in figure 7.

Temperature	Mobility ${M}_{\{100\}}^{{\rm{PF}}}$	Mobility ${M}_{\{110\}}^{{\rm{PF}}}$	Mobility ${M}_{\{111\}}^{{\rm{PF}}}$
T in K	in Å3 (eV ns)−1	in Å3 (eV ns)−1	in Å3 (eV ns)−1
$\lt 800$	0	0	0
800	8.8818 × 10−14	8.8818 × 10−14	2.93 × 10−1
850	9.7105 × 10−8	1.7795 × 10−9	5.6933 × 10−1
900	5.8741 × 10−3	8.4051 × 10−4	1.0082
950	2.2923 × 10−1	6.8938 × 10−2	1.658
1000	1.425	6.2816 × 10−1	2.5636
1050	3.3476	1.8627	3.761
1100	6.8805	4.4784	5.2861
1150	11.4558	8.3461	7.1611
1200	16.737	13.2755	9.4112
1250	22.4071	18.9898	12.0447
1300	28.2093	25.213	15.06
1350	33.9597	31.7057	18.4506
1400	39.5291	38.274	22.2025
1450	44.8397	44.7766	26.3002
1500	49.7514	51.0148	30.6662
1550	54.1681	56.8348	35.2213
1600	58.1903	62.2964	39.9958
1650	61.8892	67.4445	45.0044
1700	63.6398	70.4575	48.9622
1750	66.429	74.5941	54.0443
1800	68.8907	78.3509	59.2121
1850	70.7543	81.4046	64.1717
1900	71.9325	83.6335	68.749
1950	72.9413	85.6235	73.3545
2000	73.6987	87.2768	77.8687