Phase diagram of grain boundary facet and line junctions in silicon

The presence of facets and line junctions connecting facets on grain boundaries (GBs) has a strong impact on the properties of structural, functional, and optoelectronic materials: They govern the mobility of interfaces, the segregation of impurities, as well the electronic properties. In the present paper, we employ density-functional theory and modiﬁed embedded atom method calculations to systematically investigate the energetics and thermodynamic stability of these defects. As a prototype system, we consider (cid:2) 3 tilt GBs in Si. By analyzing the energetics of different faceted GBs, we derive a diagram that describes and predicts the reconstruction of these extended defects as a function of facet length and boundary inclination angle. The phase diagram sheds light upon the fundamental mechanisms causing GB faceting phenomena. It demonstrates that the properties of faceting are not determined solely by anisotropic GB energies but by a complex interplay between geometry and microstructure, boundary energies as well as long-range strain interactions.


I. INTRODUCTION
Crystal defects such as grain boundaries (GBs) severely impact materials properties [1]. For structural materials, the Hall-Petch strengthening of metals [2,3] and the brittle fracture induced by preferential segregation of impurities at GBs [4] are characteristic examples of how GBs affect mechanical properties. The role of GBs is also substantial in multi-and polycrystalline optoelectronic materials. Today, multicrystalline Si (mc-Si) dominates the Si-based photovoltaics industry thanks to its lower production cost than single-crystal Si. However, mc-Si has a density of 10 4 cm −1 GBs [5] and these defects as well as their interaction with impurities constitute one major limiting factor in the efficiency of the aforementioned devices. A prominent example is the recombination of electrons at GBs. This mechanism constitutes the major electrical losses channel in high-performance mc-Si solar cells [6]. Indeed, GBs often act as strong gettering centers for metal impurities due to the presence of over-or undercoordinated atoms and the different-than-bulk strain. These impurities, the over-or undercoordinated and/or highly strained host atoms may introduce deep states into the fundamental band gap. Such states would act as recombination centers and reduce the efficiency of solar cells [7].
The 3 tilt GBs with 110 rotation axis in Si constitute a system of GBs with special fundamental and technological interest: 3 GBs constitute up to 80% of GBs in Si * lymperakis@mpie.de ingots grown by dendritic casting [8]. These boundaries are commonly observed to facet toward the low index {111} and {112} boundaries. The {111} and {112} boundaries are also formed by the dissociation of higher value GBs [9]. The 3{111} twin boundaries are coherent and can be viewed as hexagonal inclusions into the cubic lattice. These boundaries have negligible boundary energy, i.e., less than 0.5 meV/Å 2 and are electrically inactive. Hence they are preferred over other electrically active boundaries. On the other hand, the 3 {112} GBs have two stable/metastable states, the symmetric (S-3) and the asymmetric (A-3), with the latter having lower energy. The two states of the 3 {112} GB have different electronic properties. They also show potential as impurity gettering centers: The presence of a fivefold coordinated atom and highly strained bonds at the S-3 boundary introduces deep states into the fundamental band gap while the reconstructed A-3 GB has no deep states in the gap [10]. Furthermore, density-functional theory (DFT) calculations show the S-3 GB to favor interstitial Fe [11] as well as P and As substitutional [10] segregation. Both {112} boundaries also show very different site selectivity for C substitutionals [12].
A striking difference between the A-3 and S-3 GBs is that the latter is commensurate to the underlying coincidence site lattice (CSL) while the former requires removal of atoms from the interface area as well as a rigid shift of the one grain with respect to the other. Recently, a unique anisotropic segregation mechanism going beyond the classical planar McLean-type segregation has been demonstrated at faceted Si GBs [13]. Specifically, by combining DFT with modified embedded atom method (MEAM) potential calculations, it was shown that differences in the local atomic geometry at the 3 {111} and A-3{112} line junctions associated with the strain arising from the above-mentioned partial dislocations are the origin of the aforementioned preferential impurity segregation. However, the driving force that causes faceting in these system is still an open issue.
Following Herring's thermodynamic arguments, faceting of interfaces is primarily driven by a high anisotropy in the boundary energies and results from the minimization of the interfacial free energy [14]. The facet length depends not only on the energy of the facets but on the interaction energy between dislocations and line forces as well [15,16] (see Fig. 1). Faceting may also be caused by impurities and/or changes in composition [17][18][19]: It has been proposed and demonstrated that faceting and segregation in alloys are strongly interrelated and are the result of two coupled energy-reduction mechanisms.
As has already been mentioned, the 3{111} and both the A-3 and S-3 GBs have very different energies (i.e., 0.4 vs 55 and 69 meV/Å 2 , respectively). This is consistent with experimental evidence that flat 3 tilt boundaries with the 110 rotation axis facet toward {111} and {112} boundaries [13]. However, it is not clear whether this faceting behavior of Si GBs is an intrinsic property of the GBs, i.e., driven only by the anisotropy in the boundary energy with respect to the inclination angle or mediated by impurities/solute atoms. The presence of stable/metastable reconstructions of the 3 {112} GBs increases the complexity of this system: For small inclination angles with respect to the {111} plane, the boundary area of the low-energy {111} facets will be considerably larger than the area of the {112} facets. At these angles, the interaction between the partial dislocations inherent to the junctions with A-3 facets may dominate and constitute the S-3 energetically preferential facets. The latter implies that transitions between facet junction reconstructions may occur, depending on the inclination angle and the facet period. This would drastically influence the electronic properties and the impurity segregation, as well as the GB mobility of faceted 3 GBs. Therefore, a phase diagram of this boundary system that would describe the facet junction reconstruction as a function of the inclination angle and the facet period is critical to understand and control the properties of these boundaries.
To address under what conditions faceting occurs, we investigate the energetics, atomic structure, and strain distribution of flat and faceted 3 tilt GBs in Si with the 110 rotation axis by employing DFT and large-scale MEAM potential calculations. In a first step, we have therefore parameterized a MEAM Si potential. For this, we use a material properties database derived from our as well as previously reported DFT calculations as well as reported experimental data. Our calculations reveal that faceting of tilt 3 GBs is an intrinsic property of Si and can occur even in the absence of impurity segregation. However, line junction and facet reconstructions are found to sensitively depend on the inclination angle and the facet length.
The paper is organized as follows. In Sec. II, the methodology is presented. Section III presents the results on the energetics, the atomic structure, and the strain associated with faceted 3 tilt GBs with 110 rotation axis. Based on these calculations, a phase diagram is derived that describes and predicts the facet and line junction reconstructions as a function of facet period and inclination angles. This phase diagram is shown to provide implications on GB faceting and mobility, impurity segregation, and electronic properties. Finally, Sec. IV summarizes the results. Appendix describes the parametrization of the MEAM potential.

