Free‐Energy Parameterization and Thermodynamics in Si–Ge–Sn Alloys

Different models are used to describe and parameterize enthalpy and free energy in GeSnSi alloys over the entire composition range from density‐functional theory calculations and are used to explore if concentrated alloys in the equiatomic concentration range are thermodynamically feasible. It is found that the energetics are very well behaved and follow the expressions of a ternary regular solution rather closely. While a small asymmetry is found in the ideal solution energetics for the three binaries, it is found that the assumption of a symmetric ideal solution as underlying the quasichemical model does not change the general findings. All free energy contributions within this work are parameterized in different approximations, providing a complete parameter set that can be used for further thermodynamic studies of the GeSnSi alloy system. In the obtained results, the enthalpy of mixing of the SiSn system has rather large values, which is maintained more or less unchanged in the ternary system, making a concentrated ternary alloy thermodynamically unfavorable. Calculating the full free energy from adding configurational and vibrational entropy contributions does not change the overall situation.

hoped to extend the accessible bandgap range. Furthermore, while the Hume-Rothery rules gave little hope for stable equilibrium alloys, the addition of carbon in small, meta-stable concentrations was hoped to compensate for the tensile strain caused by the addition of the larger Ge to Si in the extensively studied pseudomorphic SiGe/Si system. [8][9][10] Subsequently, it was found that neither the lattice constant nor the bandgap followed a simple Vegard's law like interpolation scheme. [6,11] Nevertheless, the research on SiGe:C alloys with small amounts of C led to the development of a new device, the heterojunction bipolar transistor, by Motorola in early 2000, still in use in high-performance wireless applications.
SiGeSn alloys have only recently become the focus of research, starting with the density-functional theory (DFT) prediction that holes in GeSn alloys may have even higher mobility than their peers in pure Ge. [12] However, even more exciting was the prediction that a direct bandgap with a conduction-band minimum at the Γ-point could be achieved in GeSn alloys once a minimum concentration of Sn between 6 and 11 at.% (depending on the prediction) [13][14][15][16][17] was reached in the alloy. This is significant since all other crystalline group-14 semiconductors and alloys have indirect band gaps (unless they are turned into functionalized two-dimensional crystals [18] ). The direct bandgap would provide an exciting new convenient materials system for optoelectronic applications with opportunities in the small-gap infrared range. Indeed, recent work is suggesting that the bandgap of GeSn alloys displays a continuous evolution from indirect to direct through band mixing effects with a significant direct component already at Sn concentrations around 6%, as opposed to a distinct indirect-to-direct transition from L to Γ conduction-band minimum such as seen in III-V alloys. [19] Besides the direct-gap character being already accessible through these band-mixing effects in GeSn, it has also been speculated that there may be a "magical" composition range in the SiGeSn ternary system that may be stable (due to strain compensation between the large Sn and the small Si atoms), especially for strained thin films, where a direct, nonzero band gap is expected. [20] While DFT calculations for ternary alloys with a Si:Sn ratio of ∼4 that lattice-matched the film to Ge suggested that the optical onset had a tail-like shape reminiscent of the Urbach tail observed in amorphous materials, [21] recent experiments on a Ge 0.90 Sn 0.04 Si 0.06 found a reduction of the indirect bandgap absorption at lower energies and of the Urbach energy. [22] However, very recent theoretical suggestions of high-entropy alloys like equiatomic alloys [23] have directed attention more into the center of the ternary phase diagram, supported by recent work that has demonstrated that films with Sn concentrations up to 30% seem indeed possible to grow. [24] Lastly, while the ordering effect of C-atoms on third-neighbor positions in SiC alloys, which had been predicted by classicalforce field simulations and DFT calculations, [25][26][27] and DFTbased Monte Carlo (MC) simulations [11] could be unambiguously verified experimentally in MBE-grown samples, [25,26] an unambiguous experimental proof of predicted ordering effects in GeSn alloys [28] is still missing. Since the identified driving energies for ordering in GeSn alloys are much smaller than in SiC and similar to the entropy contributions at growth temperature, adding those into the calculation may be much more important than for SiC alloys. The ordering predictions for SnGe [28] are based on MC equilibration, which inherently includes configurational entropy, but no vibrational entropy. Calculating the vibrational contributions is, however, rather tedious and would slow the already demanding MC simulations down to impracticability, which makes a fast method to also include the energy contribution of vibrational entropy desirable.
In this article, we thus want to introduce a number of computationally efficient parameterizations for enthalpy and entropy contributions to the free energies of GeSnSi alloys based on full DFT calculations, which will close this methodological gap. We then will use the developed model to explore the solubility and energetic stability of SiGeSn alloys with a special focus on the equiatomic composition range. We will validate our approach in comparison to structural data, and finally examine mixing enthalpies, equilibrium phases, ordering effects, and full free energies of mixing to understand the thermodynamic stability of SiGeSn alloys to answer the question of if specifically equiatomic Si-Ge-Sn alloys, where a "magical" composition has been frequently speculated to reside, [23] have a thermodynamic chance to exist. Since recent work suggests that strain relaxation through plastically relaxed buffer layers is conducive to higher Sn solubility, [20,[29][30][31] we focus here on the unstrained bulk system.

Computational Approach
In this work, we use DFT-based calculation as implemented within the Vienna Ab-initio Simulation Package (VASP). [32,33] After careful validation against experimental results, we found that PAW potentials [34] within the local density approximation (LDA) [35] successfully reproduced the lattice constants as well as the phonon densities of states [36] needed for the vibrational free energy, [37] while we will show that the mixing energetics are better described by Perdew-Burke-Ernzerhof (PBE) functional calculations. For the kinetic energy cut-off, we used a value of 250 eV. For the 64-atom conventional supercells used in this work, we use 4 Â 4 Â 4 Monkhorst-Pack meshes [38] with Methfessel-Paxton smearing for everything except the final energy calculations, where we employed the tetrahedron method with Blöchl corrections. All atomic positions and lattice constants were fully relaxed.
Since different DFT codes use different baselines for the total energies they calculate, we report here cohesive energies E coh for easier comparison as well as since it is a more sensible quantity when determining bond energies, since the sum of all bond energies in a system is equal to its cohesive energy. E coh is calculated from the raw total energies for the 64-atom special quasirandom structures (SQS) cells minus the energies of the isolated atoms in vacuum where we find with our PBE settings E Ge ¼ À0.943 eV, E Sn ¼ À0.821 eV, and E Si ¼ À0.996 eV, while LDA results in E Ge ¼ À0.759 eV, E Sn ¼ À0.661 eV, and E Si ¼ À0.801 eV.
For the energetics and lattice constants, we tested the influence of spin-orbit coupling (SOC), which is especially important for the band structure due to the low position of Sn in the periodic system. However, for both energetics and lattice constants, we will show that the influence of SOC is negligible. Because of that, we perform all subsequent calculations without spin-orbit coupling.
www.advancedsciencenews.com www.pss-b.com In this work, we want to foremost study random alloys, their energetics, and eventually also their electronic properties. To do that with DFT, which requires small cells, the currently most popular strategy is to use SQS, [39] which we produced with the mcsqs code of the ATAT package. [40] One of the most important things that need to be reproduced when simulating a random configuration is the bond distribution. Analytically, one would expect a parabolic dependence of the bond fractions on the composition, specifically, for a system with a total number N tot of bonds which is also valid for systems with more than two components. [41] Figure 1a shows the ideal parabolic distributions along with symbols that represent the bond count in the 64-atom SQS cells used in this work, showing that the bond distribution follows nearly perfectly in the ideal case. This allows us to examine the energetics of the alloys with ensemble-averaged bond energies within the quasichemical model described in Sec. 3.4.
The SQS cells examined include for all three binary systems 39 structures with compositions A 64Àn B n , n ¼ 0, 2, 4, · · · , 64 (33 structures), plus n ¼ 1, 3, 5, 59, 61, 63 to have a finer mesh at the elemental ends. Furthermore, there are 75 ternary systems. These consist, on the one hand, of small-concentration additions of one component to the binary system of the other two components, specifically 7 Ge n Sn 2 Si 62Àn structures with n ¼ 8, 15, 23, 31, 39, 47, 54, and 7 Ge n Sn 4 Si 60Àn structures with n ¼ 8, 15,23,30,37,45,52, as well as their permutations over all components, resulting in 42 structures. Additionally, we use ternary alloys from one corner of the ternary phase diagram to the midpoint of the opposite side, i.e., 11 structures Ge n Sn 64À2n Si n with n ¼ 1, 2,4,8,12,16,20,24,28,30,31 and permutations for additional 33 structures. To demonstrate that the SQS cells used reproduce the familiar "negative" bandgap behavior in GeSn alloys, we show their calculated band gaps versus Sn concentration in Figure 1b, which like previous calculations by e.g. Polak et al. [42] do not really follow a parabolic concentration dependence, but a more complex one which we model here with a 6 th -degree polynomial.

Lattice Constants
As the first validation of our chosen methodology, we determine the lattice constants for SnGe alloys on the Ge-rich side, since earlier calculations had somewhat underpredicted the experimental values. Figure 2a shows the LDA-calculated lattice constants in comparison to experiments [43] and previous DFT calculations, [44] where our settings produce an excellent agreement. We also find a theoretical linear dependence of the lattice constant on Sn-concentration of aðGe 1Ày Sn y Þ ¼ ð5.648 þ 0.8603yÞ Å. Figure 2b shows the deviation of the calculated lattice constants from Vegard's law, [45] which is the straightline interpolation between the Ge and Sn lattice constants. For concentrations up to 32%, we find small overbending, while the deviation is negative for larger Sn concentrations. Overall, up to 70%, the deviation can be described by a cubic fit, ΔaðyÞ ¼ ð0.057y 3 À 0.099y 2 þ 0.026yÞ Å. For this calculation, we also checked the agreement between calculations without and with SOC (which slows down the calculations by a large factor). Within some small numerical oscillations, the symbols follow the same dependence, and the fitted cubic functions are nearly on top of each other. Therefore, for structural properties, we do not find it necessary to include SOC (and spin polarization in general) in the calculations.  [65] in comparison to previous modified Becke-Johnson potential results from Ref. [42] (red squares).

Mixing Enthalpy in Random Alloys
www.advancedsciencenews.com www.pss-b.com from the actual energy of an alloy, normalized per atom, which then is determined by where N ¼ l þ m þ n is the number of atoms in the supercell. Since mixing enthalpies are crucial for alloy thermodynamics and need to be determined accurately, we examined again the significance of SOC as a first test. As shown in Figure 3a, SOC plays a negligible role here as well, allowing to perform all energetic calculations without spin. For the mixing enthalpy in all three binary alloys at zero external stress in Figure 4, we find a slight non-parabolic dependence on composition, which can be well described by a regular solution expression with a regular solution parameter H r ðxÞ ¼ H 0 ð1 þ H 1 xÞ with linear concentration dependence The magnitude of H 1 determines the asymmetry of the enthalpy of mixing, which for H 1 ¼ 0 has its extremum (here, maximum) at x max ¼ 0.5 while for H 1 6 ¼ 0 Equation 5 can be extended to multicomponent alloys in an additive fashion, which reads for the case of SiGeSn ternary alloys where we have introduced the composition vector The concentration variables in the brackets have been chosen with respect to minimizing the object function of the regression. While for LDA, we just considered the separate binary alloys for fitting, for PBE, we tried fits both for binary-only structures and for all 192 structures with ternaries included.
As can be seen from the results for the fitted regular solution parameters, energies, and maximal concentrations shown in Table 1, we found little difference between binary-only and comprehensive PBE fits. Equally, as can be seen in Figure 5a,b, the parity plots for the fitted cohesive energies as well as the enthalpies of mixing show excellent fitting quality for the data set that includes both binary and ternary SQS cells, with E fit coh ðE DFT coh Þ ¼ ðÀ0.1 AE 0.2Þ þ ð0.9995 AE 0.0009ÞE DFT coh and ΔH fit ðΔH DFT Þ ¼ ð0.0002 AE 0.0007Þ þ ð0.996 AE 0.009ÞΔH DFT , where the errorbars bracket the 95% confidence interval of the regression result. The corresponding mixing enthalpy plots for the ternary systems as shown in Figure 5c,d show consequently equal fitting fidelity as found for the binary alloys shown in Figure 4. Summarizing these points, we find that the ternary mixing enthalpy can indeed be split into the pairwise terms as assumed in Eq. (7) and that (a) (b) Figure 2. a) Calculated lattice constants versus experimental lattice constants (black diamonds [43] ) and previous calculations (blue triangles [44] ) of a GeSn alloy. b) Deviation from Vegard's law [45] with and without spin-orbit coupling (SOC). Red solid and blue dashed lines are cubic fits. there are no zones in the ternary phase diagram that do not follow the (asymmetric) regular solution energetics. The concentrations for the maxima of the different curves are very similar and differ by less than 2% between LDA and PBE. Also, the energy sequence between the different binaries is similar in both approximations, yet with some quantitative differences. SiSn alloys have in both approximations a very high mixing enthalpy, with LDA finding a 9% lower value, while GeSn alloys are higher in LDA by 10%. Although the maximum value of the mixing enthalpy in SiGe differs by 12 meV between LDA and PBE, both values are below room-temperature k B T of 25 meV, underlining the known thermodynamic stability of the alloy across the entire composition range. Thus, one would qualitatively expect similar thermodynamic predictions from both exchange-correlation functionals, although PBE typically produces mixing enthalpies and formation energies in closer agreement with experiments. [46] The calculated mixing enthalpies are overall in line with the expectations from the size mismatch between the three atoms discussed in Sec. 1. While the binary mixing energies remain basically unchanged when going to ternary alloys, which we will also confirm using the quasichemical bond-energy approach from Ref. [41] in Sec. 3.4, the rather high mixing enthalpy of SiSn inhibits ternary alloy formation with high Si and Sn concentrations, a situation that will also not change when temperature-dependent free-energies are considered instead of zero-temperature enthalpies as will be shown in Sec. 3.5.

