Molecular Dynamics Study of the Effect of Carbon Atoms on the Surface Tension of Silicon–carbon Alloy

Abstract

We conducted molecular dynamics simulations of Si–C alloy to understand the atomistic behavior of solute C atoms near the melt surface and to estimate the surface tension. The surface tensions of liquid Si and C were first evaluated and compared with experimental values and those for other metals. The composition dependence of the surface tension of Si–C alloy was then evaluated, and compared with estimates obtained using the modified Butler’s model. The behavior of C atoms at the surface of liquid Si–C alloys is also discussed.

1. Introduction

High-temperature solution growth using silicon-based liquid alloys as solvents is a promising method for producing high quality SiC single crystals.¹⁾ The solution growth rate was improved by using alloy solvents with high carbon solubilities such as Si–Cr,²⁾ Si–Ti³⁾ and Fe–Si⁴⁾ alloys. These have reportedly accounted for growth rates of 100‒300 μm/h which are comparable to those obtained using the conventional physical vapor transport (PVT) method. Crystals grown by the solution growth method reportedly also have smaller dislocation densities⁵⁾ than those grown by PVT. Diodes fabricated from such crystals exhibit the equivalent or greater withstand voltages than those fabricated from crystals grown by PVT.⁶⁾

To improve the process by growing larger diameter crystals or by increasing the growth rate, precise control of the temperature and flow in the solution is essential. Computational fluid dynamics (CFD) simulations enable the estimation of temperature and flow distributions in the solution, and are frequently employed for the process design.^7,8,9) However, data on the melt properties of such solvents is limited, so the properties of pure Si melt are often used instead in CFD simulations. Understanding the melt properties of Si-based alloys containing carbon is important for precisely estimating the growth conditions in various solvent systems. Especially during growth without applying forced convection in solution, the remarkable growth of SiC was found at the edge of the crystal/solution contact area with cylindrical liquid bridges.¹⁰⁾ This is presumably due to the Marangoni effect dominating the convection around the meniscus of the liquid bridge. The surface tension of the alloy solution and its temperature and composition dependence should therefore be known for controlling the Marangoni flow in the solution.

Studies on the surface tension of alloys containing carbon with the tendency to form stable metal carbides have shown varying results. Belton¹¹⁾ measured the surface tension of Fe–Si–C and Fe–Cr–C alloys using a sessile drop method at a constant C content, and found a decrease in surface tension with increasing Si and Cr contents. This suggested a surfactive effect of the associative adsorption of SiC and CrC. Kawai et al.¹²⁾ also used the sessile drop method to measure the surface tension of Fe–(1‒10 mol%)Si–(1‒3 mol%)C alloys, and found that the maximum surface tension occurred at 2 mol% Si at constant C content. To elucidate the surfactive effect of associative SiC at the melt surface, one of the authors¹³⁾ measured the surface tension of Fe–(30‒40 mol%)Si alloys using the maximum bubble pressure method. The author concluded insignificant surface segregation of C or associative carbide at the surface, because of little observed difference in the surface tension with and without the addition of carbon into the Fe–Si alloy. Recently, Eustathopoulos et al.¹⁴⁾ reviewed the surface tension of Si, and pointed out that C in the Si melt had no effect on the surface tension within the maximum carbon solubility. However, the behavior of C atoms at the melt surface has not yet been confirmed.

Molecular dynamics (MD) simulations are a useful method for investigating the behavior of atoms. Recent advances in MD have allowed the estimation of surface tension^15,16) and other physical properties.^17,18) In the current study, we conducted MD simulations of stable and metastable Si–C melts with free surfaces, to analyze the microscopic behavior of Si and C atoms at the surface and to estimate the surface tension. The surface tension of the Si–C melts was also studied using the modified Butler’s model, and the effect of carbon on the melt surface of Si is discussed.

2. Simulation Method