II. METHODOLOGY
All first-principles calculations were performed by employing DFT and the projector augmented-wave method [20,21]. Both local density approximation (LDA) and generalized gradient approximation (GGA) were employed. A kinetic energy cutoff of 450 eV was used for the expansion of the plane-wave basis set. An equivalent of a 6 × 6 × 6 Monkhorst-Pack k-point mesh for the bulk unit cell was used to sample the Brillouin zone. Atomic positions were relaxed until the absolute value of the maximum force on all atoms was less than 0.01 eV/Å.
The generalized stacking fault energies (GSFE) were calculated using a diamond crystal with the three primitive vectors oriented along 110 , 111 , and 112 . LDA was used for exchange and correlation. Supercells consisted of ten unit cells along the surface normal and a 1 × 1 interface cell. All DFT calculations of the {111}, S-, and A- {112} 3 GBs were performed using supercells consisting of two mutually compensated GBs and with a 2 × 2 interface cell. The separation distance between the two interfaces was larger than 30 Å. This separation was validated to be large enough to decouple the two boundaries and to attain boundary energies with an accuracy better than 0.1 meV/Å 2 .
All interatomic potential calculations were performed using the molecular dynamics simulator LAMMPS [22]. The atomic geometries and energies of flat and faceted GBs were obtained by simulated annealing within the isothermalisobaric ensemble ( NPT). The temperature of the system was first raised to 700 K. It was then cooled to 0.1 K with a cooling time of 750 ns. We have explicitly checked that an order of magnitude slower cooling rate had no effect on both the energetics and the atomic geometries. The final structures and energies were then obtained by performing a conjugate gradient relaxation of the annealed structures. The conjugate gradient minimization was terminated either when the energy change was below 10 −6 eV or when all atomic forces were below 10 −6 eV/Å.