Ordering Effects
The large size difference between Sn on the one side and Ge and Si on the other side, being 15 and 19% larger, respectively, reminds us of investigations on Si-C alloys around the end of the previous century, where ordering effects of the considerably smaller C atom in the Si matrix have been predicted and observed. Figure 6 shows the redrawn pair distribution function between the C atoms for a Si 0.8 C 0.2 system before and after MC annealing for 6000 steps at 1400 K from Ref. [11]. It was found that after annealing, there were no longer C atoms in nearest neighbor positions. At the same time, the peak at the third neighbor distance contained significantly more atoms than the one at the second neighbor distance, indicating third-neighbor ordering. The third neighbor positions are at opposite sides of one of the hexagonal rings in the diamond structure in [110] projection and allow maximum strain relaxation as discussed previously. [26] The ordering hypothesis was confirmed in annealed samples by Raman spectroscopy through the presence of a third-neighbor peak, whose nature was identified by DFT-based modeling of the spectra. [26] Taking the third-neighbor arrangement to the limit, one can construct a supercell on the diamond lattice with A 4 B composition, where all B atoms are third-nearest neighbors to each other  www.advancedsciencenews.com www.pss-b.com as shown in Figure 8. In the late 1990s, it was speculated that such a structure may allow to create a metastable high-C content phase in the Si-C phase diagram, which shows hardly any miscibility and a very narrow line compound at 50% composition, SiC. Kouvetakis et al. [47] developed a chemical vapor deposition-based synthesis route for Si 4 C with a precursor containing 4 Si and 1 C atoms, which was theoretically examined by Windl and Sankey, [48] and found that the completely ordered Si 4 C structure is 0.1 eV lower per C atom than the MC optimized semiordered alloy from Figure 6. Considering this previous work, the question arises if something similar might happen for GeSn and SiSn alloys. Indeed, previous work [28] has examined ordering effects in Ge 1Àx Sn x alloys. With cells too small to allow for full phase separation as should happen according to the phase diagram and thus enforcing a solid-solution, short-range ordering with strongly decreased Sn-Sn nearest neighbor and enhanced third-neighbor frequency was found. However, these simulations were run at 300 K, whereas nearly all Ge 1Àx Sn x alloys are grown between 300 and 400°C with a much higher thermal budget. While it can be left to debate if the kinetic energy of 25 meV per atom at 300 K would be high enough for sufficient diffusion during growth to reach thermodynamic equilibrium, it is worth noting that the short-range ordering is strongly decreased at actual growth temperatures. To demonstrate that, we have repeated the MC simulations of Ref. [28] with a parameterized energy expression based on Sn-Sn pair interaction energies up to third neighbors calculated from DFT in a 216-atom supercell shown in the inset in Figure 7, and modeling the energy from the Sn-Sn pair energies up to third neighbors while Ge-Ge and Si-Ge energies are set to zero. With this setup, we can model the results from Ref. [28] well, as we demonstrate here for the example of 14.01% Sn, where it was found that first, second, and third neighbor shells are occupied by 9, 91, and 130%. Figure 7 shows our results from MC simulations for 8000 atom cells for 10 9 steps. While the second neighbor occupation stays approximately the same up to typical growth temperatures around 600 K, the first-neighbor shell fills up to about 40%, while the third neighbor shell gets close to the occupation for a random alloy. In line with that, there is no evidence of Sn-Sn short-range ordering effects in the experimental literature. For example, an extended X-ray absorption fine structure investigation of the Sn local environment in Ref. [49] identifies the studied samples, grown at 320°C, as homogeneous random substitutional alloys. Equally, www.advancedsciencenews.com www.pss-b.com a neighbor analysis from atom-probe tomography on SiGeSn films in Ref. [50] finds a distinct repulsive interaction between Si and Sn, but perfectly random Sn-Sn coordination in samples grown at 350-475°C. Two other studies find from atom-probe tomography neighbor analysis perfectly random Sn-Sn coordination on the first four neighbor shells in both nanowires [29] and strain-relaxed films on a multilayer stack with graded Sn content, [30] both grown at low temperatures around 300°C, where it is argued that unstrained systems were conducive to random-alloy formation with high Sn concentration.
Although it is thus unlikely that ordering effects play an actual role, since the goal of our paper is to study if in the best-case scenario, a ternary concentrated Si-Ge-Sn alloy would be feasible, we calculated the formation enthalpy of the 80:20 structure shown in Figure 8 for Si 4 Ge, SiGe 4 , Ge 4 Sn, GeSn 4 , Si 4 Sn, and SiSn 4 . As can be seen in Figure 4, the energy gain with respect to a random alloy is maximal for the binary system with the largest size mismatch, SiSn, where ordering lowers the formation enthalpy per atom on the Si-rich side by 22 (24) meV and on the Sn-rich side by 14 (12) meV (numbers before (inside) brackets indicate LDA (PBE) results). For GeSn, the energy gain from the order is smaller with 14 (12) meV on the Ge-rich side and 7 (7) meV on the Sn-rich side, while for GeSi, which forms a nearly ideal solution, the energy gain is always less than 1 meV, specifically 0.3 (0.4) meV on the Si-rich and 0.7 (0.3) on the Ge-rich side.
Since SiSn alloys, which form the only binary system where ordering results in energy gain at least comparable to roomtemperature k B T of 25 meV, are found to have a mixing enthalpy high enough to make their formation rather unlikely as discussed in Secs. 3.2 and 3.5, while GeSn and SiGe have rather small energy gains from ordering which are not competitive with entropy contributions even at low temperatures, we do expect ordering effects to play a negligible role in the SiGeSn alloy system.