A classical MD simulation was performed using LAMMPS.¹⁹⁾ The simulation time steps were set to 0.4 fs for Si and 0.1 fs for C and the Si–C alloy. A Nosé-Hoover thermostat and barostat were employed to control the temperature and pressure, respectively.^20,21) A periodic boundary was applied to all boundaries of the simulation cell. OVITO²²⁾ was used to perform post-processing for atom visualization. According to our previous study,²³⁾ Tersoff_94 potential²⁴⁾ was used for interatomic interactions of the Si–C system. A total of 15552 atoms (cell volume: 18 × 6 × 18 of the unit cell) of Si and C were placed in a cell with a diamond structure. The cell of the Si and Si–C alloy was heated to 2600‒4000 K and the cell of C was heated to 6700‒7100 K, and then kept at 1 atm for 100 ps with the isobaric-isothermal ensemble to make a uniform liquid phase. Here it should be noted that the melting point of Si and diamond derived from Tersoff_94 potential²⁴⁾ were found to be 2560 K and 6600 K, respectively, in a previous study.²³⁾ A cell of the liquid phase possessing free surfaces (vapor/liquid cell) was prepared by combining the same size of vacuum in one direction (18 × 18 × 18 of the unit cell), as shown in Fig. 1. To confirm the effect of the surface area on the simulation results, the number of Si atoms was varied from 768 to 27648 (6n × 18× 6n of the unit cell, n = 2/3, 1, 2, 3, 4). Some of Si–C alloy employed in this study was oversaturated in comparison to the carbon solubility in the Si melt.²⁵⁾ We confirmed that Si–(10‒50 mol%)C alloy formed the uniform liquid as a stable or metastable phase for a calculation temperature range by the identify diamond structure (IDS) modifier²⁶⁾ implemented in OVITO. The vapor/liquid cell was heated to various temperatures and kept at constant volume for 500 ps. The cells were sliced in 400 slabs for the Si melt and Si–C alloys and in 280 slabs for the C melt toward the y axis to satisfy the thickness of the single slab smaller than the monoatomic layer thickness. The number of atoms, the number of the nearest neighbors and the potential energy in each slab were extracted to discuss the surface structure and to estimate the surface tension. The average potential energy of the bulk slabs was subtracted from the potential energy of each slab and divided by the surface area to calculate the surface tension, σ, as follows;   

$$\sigma =\frac{1}{A}{\sum }_{i=1}^{400}\mathrm{\Delta }{E}_i$$(1)

  

$$\mathrm{\Delta }{E}_i=E_i-E_{Bulk} if\mathrm{\hspace{0.17em} }{E}_i <0$$

  

$$\mathrm{\Delta }{E}_i=0 if\mathrm{\hspace{0.17em} }{E}_i \ge 0$$

where A is the surface area of the liquid surface determined by constructing surface mesh²²⁾ implemented in OVITO; ΔE_i is the excess energy of the ith slab; E_i is the potential energy of the ith slab; E_Bulk is the average potential energy of the bulk liquid determined from the inflection point of the first derivative of the number density distribution of atoms. Because the melting point of Si, T_{m.p. of Si}, derived from Tersoff_94 potential²⁴⁾ (2560 K) deviated from the experimental value (1687 K), the normalized temperature (T/T_{m.p. of Si}) was used hereafter.

Fig. 1.

Cell of Si melt including free surface (vapor/liquid cell) at T/T_{m.p. of Si} = 1.02. (Online version in color.)

3. Results and Discussion

3.1. Surface Tension of Si and C Melts

Figure 2 shows the time-mean values of the number of atoms, the number of nearest neighbor atoms, the potential energy and the excess energy from Eq. (1) of a Si melt (15552 atoms) at T/T_{m.p. of Si} = 1.02 (2600 K in the simulation) as an example. The duration was taken from the last 80 ps of the simulation results to obtain the time-mean values. As shown in Fig. 2(a), the number of Si atoms changed sharply near the melt surface. The thickness of the melt surface and the bulk was determined from the inflection point of the first derivative of the number density distribution of atoms. The thickness of the melt surface was about 0.25‒0.30 nm, which was comparable to 1 atomic layer because the Si–Si bond length in the liquid phase was about 0.25 nm by using Tersoff_94 potential.²⁴⁾ The excess energy was distributed in both free surfaces, corresponding to the distribution of the number of nearest neighbor atoms, Z. Figures 3(a) and 3(b) show the predicted radial distribution functions (RDF) of Si atoms and an extracted histogram of the Z for the Si melt in the bulk and surface at T/T_{m.p. of Si} = 1.02 (2600 K in the simulation), respectively. The ratio of the mean Z_Bulk to Z_Surface obtained in this study was 0.84, which agreed with the value proposed by Tanaka et al.²⁷⁾ for general liquid metals. Figure 4 shows the correlation of the surface tension of the Si melt with the number of particles from 768 to 27648 at T/T_{m.p. of Si} = 1.02. The dashed lines are experimental values.^{14,28,29,30,31)} The estimated surface tension increased with increasing number of Si atoms, and then became constant, agreeing with the reported values when the number of atoms was larger than 1728. Therefore, we largely conducted simulations with 15552 atoms.