III. RESULTS AND DISCUSSION
In general, five macroscopic degrees of freedom are necessary to define and identify a GB. Three of them define the misorientation of the two neighboring grains. The other two define the boundary plane orientation. In addition, to achieve a complete macroscopic and microscopic description of the GBs, one has to also consider translations of one grain with respect to the other, the density, and the exact position of the atoms at the GB. An illustrative example that the aforementioned five macroscopic parameters alone are inadequate to provide a full description of the lowest energy structures of GBs is the 3 {112} GB in Si. The atomic geometry of this GB as derived by the CSL is shown in Fig. 3(c) and the relaxed structure in Fig. 2(a). It consists of five-and sevenatom-rings which are symmetrically aligned with respect to the (112) plane. One of the atoms at the boundary plane is overcoordinated, i.e., it has five rather than four bonds. However, high-resolution transmission electron microscopy investigations observed another interface reconstruction [23]. This reconstruction is shown in Fig. 2(b) and consists of five-, six-, and seven-atom rings. Unlike the aforementioned symmetric geometry, all atoms at the reconstructed 1 × 2 boundary plane are fourfold coordinated. It also lacks the mirror symmetry. This A-3 {112} GB in Si indeed has a smaller boundary energy than the symmetric one as shown by present and previous DFT calculations [23,24].
In the present paper, we focus on the energetics of 3 tilt GBs having the 110 rotational axis. This choice fixes four out of the five macroscopic degrees of freedom. This leaves the inclination angle of the boundary plane with respect to the {111} plane as the only free degree of freedom. To determine the boundary energy along this parameter, we follow two different approaches. In the first approach, bicrystals are constructed using the 3 CSL with the 110 rotation axis. We consider four different boundary plane inclination angles Thus, periodic boundary conditions can be applied in any or all of the three directions [25]. In the second approach, the two grains are terminated by two free surfaces, which are parallel to the GB. Periodic boundary conditions are applied in the GB plane. The thicknesses of the grains normal to the GB range from ≈100 Å to ≈150 Å. This thickness has been found to decouple the free surfaces from the boundary.
Next to the five macroscopic degrees of freedom that are uniquely defined by the above approach, we need to minimize the GB energy with respect to the microscopic degrees, i.e., GB translations and atomic densities and positions. This issue is addressed by a heuristic approach: Using the above-mentioned bicrystals as input structures, we apply rigid translations on one grain with respect to the other in both boundary plane directions. A step width of a x /20 and