Bond Energies
Everything up to now was based on cohesive energies, averaged over all the atoms. However, we found it instructional to examine the values of the average bond energies in the alloys, which allow examining solubility and stability in multicomponent systems in Figure 7. The ratio of first-(blue), second-(green), and third-neighbors (red) between Sn atoms relative to a random alloy for a Ge 0.86 Sn 0.14 alloy in an 8000-atom cell that has been equilibrated with Monte Carlo (MC) simulations for 10 9 steps using the pair-potential approximation and Sn-Sn pair interaction energies up to third neighbors from the inset, calculated with DFT in a Ge 214 Sn 2 cell. The solid lines are second-order polynomial fits. Figure 8. Si 4 Sn in the ordered Si 4 C structure that was proposed in Ref. [48]. Blue (small) balls represent Si, gray (large) balls Sn atoms. Figure 6. Radial distribution function gðrÞ for the C atoms in a 40-atom supercell with 8 C and 32 Si atoms. The triangles denote the distances for first, second, and third neighbors, respectively, in an ideal structure. The dashed line is the distribution of the randomly distributed C atoms at the beginning of the ab initio Monte Carlo simulation, the solid line shows the distribution after 6,000 steps. Reproduced after Ref. [11]. The inset shows two C atoms in the third-neighbor position in Si, [110] plane.
www.advancedsciencenews.com www.pss-b.com comparison to the underlying binary systems from the angle of the well-known quasichemical solution theory, which we have recently expanded into a bond-order bond energy model based on DFT calculations. [41] Within the quasichemical solution theory, for the example of an AB binary alloy, specific average bond energy values are assumed to exist for bonds between A-A, A-B, and B─B atoms. With these three bond energies and the number of A-A, A-B, and B-B bonds in the alloy, the cohesive energy of a specific alloy system can be calculated. If one is interested in the energy of a specific atom, which is the enthalpy contribution to its chemical potential, it can then be calculated by adding half of the bond energies of all bonds connecting to that atom, where half of each bond energy is attributed to each atom at the end of the bond. [41] The quasichemical solution theory is based on ensemble averages of the bond energies, which scatter around the mean value in the real system based on local arrangement and strain. If one is interested in the properties of ideally random alloys, the classical random fractions of bond energies can be used as given in Eq. (2) and shown in Figure 1.
In our previous work, we have shown that by calculating DFT energies for a series of supercells with changing compositions and counting the bonds between the different species in them, the different cohesive energies can be well fitted. [41,51] We have also shown that in more polar or ionic structures where charge transfer significantly affects the bond strength and the bond energy between two elements changes when other elements are present, the introduction of a bond synergy parameter can appropriately modify the bond strength as a function of composition, [52] which we, however, found to not be necessary here.
The cohesive energy for Si-Ge-Sn alloys based on the quasichemical solution model [41] is given by where ε AB is the bond energy between atoms A and B. In case there are mutual solubility issues, we have also shown that this can be addressed by introducing a bond-order sigmoid function that allows using different bond energies on A-rich and B-rich ends for non-symmetric nonregular enthalpies of mixing. [41] As a first step, we will try fitting a straight quasichemical solution expression (which assumes a regular solid solution with parabola-like enthalpy of mixing) and then investigate if more elaborate bond-energy schemes are necessary. We start from the elemental bond energies, which we also visualize in Figure 9 for the LDA case, where the bond energy relative to the asymptotic value at very large bond length is plotted for Si, Ge, and Sn. Their LDA ground state bond energies from that plot are À2.65, À2.29, and À1.99 eV, respectively, while the corresponding PBE values are À2.21, À1.77, and À1.51 eV. Although the mixing enthalpies were comparable, our average cohesive or formation energy for the three elemental solids, given by twice the bond energy, is significantly lower for LDA in comparison to PBE by À0.99 AE 0.05 eV/atom, which is a well-known fact and within the standard deviation of typical benchmarking results, where LDA has been found to find lower formation energy than PBE by À0.58 AE 0.55 eV/atom. [46] Since the average deviation from experimental values for PBE with a value of À0.1 eV/atom is much smaller than for LDA, we will base the further discussion of energetics only on the PBE values.
Fitting the DFT cohesive energies for the three binaries, we find excellent agreement between fit and original DFT energies with a dependence of E fit coh ðE DFT coh Þ ¼ ðÀ0.2 AE 0.5Þ þ ð0.999 AE 0.002ÞE DFT coh , with only a small increase in uncertainty as compared to the results for the regular solution model from Sec. 3.2. The resulting parity plot is shown in Figure 10a. An even more sensitive test of the quality of the bond energies is a look at the mixing enthalpies from the bond energies in comparison to the straight DFT values, for which we again find excellent agreement for both PBE and LDA (Figure 10b,c, the latter shown for completeness). There is also excellent agreement between the bond-energy derived maximum mixing enthalpy, given by the difference between the average of the homobond and the bond energy of the heterobond and shown as ΔE in Table 2, and the values that were determined in Sec. 3.2 and shown in Table 1. Also, the elemental values for the same elements in different alloys (e.g., the Si-Si bond energy in SiGe vs SiSn alloys) agree with each other nearly perfectly within just a few meV.
Thus, while the quasichemical model sacrifices a few meV of accuracy in the description of DFT energies, its true value is in the fact that it allows in a very intuitive way to answer the question about the existence of a "magical" composition region in the ternary GeSnSi phase diagram where the local-strain compensation of Si and Sn or an electronic effect may result in a lowered mixing enthalpy. If that were true, the parity plot would show discrepancies, and calculation of the ternary mixing enthalpies from the unchanged binary bond energies should result in values that are significantly different from the straight DFT results. The quasichemical mixing enthalpies from the binary bond energies are compared to the DFT values in Figure 10e,f for two cases of interest, the one for ternary alloys with small Si concentration Figure 9. Bond energy as a function of bond length relative to the saturation large-bond length value for Si, Ge, and Sn within LDA, calculated from the total energy of a perfect elemental unit cell with varying lattice parameters.
www.advancedsciencenews.com www.pss-b.com (here we chose(Ge 1Àx Sn x ) 0.98 Si 0.02 and (Ge 1Àx Sn x ) 0.96 Si 0.04 alloys in comparison to Ge 1Àx Sn x ) and the three cuts from the corners of the ternary phase diagram through its center, consisting of (Ge 0.5 Sn 0.5 ) 1Àx Si x , (Ge 0.5 Si 0.5 ) 1Àx Sn x , and (Si 0.5 Sn 0.5 ) 1Àx Ge x alloys. Looking at the parity plot for the 75 ternary PBE bond energies in Figure 10d, we find again excellent linearity between the bondenergy fit (which was only performed for the binaries but included none of the ternary allots) and the DFT values with a dependence of E fit ðE DFT Þ ¼ ð0.71 AE 0.71Þ þ ð1.004 AE 0.003ÞE DFT Þ, with an only marginally larger 95% confidence interval than for the binary alloys. Similarly, the bond-energy derived mixing enthalpies for low-Si and central composition lines in Figure 10e,f look like the best fits of the DFT data. Thus, it looks like the bond energies in the ternary alloys stay basically unchanged as compared to the binary pairs, and we do not find any synergistic effects that would strengthen or weaken specific bonds when moving from the binaries to the ternary alloy as, for example, observed in oxides, [52] including for concentrated alloys around the center of the phase diagram, further confirming our previous conclusions. With that, we can use the bond energy values from Table 2 to closely approximate the mixing enthalpy of the ternary alloys per atom by a rather simple expression  Table 2. Fitted LDA and PBE bond energies for the three binary alloy systems Si-Sn, Ge-Sn, and Ge-Si, along with the average of the homobond energies (indicated by · · · h i) and the difference between the heterobond energy and the average of the homobond energies (ΔE), which are equal to the maximum enthalpies of mixing seen in Figure 4. All values in eV.
As before, based on the large enthalpy penalty for Si-Si bonds that results in the high mixing enthalpy for this system, a ternary alloy with higher Si concentrations should be thermodynamically unfavorable.