Fig. 2.

Distribution of (a) the particle number and its first derivative, (b) the potential energy, (c) the excess energy and (d) the first coordination number at T/T_{m.p. of Si} = 1.02.

Fig. 3.

(a) Radial distribution functions (RDF) of Si atoms and (b) histograms of the first coordination number of the bulk and surface for the Si melt at T/T_{m.p. of Si} = 1.02. (Online version in color.)

Fig. 4.

Correlation of the number of Si atoms with estimated surface tension at T/T_{m.p. of Si} = 1.02 together with (1) Eustathopoulos et al.,¹⁴⁾ (2) Keene,²⁸⁾ (3) Fujii et al.,²⁹⁾ (4) Millot et al.³⁰⁾ and (5) Zhou et al.³¹⁾

Figure 5 shows the surface tension of the Si melt as a function of temperature along with experimental values,^{14,28,29,30,31)} which are summarized in Table 1. Surfactive-contamination such as oxygen affects the surface tension, so several experimental values obtained by the container-less method^29,30,31) were chosen for comparison. The surface tension obtained in the current study agreed with data reported by Fujii et al.,²⁹⁾ Millot et al.,³⁰⁾ Zhou et al.³¹⁾ and reviewed by Keene.²⁸⁾ Although the estimated values were slightly smaller than the values suggested by Eustathopoulos et al.,¹⁴⁾ the estimated temperature coefficient agreed with the experimental values. The densities of liquid Si and C are listed in Table 2. The derived density of liquid Si was slightly smaller than the experimental value.³²⁾ We previously showed that the derived latent heat of fusion by MD simulation agreed fairly well with the experimental value.²⁴⁾ Figure 6 shows the proportional correlation of the surface tensions of pure metals at their melting points to their heats of evaporation derived under the assumption of Z_Bulk/Z_Surface = 0.83.²⁷⁾ The surface tension of liquid C shows the same tendency although those of liquid Si by the present prediction and of the reported experimental value are smaller than the estimated line in spite of Z_Bulk/Z_Surface = 0.84 as described above. We speculate that this tendency is attributed to the difference of the interatomic bonding energy of Si at the surface from that in the bulk although it was not clarified in the present work.

Fig. 5.

Estimated surface tension of the Si melt and comparison with experimental values from (1) Eustathopoulos et al.,¹⁴⁾ (2) Keene,²⁸⁾ (3) Fujii et al.,²⁹⁾ (4) Millot et al.³⁰⁾ and (5) Zhou et al.³¹⁾

Table 1. Estimated surface tension of the Si melt and C melt and comparison with experimental values.^{14,28,29,30,31)}

element	Surface tension/mN·m−1	Literature
Si	742 - 0.0776 (T - 2560)	This study (15552 atoms)
	840 - 0.190 (T - 1685)	Eustathopoulos et al.14)
	733 - 0.062 (T - 1687)	Keene28)
	732 - 0.086 (T - 1685)	Fujii et al.29)
	775 - 0.145 (T - 1683)	Millot et al.30)
	721 - 0.0615 (T - 1687)	Zhou et al.31)
C	3910 - 0.582 (T - 6600)	This study (15552 atoms)

Table 2. Estimated liquid density of the Si melt and C melt and comparison with experimental values.^32,33)

element	This study	Density/g·cm−3 in literature
Si	2.31 - 1.02 × 10−4 (T - 2560)	2.51 - 2.72 × 10−4 (T - 1687)32)
C	2.84 - 4.18 × 10−4 (T - 6600)	3.2 - 4.5 at melting point of graphite33)

Fig. 6.

Correlation of the surface tension of various metals with their latent heat of evaporation. Open circles are estimated values from the current study and black dots are values reported by Tanaka et al.²⁷⁾

3.2. Surface Tension of Si–C Binary System