083604-3
a y /4 is used. Here, a x and a y are the two CSL primitive vectors in the boundary plane with a x being along 110 . For each translation, we randomly remove zero to ten atoms from the 1 × 1 interface. After removing these atoms, we anneal the structure. We repeat the aforementioned step ten times. At each of these steps, different combinations of atoms are randomly removed from the initial structure. The annealing procedure effectively overcomes any kinetic barriers that may prevent the system to find its low energy structure. At the end of each annealing step, conjugate gradient-based atomic relaxation is performed.
The above approach determines for each of the inclination angles the GB energy within a grandcanonical ensemble. It is calculated as where E tot is the total energy of the bicrystal, n denotes the number of Si atoms, μ Si is the chemical potential of bulk Si, and A is the GB area. To avoid the ambiguity due to the presence of the free surfaces, we have evaluated the total energy as the sum of individual atomic energies excluding all the atoms residing within 5 Å from each surface. We have explicitly checked the effect of the free surfaces on the calculated GB energies. We find that the aforementioned approach yields GB energies with an accuracy better than 0.01 meV/Å 2 . In Fig. 4(a), the boundary energies derived by this approach for a GB with an inclination angle of 35.26 • are shown. The calculated boundary energies span a wide energy range from 0.04 eV/Å 2 to 5.0 eV/Å 2 [for the sake of clarity, the energy scale in Fig. 4(a) is truncated to 0.20 eV/Å 2 ]. The lowest energy boundary structure is commensurate to the CSL, i.e., it has no rigid shift and no atoms are removed from the initial CSL based geometry. Let us have a closer look at the relaxed and unrelaxed geometries of the lowest energy boundary [see Fig. 3(d)]. Although we started from a flat boundary, the boundary spontaneously dissociates into {111} and S-{112} facets. The final structure can be obtained from the initial one by slightly displacing six atoms per 1 × 1 interface cell. The displacement vectors are smaller than the bulk lattice constant [see the red arrows in Fig. 3(d)]. The fact that displacements are relatively small indicates that kinetic barriers that may hinder faceting can be overcome by annealing. Furthermore, by inspecting the lowest energy structures for the other three inclination angles (i.e., 43.31 • , 62.06 • , and 70.53 • ), we find the same behavior: All these boundaries are energetically unstable against dissociation into {111} and S-{112} facets. Facet periods can be as small as a single interface primitive vector. Hence, this group of GBs is intrinsically unstable against nanofaceting.
In general, faceting of GBs may be associated with dislocations and line forces. Thus, the facet length is determined by the competition between the former and the latter [1]. The character of the dislocations and their Burgers vectors at the facet junctions are determined by the rigid body translations of the two facets. The line forces arise from the different interface stresses at the two facets at the junction. Therefore, the energy of a faceted GB as a function of the facet period L can be written as [1,16] where g f1 and g f2 are the energy contributions of the two facets. In the case of faceting toward {111} and {112}, these are written as where g {111} and g {112} are the boundary energies of the flat {111} and {112} GBs and θ is the angle between the inclined boundary and the {111} plane. The last two terms on the right side of Eq. (2) describe the interaction energies between dislocation and line forces and include contributions from the core energies. A, B, and C are constants that depend on the elastic constants of the material and the inclination angle, as well as the energy and diameter of the cores. A schematic representation of a faceted GB, indicating the inclination angle and the facet period, is shown in Fig. 1. In Fig. 4(b), the energies of faceted GBs with S-{112} facets for four inclination angles alongside with the facet contributions, i.e., the sum g f1 + g f2 from Eq. (2), which account for the core energies and the long-range interactions between dislocations and line forces, are expected to be negligible. In contrast, in the case of A-3 facets these terms are expected to be considerable and dominate at small {112} facet lengths.
To shed light upon the aforementioned long-range interactions, we compute the elastic strain tensor at the facet junctions. We therefore implemented large supercells consisting of ≈10 5 atoms and a pair of faceted GBs consisting of 10 and 20 units of {111} and {112} facets, respectively. In Figs. 5(a) and 5(b), we plot the ε xx component of the elastic strain tensor for facet junctions consisting of {111} and Sand A-{112} facets, respectively. As can been seen in the case of S-{112} facets, ε xx is negligible everywhere but at the {111} and {112} boundaries. The tensile strain at the {111} boundary, which corresponds to a hexagonal inclusion in the Si diamond matrix, is due to the larger interatomic distance in hexagonal Si. In the case of A-{112} facets, the strain shows a long-range behavior. It is in qualitative and quantitative agreement with the strain field calculated by elasticity theory for two mutually compensating edge type dislocations placed at the two junctions [see contour lines in Fig. 5(c)].
To quantify the impact of the long-range interactions of the dislocations on the GB energetics, we construct faceted junctions consisting of {111} and A-or S-{112} facets. This is achieved by inclining the {112} boundaries toward {111}. The surpercells then contain a pair of symmetry equivalent and mutually compensated faceted GBs. The separation distance between the two inclined GBs is at least 100 Å. We have explicitly checked for convergence at each angle of inclination by varying the length of the both facets. As in the first approach, the atomic geometries are relaxed by performing simulated annealing and subsequent conjugate gradient atomic relaxation as described in Sec. II. This prevents the system to get trapped in high energy metastable states. The GB energies as function of the facet period are plotted in of the facet length. This is in agreement with the finding that, in this case, contributions to the boundary energy arising from dislocations and/or line force interactions are either negligible or cancel out. However, faceted GBs with A-{112} facets show strong dependence on the facet period: Larger boundary energies are calculated for small facet periods. They asymptotically converge to the g f1 + g f2 energy term at long facet lengths [see Eqs. (2) and (3)]. These findings indicate that finite length faceting toward {111} and A-{112} boundaries is energetically unfavorable. Therefore, under conditions of thermodynamic equilibrium, the facet length will reach the maximum one that is possible due to the conditions given the microstructure and grain size. This behavior is similar to the case of GBs in Al where it was found that the GB tension is not sufficient to stabilize finite facet lengths [16,26].
Another important result that can be deduced from Fig. 6 is that with decreasing facet length a transition from A-{112} to S-{112} facets takes place. Although the g f2 energy contribution of the A-{112} boundaries are smaller than S-{112} for all inclination angles, at small facet periods contributions to the boundary energy arising from the dislocation interactions dominate. The latter constitute facet junctions with S-{112} facets energetically favorable at small facet periods. This is more systematically shown in Fig. 7(a), where the energy range of faceted GBs consisting of A-{112} facets and facet periods from the minimum symmetry allowed one to the asymptotic limit of infinitely long ones, has been been marked by the yellow shaded region. As can be seen faceted GBs with A-{112} facets are energetically favorable for all facet periods only for high inclination angles, i.e., angles larger than ≈80 • . For all other inclination angles, the thermodynamically most stable reconstruction depends on the facet period.
Our results are summarized in Fig. 7(b), which shows the energy difference between the two phases, i.e., between faceted inclined GB consisting of {111} and A-{112} (g A− 3 ) and S-{112} (g S− 3 ) facets. To estimate the impact of temperature on the phase boundary, we computed the vibrational contributions to the free energy of flat S-and A-{112} GBs. We followed the methodology that has been outlined in Ref. [27]. Including vibrational contribution to the free energy of the GBs increases the boundary energy difference between A-and S-{112} GBs by ≈2 meV/Å 2 at 1000 K with respect to to the 0 K results. Including them would result only in a small upward shift of the phase boundary in Fig. 7(b). To give a specific example, at an inclination angle of 45 • a ≈2 nm upward shift toward larger facet periods is estimated. Since vibrational effects are negligible, we do not included them in the phase diagram.
The phase diagram in Fig. 7(b) reveals that the reconstruction of the facets as well as of the line junctions depend strongly on the two geometrical characteristics of the facet junctions, i.e., the inclination angle and the facet period. The transition from one phase to another has important consequences on (i) the electronic properties of faceted GBs in mc-Si, (ii) the impurity segregation at the facets as well as at the line junctions, and (iii) the GB mobility. Previous studies showed already that overcoordinated atoms and highly strained bonds at the S-{112} GBs introduce deep states in the fundamental band gap [10]. In contrast, the reconstructed A-{112} GB shows no deep states in the fundamental band gap [10]. Hence, process conditions that stabilize the S-{112} reconstruction are detrimental for the electronic properties of mc-Si. Thus, conditions that lead to low inclination angles and/or small facet lengths should be avoided. In an experimental or industrial setup, this could be achieved by long annealing that drives the system toward thermodynamic equilibrium and thus large grains.
The elastic strain of both S-and A-{112} GBs show a periodic pattern with alternating compressive and tensile regions along 111 [see Figs. 5(a) and 5(b)]. Explicitly studying segregation of impurities to a GB is beyond the scope of the present paper. However, the results of impurity-free GBs provide some qualitative insight into impurity segregation. In this respect, it is important to note that the ε xx strain extends deep into the bulk at the S-{112} facets. Thus different induced segregation profiles at S-and A-{112} GBs are expected [10,28]. Hence, the segregation profile will be a function of inclination angle and facet period and will be thus qualitatively described by the phase diagram if Fig. 7(b). The strain field at the {111} and S-{112} line junctions is short range compared to the range of the strain field associated with flat S-{112} GBs. On the other hand, the strain field at {111} and A-{112} junctions is long range and originates from the presence of edge-type dislocation at these junctions. The latter, in combination with the lack of mirror symmetry of the A facet, results in two different topologies and impurity segregation potentials at the two line junctions of these GB facets. Indeed, a recent study showed that the aforementioned difference in the topology of the facet junction in mc-Si causes asymmetric impurity segregation at {111} and A-{112} facet junctions [13].
GB disconnections are GB line defects. Generally, they have both step and dislocation character and play a crucial role on the shear coupled migration [29]. Recent TEM investigations on Au bicrystals revealed the lack of lattice defects at the core of disconnections with a height of two lattice spacings at 11(113) coherent GBs [30]. These boundaries showed a layer by layer migration. It was also reported that the shear-induced migration of these GBs was fully reversible and not affected by the presence of preexisting lattice or GB dislocations. This situation is in accord to the twin defects in Si where disconnections are {112} steps at the {111} boundary. This corresponds to small facet periods or small inclination angles. Hence, as it is shown in Fig. 7(b), the steps consist of S-{112} facets, i.e., there are no GB dislocations at the line junctions. Therefore, a similar mechanism as described is expected to govern the shear induced migration of twins in Si, i.e., migration will be fully reversible and will not be affected by the presence of other extended lattice defects.