Free Energies
While zero-temperature results such as described up to here can give hints about alloy stability, for a complete description, free energy contributions from lattice vibrations and configurational entropy need to be included to determine the full Gibbs energy and assess the alloy energetics at finite temperatures. To do that, we first determine vibrational entropy contributions within the harmonic approximation [37] from where gðωÞ is the phonon density of states for frequency ω.
The configurational entropy per atom in a multicomponent alloy with a composition vector fxg can be approximated by the Stirling equation [41] F config ðfxg, assuming an ideal mixture. The phonon density of states and integration are done for the SQS cells from Γ-point phonons using the quasiharmonic approximation (QHA) approach from Ref. [37]. The resulting free energy as a function of temperature has the usual shape as shown in the example of Ge 30 Sn 34 in Figure 11. To have a parameterized model for the entire composition and temperature range, we fit the vibrational free energy data points for the different compositions with a fifth-order polynomial in temperature whose coefficients f i ðfxgÞ are second-order polynomials of the GeSnSi composition vector fxg Since Sn in the diamond structure is only stable up to 286 K and Sn overall has a melting temperature of 505 K, we have calculated vibrational free energies up to only 400 K to stay within the physically sensible range of the harmonic approximation. Vibrational free energies were calculated for the same binary and ternary compositions as we used for the mixing enthalpies and then were fitted with the expressions from Eqs. (12) and (13).
The resulting parameters are shown in Table 3 and the parity plot for the modeled vibrational free energy contributions is shown in Figure 12a.
To determine the full free energy, one adds the configurational and vibrational free-energy contributions from Eqs. (11) and (12) to the enthalpy of the system, which in turn is calculated by adding the enthalpy of mixing from Eq. (7) to the enthalpy of the ideal solid solution from Eq. (3). Skipping the enthalpy of the ideal solution will result in the free energy of mixing.
Figures 12b-d shows the resulting free energies of mixing for Ge 1Àx Sn x , Ge 1Àx Si x , and Si 1Àx Sn x for temperatures of 0, 273, and 400 K (the zero-temperature curves include the zero-point energy). Even for the relatively moderate temperature of 400 K, we find the free energy of mixing to become negative for all concentrations in the GeSi system, which is in agreement with the experimental phase diagram where these two elements form a solid solution over the entire composition range. This does not happen for GeSn and SiSn, whose phase diagrams show indeed little and no solid solubility, respectively.
To examine the ternary system, Figure 13 shows the free energy and free energy of mixing over the entire composition range. As expected from the additive nature of the binary contributions which we found for the enthalpy of mixing ( Figure 5) and for the entropy contributions ( Figure 12a for the vibrational contribution; the configurational entropy as defined in Eq. (11) is already additive), free energy and free energy of mixing are well-behaved. The free energy of mixing (Figure 13b) is lowest along the SiGe edge, followed by the values along the GeSn edge, with high values along the SiSn edge. There is especially no local minimum around the equiatomic center of the phase diagram. Thus, it is evident that when starting from an equiatomic mixture of Ge 1=3 Sn 1=3 Si 1=3 , the system should decompose into two phases, one being Si 1Àx Ge x , the other one Ge 1Ày Sn y . Once the Ge concentration x in phase 1 is chosen, y and the phase fractions f 1 and f 2 can be determined from the lever rule and conservation of mass Figure 11. Free energy versus temperature for a Ge 30 Sn 34 random alloy with vibrational free energy contribution from the quasiharmonic approximation using the QHA approach [37] and with configurational entropy contribution from the Stirling formula. Dots are DFT results, the solid line is a fit with the fifth-order polynomial from Eq. (12).
www.advancedsciencenews.com www.pss-b.com With that and with the parameterized free energy, we can now calculate the energy per atom for the two-phase system Ge 1Àx Si x þ Ge 1Ày Sn y over the entire possible composition range, 0 ≤ x ≤ 0.5. The result for temperatures of 0, 273, and 400 K is shown in Figure 14 relative to the energy of the equiatomic ternary solution. For all temperatures, the two-phase mixture is lower than the ternary solution, with minimal value at the high-concentration end, where phase 1, Si 0.5 Ge 0.5 constitutes 2/3 of the system, and phase two, pure Sn, 1/3 of the system. Table 3. Parameters for the vibrational free-energy contribution from Eqs. (12) and (13), fitted to QHA DFT data from 0 to 400 K. Produces final energy units in eV for temperatures in K. Composition indices i for the coefficients a i are 1 ¼ Ge, 2 ¼ Sn, and 3 ¼ Si.  Figure 12. Free energy results. a) Parity plot between DFT-QHA vibrational free energies and model fit using Eqs. (12) and (13). The other panels show free energies for the: b) Ge 1Àx Sn x , c) Ge 1Àx Si x , and c) Si 1Àx Sn x binary systems.
www.advancedsciencenews.com www.pss-b.com Since the energy spread of all three curves over the composition range is with ∼16 meV rather small, two-phase systems with various compositions should be possible to result from solidification. However, the energy gain over the ternary solution is rather significant and ∼3 times the thermal energy of k B T at both 273 and 400 K, indicating the difficulty of non-equilibrium stabilization of the ternary solid solution. Figure 5d shows the three sections through the ternary phase diagram from the corners to the opposing midpoints. As they have to, they intersect at the equiatomic center of the phase diagram at a concentration of 33% of the respective corner element. In this figure, the work required to reach the center of the phase diagram can be estimated by integrating the energy from 100% concentration to the equiatomic concentration at 33%. It is evident that the least work required is for pure Ge as the starting point (green line), since it is everywhere below the other two curves. This is corroborated by looking at the ternary phase diagram in Figure 13b, where the path from the Ge corner passes through more of the "blue," lower energy area than those from the other corners.