Figure 7 shows the time-mean values of the number of Si and C atoms for the Si–10mol%C alloy at T/T_{m.p. of Si} = 1.02. It was found that C atoms did not segregate at the melt surface. The predicted radial distribution functions (RDF) of global, Si–Si, Si–C, and C–C atoms of the bulk and surface for the Si–10mol%C alloy at T/T_{m.p. of Si} = 1.02 are shown in Fig. 8. The consecutive tetrahedral structure generally seen in solid Si or in SiC and any clusters did not appear in the analysis using IDS modifier.²⁶⁾ Hence, it is reasonable to consider that Si atoms and C atoms placed randomly in the liquid of the MD cell. The surface tension of the Si–10mol%C alloy is shown in Fig. 9 together with the evaluated values for liquid Si in Fig. 5 and the C melt. The T/T_{m.p. of Si} was used for the Si melt and Si–C alloy, whereas the T/T_{m.p. of diamond} was used for liquid C. The surface tension of Si–10mol%C alloy was close to that of liquid Si. The temperature coefficient of the Si–10mol%C alloy was also comparable to that of liquid Si. The surface tension of liquid C was five times larger than those of the liquid Si and Si–10mol%C alloy. This implied that it was energetically unfavorable when C atoms existed at the melt surface.

Fig. 7.

Distribution of the particle number of Si atoms and C atoms at T/T_{m.p. of Si} = 1.02.

Fig. 8.

Radial distribution functions (RDF) of global, Si–Si, Si–C and C–C atoms of the bulk and surface for the Si–10mol%C alloy at T/T_{m.p. of Si} = 1.02. (Online version in color.)

Fig. 9.

Estimated surface tensions of the Si melt, C melt and Si–10mol%C alloy.

Figure 10 shows the estimated surface tensions of Si‒10‒50mol%C alloys (15552 atoms) compared with those of the Si melt and C melt at T/T_{m.p. of Si} = 1.56 (4000 K in the simulation). Here, the extrapolated value was used because the melting point of diamond is higher than 4000 K. Since the melting point of SiC was determined as 3850 K from a previous study,²⁴⁾ we chose 4000 K for the comparison of the surface tension of Si–(10‒50 mol%)C. Although the surface tension increased with increasing C concentration, it was smaller than the interpolated value between the surface tension of Si and C.

Fig. 10.

(a) Correlation of the estimated surface tension of the Si–C alloy with C concentration at T/T_{m.p. of Si} = 1.56 together with the estimation obtained using the modified Butler’s model. (b) Correlation of the C concentration at the melt surface with C concentration in the bulk liquid of the Si–C alloy.

Here, the estimated surface tensions of Si–C alloys were compared with estimates obtained using the modified Butler’s model,^27,34,35) where the surface tension of the Si–C binary alloy is represented by the following equations;   

$$\begin{array}{l}\sigma ={\sigma }_{Si}^{Pure}+{\displaystyle \frac{RT}{A_{Si}^{mol}}}ln{\displaystyle \frac{1-X_C^{Surface}}{1-X_C^{Bulk}}}\\ \hspace{1em}+{\displaystyle \frac{1}{A_{Si}^{mol}}}\{{\overline{G}}_{Si}^{Ex,\mathrm{\hspace{0.17em} }Surface}(T,\mathrm{\hspace{0.17em} }{X}_C^{Surface})-{\overline{G}}_{Si}^{Ex,\mathrm{\hspace{0.17em} }Bulk}(T,\mathrm{\hspace{0.17em} }{X}_C^{Bulk})\}\end{array}$$(2)

  

$$\begin{array}{l}\sigma ={\sigma }_C^{Pure}+{\displaystyle \frac{RT}{A_C^{mol}}}ln{\displaystyle \frac{X_C^{Surface}}{X_C^{Bulk}}}\\ \hspace{1em}+{\displaystyle \frac{1}{A_C^{mol}}}\{{\overline{G}}_C^{Ex,\mathrm{\hspace{0.17em} }Surface}(T,\mathrm{\hspace{0.17em} }{X}_C^{Surface})-{\overline{G}}_C^{Ex,\mathrm{\hspace{0.17em} }Bulk}(T,\mathrm{\hspace{0.17em} }{X}_C^{Bulk})\}\end{array}$$(3)

where ${\sigma }_i^{Pure}$ and $A_i^{mol}$ are the surface tension and molar surface area of pure i, respectively; $X_i^{Bulk}$ and $X_i^{Surface}$ are the mole fractions of i in the bulk and surface of the liquid, respectively; ${\overline{G}}_i^{Ex,\mathrm{\hspace{0.17em} }Bulk}$ and ${\overline{G}}_i^{Ex,\mathrm{\hspace{0.17em} }Surface}$ are the excess partial molar Gibbs energies of i in the bulk and surface of the liquid, respectively. The molar surface area is given using the molar volume of i, $V_i^{mol}$, and Avogadro’s constant, N_A, in Eq. (4).   

$$A_i^{mol}=1.091N_A^{1/3}{(V_i^{mol})}^{2/3}$$(4)