IV. CONCLUSIONS
In the present paper, we have employed DFT and largescale MEAM potential calculations to study the structure, energetics, and strain associated with 3 tilt GBs having the 110 rotation axis in Si. Based on these calculations, we derive a phase diagram of these boundaries. This phase diagram demonstrates that (i) these GBs are intrinsically unstable against faceting toward {111} and {112} facets and (ii) the properties of the facets and of the line junctions is the result of an intricate interplay between GB energies and long-range strain interactions. Specifically, we find that at low misorientation angles and/or small facet periods, longrange interactions dominate and S-{112} facets are favored.
Nevertheless, at large facet periods and inclination angles, the lower energy of the A-{112} facets compensate the strain interactions and these facets are energetically favored.
The significance of being able to construct such defect phase diagrams goes beyond the case of GBs in mc-Si. They shed light upon GB faceting phenomena, specifically in conjunction with the materials' microstructure. The picture deduced by this phase diagram contradicts the common perception that the properties of faceting are merely driven by the anisotropic GB energies. Although anistropic GB energies are a prerequisite for GB faceting, the phase diagram reveals that higher energy metastable GB phases are stabilized by thermodynamics and not kinetics when constituting the facets at line junctions. This insight provides critical information when designing electronic and structural materials: Microstructures that allow for large facet periods stabilize line junctions accommodated by dislocations. In contrast, fine granular structures that limit facet periods promote line junctions without extended line defects and long-range strain fields. This also has immediate implications on the topology of line junctions, the strain-driven impurities' segregation, and mobility of GBs. In general, GB phase diagrams that highlight the interplay between geometry, GB energies, and long-range strain fields constitute an indispensable tool in describing, predicting, and designing materials' microstructure.