Debye Temperature Results for Si-Ge-Sn
A somewhat simpler, alternative method to approximate the vibrational free energy of the SiGeSn alloys is through the Debye model. As we will showin the following, this approach is also a validation for the vibrational entropies reported in the previous section, since the Debye temperatures are calculated from the same phonon calculations as the vibrational entropy contributions.
To use this approximation, we first determine the composition-dependent Debye temperature by calculating the lattice contribution to the heat capacity, C V , from the phonon DOS of the different alloy SQS cells [53] and then its equivalent from the Debye model, [53,54] where the Debye temperature Θ D is numerically adjusted until C D V ¼ C V . This is done for a series of temperatures, until Θ D becomes temperature independent for high enough temperatures, which is then the reported value.
The Debye temperature allows approximating the vibrational entropy contribution F D vib , which we took from Ref. [55] as The Debye function [56] DðxÞ ¼ 3 in Eq. (17)    www.advancedsciencenews.com www.pss-b.com with c 1 ¼ À0.07713, c 2 ¼ 0.02433, c 3 ¼ 0.30548, c 4 ¼ 0.06513, c 5 ¼ 0.03487, and c 6 ¼ À0.00324 as proposed in Ref. [57]. After calculating the Debye temperatures for our binary and ternary SQS systems, we find that they can be well modeled by a simple linear combination of the elemental Debye temperatures Our calculated Debye temperatures can be excellently represented by this expression as shown in Figure 15 with fitted Debye temperatures Θ Ge D ¼ 380 K, Θ Sn D ¼ 253 K, and Θ Si D ¼ 673 K, which compare favorably to the experimental values of 374, 254, and 645 K, respectively. [58] The parity plots for the contributions to the free energy from the vibrational entropy for 0 and 400 K are shown in Figure 15b,c. While both plots show satisfactory correlation with the values from the harmonic approximation, both have an offset to higher values, which we determine by regression to be À3 and À26 meV for 0 and 400 K, respectively. Despite this offset, the resulting energetics show similar trends as found within the full harmonic calculations as can be seen by the Debye-model free energy of mixing in Figure 13c in comparison to the harmonic-approximation results in Figure 13b. Despite slightly different values and contour lines, the trend is rather similar and once again points towards the binary SiGe and GeSn edges as the most probable composition regimes for ternary alloys.