The excess partial molar Gibbs energy in the surface monolayer was assumed as in Eq. (5)²⁷⁾   

$$\begin{array}{l}{\overline{G}}_C^{Ex,\mathrm{\hspace{0.17em} }Surface}(T,\mathrm{\hspace{0.17em} }{X}_C^{Surface})={\displaystyle \frac{Z^{Surface}}{Z^{Bulk}}}{\overline{G}}_C^{Ex,\mathrm{\hspace{0.17em} }Bulk}(T,\mathrm{\hspace{0.17em} }{X}_C^{Surface})\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.5em} =0.83{\overline{G}}_C^{Ex,\mathrm{\hspace{0.17em} }Bulk}(T,\mathrm{\hspace{0.17em} }{X}_C^{Surface})\end{array}$$(5)

The estimation using the modified Butler’s equation was then performed for two different cases. First, the estimation was performed for the system where both the bulk and surface phases were considered to be a mixture of Si and C. The surface tension and the surface C concentration of the Si–C alloy were estimated from Eqs. (2) and (3), and are shown as solid curves in Fig. 10. Second, the surface of the Si–C alloy was assumed to be composed of only Si atoms as assumed in ref. 13, and the surface tension was estimated from Eq. (2) under the condition of $X_c^{Surface}$ = 0. This is shown as a dashed curve in Fig. 10. The difference between the two estimated curves was insignificant. As shown in Fig. 10(b), the surface C concentration was one order of magnitude smaller than the bulk C concentration in the composition region of the Si–C alloy. A discrepancy between the surface tension of the Si–C alloy obtained by MD simulation and the estimated values obtained by the modified Butler’s model occurred when the C concentration of the Si–C alloy increased. Although Tersoff_94 potential was the most suitable among reported Tersoff potentials in a previous study,²³⁾ the heat of mixing of the Si–C solution was different from that in ref.²⁵⁾ This indicated that Tersoff_94 potential was still insufficient to fully represent the Si–C interaction. Therefore, the discrepancy became significant at large C concentrations, while there was minimal impact on the simulation results at small C concentrations. When simulation results at large C concentrations are required, the Si–C interatomic potential will need to be revised.

In addition to the insignificant concentration of carbon at the surface of the Si–C alloy by MD simulation shown in Fig. 7, the modified Butler’s model also negated the surface segregation of C in the alloy. Hence, we determined that carbon does not behave as a surfactive component at the surface of the Si–C alloy.

4. Conclusions

The effect of C atoms on the surface tension of Si–C alloy was investigated by MD simulations. The following conclusions were made:

(1) The derived surface tension and temperature coefficient of Si in this study agreed well with experimental values.

(2) The surface tension of Si‒10mol%C alloy was close to that of liquid Si. From the number densities of Si and C atoms, it was found that C atoms did not segregate at the melt surface.

(3) The composition dependence of the derived surface tension of the Si–C alloy was compared with estimates obtained using the modified Butler’s model. The surface tension by MD simulation was consistent with the model at low C concentrations. The surface C concentration derived by MD simulation was estimated to be one order of magnitude smaller than the bulk C concentration. Hence, carbon did not behave as a surfactive component.

Acknowledgements

Part of this research was supported by a Grant-in-Aid for Scientific Research (B) (Grant Number 15H04166) and a Grant-in-Aid for JSPS Fellows (Grant Number 16J10300) from MEXT, Japan. We thank Aidan G. Young, PhD, from Edanz Group (www.edanzediting.com/ac) for editing a draft of this manuscript.

References