ACKNOWLEDGMENT
The authors would like to acknowledge fruitful discussions with Dr. Christian Liebscher.

APPENDIX: MEAM POTENTIAL PARAMETRIZATION
A large number of valence force field models have been developed and applied to study defects in Si-among others, Tersoff [32], Stillinger-Weber (SW) [33], bond-order potentials [34,35], and MEAM potentials [36,37]. However, application of these potentials to study important materials properties such as GSFE reveals critical shortcomings. Specifically, they fail to provide a quantitative and qualitative description of GSFE curves with respect to DFT (see Fig. 8). The GSFE is a material property which relates dislocation cores and the intrinsic ductility of the material. Smooth GSFE curves can also be considered as a benchmark of the interatomic screening under shear conditions.
The second-nearest-neighbor (2NN) MEAM potential is a modification of the original MEAM potential [31] and partially include 2NN interactions. 2NN MEAM potentials have been parametrized and employed to a wide range of elements and alloys (see Ref. [38] and NIST Interatomic Potentials Repository [39][40][41]). In the present paper, we parametrize and employ a 2NN MEAM potential to investigate faceting of GBs in Si. A detailed description of the 2NN MEAM formalism can be found in Ref. [42]. Here the parametrization of the potential is presented. Although the correlation of the parameters to physical properties is complicated, some parameters play a more dominant role on certain physical properties [43]. This allows us to apply a parametric study and systematically fit the parameters to these properties.  [31], SW [33], bond-order [34], and Tersoff potentials [32].
The 2NN MEAM potentials formalism consists of 16 parameters for a pure element. These are listed in Table I. The cohesive energy E c , nearest-neighbor equilibrium distance r e , and parameters α and δ enter the zero-temperature Rose's universal equation of state. α is also a function of the cohesive energy and bulk modulus B. The remaining 12 parameters include the four decay lengths of the atomic partial electron densities (β i , i = 0, 1, 2, 3), the three scaling factors of the background electron density (t i , i = 1, 2, 3), the scaling factor of the embedding function energy (A), the two parameters which control many-body screening (C min and C max ), and the cutoff distance (r c ) and cutoff range ( r) of the atomic electron densities. The parametrization of the first-nearestneighbor (1NN) MEAM potential for Si by Baskes is used as the basis for the present potential [31]. In the present paper, as in the original Baskes parametrization, the values of E c , r e , and of the bulk modulus B are considered from experiment [44][45][46]. C max , β 1 − β 3 , α, t 1 , r c , and r are set to the same value as in the 1NN MEAM formalism [31].
It has been shown that δ can be fitted to the pressure derivative of the bulk modulus [43]. In the present paper, the MEAM-calculated pressure derivative of the bulk modulus is in good agreement with DFT calculations, i.e., 4.15 vs 4.20, respectively, even with neglecting the δα * 3 term in the Rose's equation of state. Therefore, in the present parametrization, δ is set to 0.0.
C min and β 0 are used to fit the potential to the DFTcalculated GSFE curve for shear along 110 in the {111} shuffle plane. The GSFE curve along 112 in the {111} shuffle plane is used to benchmark the potential. As can be seen in Fig. 8, the present potential provides an excellent description, with respect to DFT, for both shear configurations. On the contrary, the original 1NN MEAM potential, the SW, and the bond-order potentials fail to provide both quantitative and qualitative agreement with DFT. Specifically, the present potential, in agreement with DFT, shows the possibility of dissociated dislocation having a stable stacking fault energy while the bond-order potential shows a flat maximum. The 1NN MEAM and SW potentials yield discontinuous profiles and/or spurious local minima in the region of the unstable stacking fault.
Our parametric study reveals that spurious oscillations in the GSFE curve found in the previous 1NN MEAM Si potential (see Fig. 8) can be removed by lowering the value of C min from the value of 2.0 in the original 1NN MEAM [31]. However, this increases the GB energies with respect to the DFT calculated values. Therefore, there is a trade-off between two important properties in the description of GB faceting: The quantitative and qualitative description of GB energies and the interatomic screening at the region of the facet junctions where the core of the partial dislocations is located. We have explicitly checked that the aforementioned increase in the GB energies caused by the lower values of C min has no effect on the description of the long range interactions of the line junctions.
Parameter t 2 was determined by fitting the elastic constants of Si. The elastic constants are calculated using volume conserving strain deformations [47]. The elastic constants calculated by the present 2NN MEAM (P-MEAM) are in excellent agreement with the DFT calculated as well as experimentally obtained values (see Table II). Parameters A and t 3 are used to fit the cohesive energies of diamond, fcc, bcc, and β-tin Si structures. The cohesive energy differences of these structures with respect to diamond Si are listed in Table II. Although, the present MEAM potential overestimates these energy differences, it provides qualitative agreement to DFT. A transition from diamond cubic to β tin crystal phase of bulk Si takes place at high compressive pressures of ≈10-13 GPa [48,49]. Although compressive strains are present at the {112} GBs [see Figs. 5(a) and 5(b)], the overestimation of the β tin cohesive energy does not affect the good qualitative description of the 3 GB energetics by P-MEAM. This has been validated by (i) ensuring good quantitative description of the   [31], bond order (E-bond order) [34], and SW [33] potentials compared to DFT and experimental data. All the properties are for diamond cubic Si unless otherwise denoted. E dia is the cohesive energy of Si in the ground-state diamond crystal and E denotes the cohesive energy differences between the diamond and the fcc, bcc, and β-tin Si structures in meV/atom. The lattice parameters a are in Å and the bulk modulus and the elastic constants in GPa. The energies of the 1 × 1 (110), (110), and (111) surfaces are denoted as E (100) , E (110) , and E (111) , respectively, and the GB energies of the {111} and S-and A-{112}, denoted as E 3{111} , E S− 3{112} , and E A− 3{112} , respectively, are in meV/Å 2 . The vacancy formation energy E v is in eV. GSFE including a large number of highly deformed structures at the interface and (ii) demonstrating a qualitative description of the energetics and the atomic volumes of the competing bulk crystal phases. These validation benchmarks show that P-MEAM provides an excellent qualitative description of 3 GB energetics and atomic structures. Using the P-MEAM, we calculate the GB formation energies of 3 {111} and {112} tilt GBs with the 110 rotation axis. The GB energies of {111} and S-{112} are 0.4 and 69.89 meV/Å 2 , respectively. The GB energy of the A-{112} is by 14.59 meV/Å 2 lower than its symmetric counterpart.
As has already been discussed, the P-MEAM calculated energies of the {112} GBs are overestimated. However, the GB energy difference between A-and S-{112} is in excellent agreement with the DFT results (17.38 meV/Å 2 ). Moreover, the P-MEAM calculated elastic constants are in excellent agreement to DFT. These agreements imply that the strain interaction between the facet junctions as well as transitions between S-and A-{112} GBs as function of inclination angle and/or facet length are correctly described by P-MEAM [see Fig. 7(b)]. Thus, the present MEAM potential provides an excellent qualitative description of faceting.