Conclusions
In summary, we have introduced different modeling approximations to study enthalpy and free energy in GeSnSi alloys over the entire composition range and have fully parameterized them with the help of DFT calculations. These models include specifically a regular solution model with concentration-dependent regular-solution parameter and the quasichemical model for the enthalpy, and the harmonic approximation versus the much simpler, but somewhat less accurate Debye model for the vibrational free energy, while the usual Stirling model was employed for the configurational entropy. With that, our work provides a full parameter set to determine the free energy that can be used within the CALPHAD approach to thermodynamics [59] for further studies or as simple, fast additions for DFT-based modeling such as the MC simulations from Ref. [28] or Multi-Cell Monte Carlo simulations [60][61][62][63][64] in case full free energies without additional overhead were desired.
We have used these models to study enthalpy and free energy in GeSnSi alloys over the entire composition range to explore if concentrated alloys in the equiatomic concentration range are thermodynamically feasible. We find that the energetics are very well behaved and follow the expressions of a ternary regular solution rather closely. We find a small asymmetry in the ideal solution energetics for the three binaries but also could show that assuming a symmetric ideal solution as underlying the quasichemical model does not result in significantly different conclusions.
The enthalpy of mixing is found to have rather large values for the SiSn system, which is maintained nearly unchanged in the ternary alloy and which makes a concentrated ternary alloy energetically unfavorable. [27] Calculating the full free energy by adding configurational and vibrational entropy contributions to the enthalpy does not change the overall situation within the accessible temperature range, which is limited by the phase transition in Sn at 286 K and its melting point at 505 K.