• 1)   D. H.  Hofmann and  M. H.  Müller: Mater. Sci. Eng. B, 61–62 (1999), 29.
• 2)   M.  Kado,  H.  Daikoku,  H.  Sakamoto,  H.  Suzuki,  T.  Bessho,  N.  Yashiro,  K.  Kusunoki,  N.  Okada,  K.  Moriguchi and  K.  Kamei: Mater. Sci. Forum, 740–742 (2013), 73.
• 3)   K.  Kusunoki,  K.  Kamei,  N.  Okada,  N.  Yashiro,  A.  Yauchi,  T.  Ujihara and  K.  Nakajima: Mater. Sci. Forum, 527–529 (2006), 119.
• 4)   T.  Yoshikawa,  S.  Kawanishi and  T.  Tanaka: Jpn. J. Appl. Phys., 49 (2010), 051302.
• 5)   K.  Kamei,  K.  Kusunoki,  N.  Yashiro,  N.  Okada,  K.  Moriguchi,  H.  Daikoku,  M.  Kado,  H.  Suzuki,  H.  Sakamoto and  T.  Bessho: Mater. Sci. Forum, 717–720 (2012), 45.
• 6)   K.  Kusunoki,  K.  Seki,  H.  Daikoku,  H.  Saito,  I.  Kobayashi and  H.  Mihara: Abstract of European Conf. on Silicon Carbide and Related Materials, European Conf. on Silicon Carbide and Related Materials 2018, Birmingham, (2018), INV.P.02.
• 7)   F.  Mercier,  J. M.  Dedulle,  D.  Chaussende and  M.  Pons: J. Cryst. Growth, 213 (2010), 155.
• 8)   T.  Yamamoto,  N.  Adkar,  Y.  Okano,  T.  Ujihara and  S.  Dost: J. Cryst. Growth, 474 (2017), 50.
• 9)   K.  Ariyawong,  Y. J.  Shin,  J. M.  Dedulle and  D.  Chaussende: Cryst. Growth Des., 16 (2016), 3231.
• 10)   S.  Kawanishi,  T.  Yoshikawa,  K.  Morita,  N.  Okada,  K.  Kusunoki and  K.  Kamei: J. Cryst. Growth, 381 (2013), 121.
• 11)   G.  Belton: Metall. Trans., 3 (1972), 1465.
• 12)   Y.  Kawai,  K.  Mori,  M.  Kishimoto,  K.  Ishikura and  T.  Shimada: Tetsu-to-Hagané, 60 (1974), 29 (in Japanese).
• 13)   T.  Yoshikawa: Mater. Trans., 54 (2013), 1968.
• 14)   N.  Eustathopoulos and  B.  Drevet: J. Cryst. Growth, 371 (2013), 77.
• 15)   M. J. P.  Nijmeijer,  A. F.  Bakker,  C.  Bruin and  J. H.  Sikkenk: J. Chem. Phys., 89 (1988), 3789.
• 16)   H.  Daiguji and  E.  Hihara: Heat Mass Transf., 35 (1999), 213.
• 17)   F.  Müller-Plathe: J. Chem. Phys., 106 (1998), 6082.
• 18)   F.  Müller-Plathe: Phys. Rev. E, 59 (1999), 4894.
• 19)   S.  Plimpton: J. Comput. Phys., 117 (1995), 1.
• 20)   S.  Nosé: J. Chem. Phys., 81 (1984), 511.
• 21)   W. G.  Hoover: Phys. Rev. A, 31 (1985), 1695.
• 22)   A.  Stukowski: Modell. Simul. Mater. Sci. Eng., 18 (2010), 015012.
• 23)   T.  Narumi,  Y.  Shibuta and  T.  Yoshikawa: J. Cryst. Growth, 494 (2018), 36.
• 24)   J.  Tersoff: Phys. Rev. B, 49 (1994), 16349.
• 25)   S.  Kawanishi,  T.  Yoshikawa and  T.  Tanaka: Mater. Trans., 50 (2009), 806.
• 26)   E.  Maras,  O.  Trushin,  A.  Stukowski,  T.  Ala-Nissila and  H.  Jonsson: Comput. Phys. Commun., 205 (2016), 13.
• 27)   T.  Tanaka,  K.  Hack,  T.  Iida and  S.  Hara: Z. Metallkd., 87 (1996), 380.
• 28)   B. J.  Keene: Int. Mater. Rev., 38 (1993), 157.
• 29)   H.  Fujii,  T.  Matsumoto,  S.  Izutani,  S.  Kiguchi and  K.  Nogi: Acta Mater., 54 (2006), 1221.
• 30)   F.  Millot,  V.  Sarou-Kanian,  J. C.  Rifflet and  B.  Vinet: Mater. Sci. Eng. A, 495 (2008), 8.
• 31)   Z.  Zhou,  S.  Mukherjee and  W. K.  Rhim: J. Cryst. Growth, 257 (2003), 350.
• 32)   Y.  Inatomi,  F.  Onishi,  K.  Nagashio and  K.  Kuribayashi: Int. J. Thermophys., 28 (2007), 44.
• 33)   A. I.  Savvatimskiy: Carbon, 43 (2005), 1115.
• 34)   J. A. V.  Butler: Proc. R. Soc. A, 135 (1932), 348.
• 35)   K. S.  Yeum,  R.  Speiser and  D. R.  Poirier: Metall. Trans. B, 20 (1989), 693.