Epitaxial Stabilisation of ${\bf \mathrm{Ge_{1-x}Sn_x}}$ Alloys

TL;DR

This study uses density functional theory and cluster expansion to analyze the thermodynamic stability of GeSn alloys, revealing volume mismatch as a key factor in immiscibility and exploring epitaxial stabilization strategies.

Contribution

It provides a detailed thermodynamic analysis of GeSn alloys and quantifies how epitaxial strain can stabilize specific compositions during growth.

Findings

Volume mismatch drives immiscibility at relevant temperatures.

Epitaxial strain reduces free energy of mixing for certain compositions.

Guidelines for substrate selection to stabilize GeSn alloys.

Abstract

The thermodynamic stability of germanium tin $\mathrm{Ge_{1-x}Sn_x}$ alloys is investigated across the composition range $0 \le x \le 1$ by applying density functional theory (DFT) together with the cluster expansion formalism (CE). It is known that GeSn alloys are immiscible and that non-equilibrium growth techniques are required to produce metastable films and nanostructures. Insight into the driving forces behind component segregation is gained by investigating the equilibrium thermodynamics of GeSn systems. The alloy free energy of mixing is computed by combining enthalpies from CE with entropy terms for configurational and vibrational degrees of freedom. Volume deformations due to the large mismatch in ionic radii are readily found to be the key driving force for immiscibility at all temperatures of relevance. This leads to a study of epitaxial stabilisation by employing latticed matched substrates to favour growth of alloys with fractional compositions of $\mathrm{x=0}$, approximately $\mathrm{x=0.5}$ and $\mathrm{x=1}$. Reduction of the free energy of mixing due to epitaxial strain in thin films is quantified for each substrate leading to indicators for growth of kinetically stable films.

Tables

Table 1. Table 1: Structural parameters of germanium and α 𝛼 \alpha -tin calculated with DFT LDA compared to reported experimental values.

	$a_0$ (Å)		$B_0$ (GPa)
	This work	Exp.	This work	Exp.
Ge	$5.64$	$5.657$Baker and Hart (1975)	$74$	$75.8$Bruner and Keyes (1961)
Sn	$6.47$	$6.489$Farrow et al. (1981); Thewlis and Davey (1954)	$47.16$	$42.5-53.1$Price et al. (1971); Buchenauer et al. (1971)

Table 2. Table 2: Formation energy decomposition into contributions listed in eq.   8 for bulk alloys at selected compositions computed from 64-atom SQS.

Composition	$\Delta E_{\mathrm{conf}}$	$\Delta E_{VD}$	$\delta E_{UR}^{chem}$	$\delta E^{int}$
(x)	(meV/atom)
0.25	36	192	-88	-67
0.50	41	207	-90	-76
0.75	28	128	-52	-48

Table 3. Table 3: Epitaxial formation energy decomposition for alloys with composition x = 0.5 x 0.5 \mathrm{x=0.5} grown on each of the substrates included in this study.

Substrate	$\Delta E_{\mathrm{conf}}^{\mathrm{epi}}$	$\Delta E_{VD}^{\mathrm{epi}}$	$\delta E_{UR}^{chem,\mathrm{epi}}$	$\delta E^{int,\mathrm{epi}}$
	(meV/atom)
Ge	-10	195	-107	-97
ZnTe/GaSb	-32	131	-89	-75
CdTe/InSb	-24	114	-76	-62

Equations

$$F(\mathrm{x})=E(\mathrm{x})-TS(\mathrm{x})$$

$$E_{\mathrm{config}}=\sum _{i,j=\{\mathrm{Ge},\mathrm{Sn}\}}{V}_{ij}{N}_{ij},$$

$$E_{\mathrm{config}}(\mathrm{x})=$$

$$\displaystyle+2\mathrm{x}(1-\mathrm{x})V_{\rm GeSn}\big{]},$$

$$V=V_{\mathrm{GeSn}}-\frac{1}{2}(V_{\mathrm{GeGe}}+V_{\mathrm{SnSn}}),$$

$$E_{\mathrm{config}}(\sigma )=\sum _{\alpha }{J}_{\alpha }{\hat{S}}_{\alpha }(\sigma ),$$

$$a_0^{\mathrm{G}{\mathrm{e}}_{1-\mathrm{x}}\mathrm{S}{\mathrm{n}}_{\mathrm{x}}}=a_0^{\mathrm{Ge}}(1-\mathrm{x})+a_0^{\mathrm{Sn}}\mathrm{x}+b_a\mathrm{x}(1-\mathrm{x}),$$

$$S_{\mathrm{conf}}(\mathrm{x})=-k_{B}N[\mathrm{x}\ln \mathrm{x}+(1-\mathrm{x})\ln (1-\mathrm{x})],$$

$$\Delta E_{\mathrm{conf}}=\Delta E_{VD}+\delta E_{UR}^{\mathrm{chem}}+\delta E^{int},$$

$$\Delta E_{VD}=$$

$$\displaystyle+\mathrm{x}\big{[}E_{\mathrm{Sn}}(a_{\mathrm{Ge_{1-x}Sn_{x}}})-E_{\mathrm{Sn}}(a_{\mathrm{Sn}})\big{]},$$

$$\delta E_{UR}^{\mathrm{chem}}=$$

$$\displaystyle-\big{[}(1-\mathrm{x})E_{\mathrm{Ge}}(a_{\mathrm{Ge_{1-x}Sn_{x}}})+\mathrm{x}E_{\mathrm{Sn}}(a_{\mathrm{Ge_{1-x}Sn_{x}}})].$$

$$\delta E^{int}=E_{\mathrm{G}{\mathrm{e}}_{1-\mathrm{x}}\mathrm{S}{\mathrm{n}}_{\mathrm{x}}}(a_{\mathrm{G}{\mathrm{e}}_{1-\mathrm{x}}\mathrm{S}{\mathrm{n}}_{\mathrm{x}}}){|}_R-$$

$$E_{\mathrm{G}{\mathrm{e}}_{1-\mathrm{x}}\mathrm{S}{\mathrm{n}}_{\mathrm{x}}}(a_{\mathrm{G}{\mathrm{e}}_{1-\mathrm{x}}\mathrm{S}{\mathrm{n}}_{\mathrm{x}}}){|}_{UR},$$

$$q(a_S,\hat{G})=\frac{\Delta E^{epi}(a_S,\hat{G})}{\Delta E^{bulk}(a_S)},$$

Full text

I Introduction

Germanium tin (GeSn) alloys have received significant attention over the past decade due to their potential applications in optoelectronic devices.Doherty et al. (2020) At low tin content these alloys have been shown to exhibit a direct energy band gap, making them particularly suitable for photonics devices that promise the possibility of being integrated into silicon-based technologies, as GeSn has been grown employing CMOS-compatible processes.Wang et al. (2007); Conley et al. (2014) Although GeSn alloys with low Sn concentrations have been reported in the literature for some time, progress towards general technological applications has been much slower mainly due to two of its key characteristics: i) Ge and Sn atoms exhibit a large difference in ionic radii which results in a solid solubility of tin in germanium of less than $1\%$,Olesinski and Abbaschian (1984) and ii) tin’s stable phase under standard conditions is the tetragonal $\beta$-Sn structure instead of the preferable cubic $\alpha$-Sn structure, which only becomes stable at temperatures below $13^{\circ}$C at standard pressure. Raynor and Smith (1958) These characteristics result in significant difficulty in incorporating Sn into a cubic Ge lattice without the precipitation of metallic $\beta$-Sn clusters. Significant progress has been made in recent years towards fabricating GeSn alloys by employing non-equilibrium epitaxial growth techniques such as molecular beam epitaxy (MBE) and chemical vapour deposition (CVD), allowing realization of thin films with thicknesses up to hundreds of nanometres. Piao et al. (1990); Gencarelli et al. (2013); Conley et al. (2014); Wirths et al. (2016); Dou et al. (2018); Wei et al. (2021); Dascalescu et al. (2020); Tai et al. (2020) Moreover, recent studies have reported lasers comprised of $\mathrm{Ge_{1-x}Sn_{x}}$ thin films with low tin content ($0.054<\mathrm{x}<0.16$), indicating significant advances towards fabrication of alloys with device-grade crystallinity. Wirths et al. (2015); Reboud et al. (2017); Rainko et al. (2019); Elbaz et al. (2020); Chrétien et al. (2019)

Although recent efforts have been focused on low tin content alloys with a direct band gap for applications in photonics, the electronic structure of germanium-tin alloys has been predicted to vary significantly with composition: while the addition of a few atomic percent of tin into germanium induces a transition from indirect to a direct band gap, alloys with larger tin contents exhibit semimetallic behaviour. Polak et al. (2017); O’Donnell et al. Integration of direct band gap $\mathrm{Ge_{1-x}Sn_{x}}$ into electronic devices promises advantages such as enhanced carrier mobilities, optical interconnects, and improved efficiency in tunnelling-based devices such as tunnelling field-effect transistors (TFETs) or tunnel diodes. Sau and Cohen (2007); Kao et al. (2012); Schulte-Braucks et al. (2015) However, the possibility of increasing tin composition to induce semimetallic regions or to counteract confinement-induced band gap widening in nanoscale devices further increases the technological value of these alloys. With critical dimensions for modern electronics already below $10$ nm, designs capable of exploiting semimetals using confinement- and surface-assisted band gap engineering are being considered, and the ability to further modify band gaps by varying alloy composition introduces another parameter to achieve new nanoelectronic device designs. Ansari et al. (2012); Sanchez-Soares and Greer (2016); Ansari et al. (2016); Sanchez-Soares et al. (2016); Gity et al. (2017, 2018) The growth of crystalline $\mathrm{Ge_{1-x}Sn_{x}}$ alloys with intermediate tin contents has also been experimentally achieved; recently, atomically flat epitaxial films grown on $\mathrm{Ge}(100)$ with tin content as high as $\mathrm{x=0.46}$ were reported in the literature. Suzuki et al. (2016) In their study, compressive strain due to high Sn content of films pseudomorphically grown on germanium substrates resulted in crystalline films with a critical thickness of 3 nm. In contrast, the use of lattice-matched substrates has been reported to allow growth of alloys with compositions $\mathrm{0.26<x<0.99}$ and even stabilisation of pure $\alpha$-Sn films with thicknesses in excess of 100 nm. Farrow et al. (1981); Hochst and Hernandez-Calderon (1983); Asom et al. (1989); John et al. (1989); Bowman et al. (1990); Piao et al. (1990) Strain reduction by growth on lattice-matched substrates thus allows fabrication of semimetallic $\mathrm{Ge_{1-x}Sn_{x}}$ with thicknesses on the order of hundreds of nanometres, which is well above the lengths required for nanoelectronic applications.

This work investigates the thermodynamic stability of bulk and epitaxial $\mathrm{Ge_{1-x}Sn_{x}}$ alloys. Density functional theory (DFT) calculations are performed in conjunction with the cluster expansion (CE) formalism to quantify the magnitude of various contributions to the Helmholtz free energy of mixing that lead to immiscibility in bulk GeSn alloys. The case of epitaxially grown alloys is then explored by analysing the impact of constraining the lattice constant of alloys along a plane on the structural properties and free energy contributions. Results suggest that high tin content alloys may be epitaxially grown by appropriate choice of substrate material and orientation.

II Methods

Computational methods are applied to study the thermodynamic stability of $\mathrm{Ge_{1-x}Sn_{x}}$ alloys for fractional compositions $0\leq\mathrm{x}\leq 1$. The Helmholtz free energy is given by

[TABLE]

with $E$ the internal energy, $T$ the temperature, and $S$ the entropy; electronic, vibrational and configurational contributions to the free energy are considered. The configurational contribution to the internal energy is considered using two approximations: the Bragg-Williams (BW) model and the cluster expansion (CE) formalism. Bragg and Williams (1934); Williams (1935a, b); Sanchez et al. (1984) Both methods define the configurational energy in terms of interaction parameters that can be extracted from density functional theory (DFT) calculations using periodic supercells. In the BW model, the bond energies do not account for local variations in the alloy composition, whereas the CE is designed to describe all bonding configurations arising in a random alloy. The DFT calculations applied here to determine the alloy interaction energies employ the usual Kohn-Sham (KS) framework in an implementation using norm-conserving pseudopotentials and linear combination of numerical atomic orbitals (NAO) basis sets.Soler et al. (2002); Ozaki (2003); Ozaki and Kino (2004); Atomistix Toolkit version 2016.4 Pseudopotentials include 4s and 4p electrons to describe germanium and 4d, 5s, and 5p electrons for tin. Basis sets used to expand wavefunctions include s4p4d3f2 NAO for germanium and s2p3d3f2 NAO for tin, where the notation indicates the number of $s$-type, $p$-type, $d$-type, and $f$-type orbitals centred about atoms of each species. Brillouin zone integrations are performed on a $k$-point grid generated according to the Monkhorst-Pack scheme with a density of at least 7 $k$-points / Å*-1*. Monkhorst and Pack (1976) Real-space quantities are discretised on a grid with a minimum cut-off energy of 100 Hartree. A parametrisation of the local density approximation (LDA) is employed for the exchange-correlation potential. Perdew and Zunger (1981)

Minimisation of the total energy with respect to ionic positions and periodic cell dimensions (geometry optimization) is performed on all supercells to determine relaxed geometries, which are defined as having a maximum force acting on an atom less than a threshold of $|\,5\times 10^{-2}\,\mathrm{eV/}$Å$\,|$ and with all elements of the corresponding stress tensors less than $|\,0.1~{}\mathrm{GPa}\,|$. The inclusion of spin-orbit interactions in the calculations is found to have a relatively small impact on the final atomic positions and total energies for a reference set of alloy supercells and thus is neglected, a finding consistent with a recent study for germanium tin alloys. Polak et al. (2017) However, it is noted that spin-orbit coupling can have a significant influence on KS band dispersions. O’Donnell et al.

The structural properties of tin and germanium in the diamond structure computed with this approach are given in table 1. The calculated equilibrium lattice parameters show agreement with experimentally reported values to within 1%. Farrow et al. (1981); Baker and Hart (1975); Thewlis and Davey (1954) The bulk modulus as calculated for germanium exhibits a deviation of less than 3% with respect to the experimental values obtained from ultrasound techniques, and the calculated bulk modulus for tin in the $\alpha$-phase lies within previously reported values obtained from neutron scattering experiments. Bruner and Keyes (1961); Price et al. (1971); Buchenauer et al. (1971)

The BW model is a mean-field model resulting from a nearest neighbour description for the configurational energy of a substitutional alloy

[TABLE]

with $V_{ij}$ and $N_{ij}$ giving the bond energies and number of bonds, respectively, between atoms of type $i$ and $j$ located on lattice sites. The BW model describes the atom-atom interactions in GeSn alloys solely in terms of the three possible nearest neighbour interaction or bond energies $V_{\rm GeGe},V_{\rm GeSn}=V_{\rm SnGe}$, and $V_{\rm SnSn}$ yielding an expression as

[TABLE]

where $z$ is the nearest neighbour coordination number and $N$ the number of atoms in a considered cell. Germanium tin alloys under typical experimental growth conditions are random alloys. Biswas et al. (2016); Mukherjee et al. (2017) Topologically, each lattice site is given an occupation Ge or Sn with probability ${\rm(1-x)}$ and ${\rm x}$, respectively, with each lattice site having four nearest neighbours as per bonding in non-amorphous group IV materials. To determine the nearest neighbour interactions in the BW model, energies for the homonuclear bonds Ge-Ge and Sn-Sn are obtained from the DFT total energies for relaxed diamond cells and for isolated Ge and Sn atoms, respectively. For the heteronuclear bond energy, a DFT total energy calculation is performed for a fictitious zinc-blende germanium tin alloy and taking the aforementioned isolated atom simulations as reference. This parametrisation neglects bond energy changes due to local variations in the occupation of nearest neighbour atom sites and from higher coordination shells. However, the BW model yields a useful indicator of alloy stability through the relative binding energy

[TABLE]

which provides a comparison of the energy of a heteronuclear bond to the average energy of the homonuclear bonds. The magnitude and sign of the relative binding energy provide a clear and intuitive indication for a preference to segregate or to form a solid solution.

The most severe limitations inherent to the BW model can be overcome by including interactions beyond nearest neighbours through the cluster expansion formalism. By characterizing lattice site occupations $i$ using spin variables corresponding to the atom type as $\sigma_{i}$, analogous to the Ising model, the configuration within a supercell representing an alloy can be represented by a vector $\sigma=\{\sigma_{i}\}$. The internal configurational energy may then expanded as

[TABLE]

whereby lattice sites are grouped into clusters denoted $\alpha$ with associated pseudospin $\hat{S}_{\alpha}$, which is the average product of lattice site spins over each set of unique clusters, and $J_{\alpha}$ are effective cluster interactions (ECI). By including all clusters, the internal energy of any alloy configuration can be exactly reproduced. However the main effects of the chemical bonding are generally speaking localized and the magnitudes for the ECIs are found to rapidly decrease with increasing cluster size and range, thus the internal energy of an alloy can be well approximated by truncating the sum in eq. 5 to include relatively small and short range clusters. The ECI in the truncated sum can then be extracted from explicit total energy calculations performed on a relatively small set of configurations in an approach known as the Connolly-Williams method or the structure inversion method. Connolly and Williams (1983) In order to define an optimal set of clusters to include when building the CE for the germanium tin alloys, a set of 47 crystalline test structures with a maximum of 8 atoms per supercell are generated following an algorithm designed to maximise the predictive power of the fit whilst minimising the computational effort of the required first principles calculations. (van de Walle and Ceder, 2002a) The total energies for the test structures are used to extract the ECI coefficients for two atom clusters with a maximum range of 1 nm and for three atom clusters with a maximum range of 0.52 nm. The selection of clusters to include in the CE applies a statistical measure relying on a cross-validation score and variance minimisation techniques to reproduce the supercell energies for the test structures. Stone (1974) The resulting cross validation score is 5 meV/atom, which is comparable to the error expected from the DFT energies used to determine the ECI values.

The temperature dependence of the free energy requires determination of the configurational, electronic and vibrational contributions to the entropy. The configurational entropy is estimated in a first approximation by enumerating unique atomic configurations for a random alloy at each composition. Electronic contributions to the free energy are computed from a temperature-independent, one-electron band model employing the density of states calculated from DFT. Wolverton and Zunger (1995); van de Walle (2009) These electronic band calculations are performed with a set of supercells designed to mimic the radial correlation functions of random alloys known as special quasirandom structures (SQS); Zunger et al. (1990); Hass et al. (1990) 64-atom cubic supercells ($2\times 2\times 2$ simple cubic) have been generated stochastically using a simulated annealing procedure. Van de Walle et al. (2013) The dimensions of generated SQS allow targeting the site-site correlations of clusters found to be most relevant by the CE whilst maintaining DFT calculations tractable. Vibrational contributions to the free energy are estimated using a bond stiffness versus bond length model. Within this model, force constant tensors are obtained by computing reaction forces on crystalline structures with slight atomic displacements and at varying levels of strain. van de Walle and Ceder (2002b); Fultz (2010) Three crystalline structures are chosen: i) elemental Ge; ii) elemental Sn; and again iii) a fictitious GeSn zinc-blende structure. Cubic supercells with 64 atoms are employed such that the minimum distance between displaced atoms is at least 1.2 nm. Reaction forces are computed at 10 different levels of strain and for atomic displacements limited to 0.02 nm away from ideal lattice sites. The force constant tensors are used to perform a polynomial fit relating bond lengths to bond stiffness; this parametrisation can then be used on arbitrarily large supercells to calculate phonon frequencies and determine vibrational contributions to the free energy $F_{vib}$. Bozzolo et al. (2007) Additional details on the bond stiffness vs. bond length model and comparisons of extracted vibrational properties with experimental literature are provided as Supplementary Information.

Two cluster expansions are defined: the cluster expansion of the configurational energy from a set 47 crystalline cells as described previously, and a second expansion for the temperature-dependent vibrational contributions to free energy computed from 64-atom SQS. These two expansions can then be merged to produce a set of temperature dependent ECI. In the following, it will be shown that the electronic contributions to the free energy are negligible, and as such are ignored. Once the TECI are determined, they can be implemented in semi-grand canonical ensemble Monte Carlo simulations to calculate a phase diagram. Configurational entropy contributions are readily included in the MC method as the simulation samples accessible configurations at a given temperature. In the semi-grand canonical ensemble, the number of atoms is fixed, with the energy and composition allowed to fluctuate based on an externally imposed difference in chemical potential between the two atom types as well as an externally imposed temperature. van de Walle and Asta (2002) Supercells based on a face-centred cubic primitive cell ($84\times 84\times 84$) and with a corresponding lattice constant of 26 nm and 1,185,408 atomic sites are employed for an improved description of long-range order and ensemble averages. The number of equilibration passes and the averaging passes are determined based on a target precision of $10^{-3}$ for the atomic composition of the phase. The phase diagram across the full alloy composition range is then determined using the procedure described in ref. van de Walle and Asta, 2002.

III Results and Discussion

III.1 Bulk alloys

The configurational energy of mixing or formation energy is evaluated in a first approximation using the BW model for the composition range $0\leq x\leq 1$. Computed nearest neighbour interactions result in a relative binding energy value of $V=+11~{}\mathrm{meV}$: although its positive sign indicates phase segregation is energetically favourable, its magnitude is near the precision limit of DFT simulations employed in its calculation. Thus while the BW model predicts a slight preference for the system to segregate, kinetically grown, strain-free metastable alloys could exhibit random species distribution due to the relatively large energetic cost of rearranging bond types when compared to predicted stability gains. The magnitude of the relative binding energy also suggests that entropy contributions to the free energy may have a substantial impact on reducing the the free energies for mixing tending to stabilise metastable structures at typical growth temperatures.

To provide a more accurate treatment of the free energies and to identify the dominant stabilising effects for germanium tin alloys, the CE is applied to account for configurational variations beyond the nearest neighbour description provided by the BW model. Figure 1(a) shows the extracted ECI coefficients for the CE. The magnitude of two-atom cluster ECIs are seen to decay rapidly with distance from a maximum of approximately 20 meV/atom for the nearest-neighbour pair cluster to less than 5 meV/atom for longer range pairs. The ECIs of the three-atom clusters are an order of magnitude smaller than typical values for pair clusters, with a maximum value of 0.5 meV/atom. Consistent with the simple BW analysis and in agreement with experiment, the CE results also predict germanium and tin to be immiscible at zero temperature and thermodynamic equilibrium as shown in fig. 1(b).

Results obtained with lattice models are compared directly to DFT calculations over the set of 64-atom SQS. Atomic positions and lattice constants in the SQS are relaxed using DFT by constraining the cell shape to cubic while allowing the cell volume to vary, an approach previously argued to provide total energies representative of disordered states. Ghosh et al. (2008) Figure 2 shows the evolution of Ge-Ge, Ge-Sn, and Sn-Sn bond lengths with alloy composition as determined from SQS. Bond lengths calculated for low Sn composition deviate less than 2% with respect to values reported in previous theoretical and experimental literature.O’Halloran et al. (2019); Gencarelli et al. (2015) The difference in atomic radii between Ge and Sn and the local variations in chemical environments result in bond lengths with standard deviations on the order of 1% at intermediate compositions. The alloy lattice constants are compared to literature values by fitting the optimized SQS lattice constants $a_{0}$ at varying compositions using

[TABLE]

from which a bowing parameter value of $b_{a}=0.0056$ nm is obtained, consistent with values reported in other recent theoretical and experimental works. Gencarelli et al. (2013); Xu et al. (2017); Polak et al. (2017)

Figure 3 compares the formation energy of random alloys across the full composition range as predicted with the BW model, the cluster expansion, and from DFT calculations of SQS. It is observed that although the BW model correctly predicts immiscibility, it underestimates the formation energy of GeSn random alloys by an order of magnitude with respect to results obtained from the CE fit and from 64-atom SQS. Hence there is a stronger tendency for germanium and tin to segregate than the local bond analysis from the BW model indicates. The improved treatment by the CE is seen through its ability to reproduce the energies of relaxed SQS with an RMS error of 2 meV/atom (green dashed curve).

The free energy at finite temperature has been computed by including the electronic, vibrational, and configurational contributions. Electronic contributions are computed from the density of states using the KS band structures obtained for the SQS. The vibrational contributions are calculated by the bond stiffness versus bond length model, and configurational entropy is estimated by enumerating unique configurations at a given composition for a cell containing a large number of atoms, using the Boltzmann entropy expression and Stirling’s approximation to yield

[TABLE]

with $k_{B}$ Boltzmann’s constant. We note that while this simplified approach assumes all alloy configurations at a given composition are equally accessible and thus represents an upper limit to configurational entropy, it provides estimates useful for comparing the magnitude of various contributions to free energy. Figure 4 shows the resulting temperature dependence of the total free energy of mixing and its various contributions for the representative compositions of $\mathrm{x}=0.25$, $0.50$, and $0.75$; it is seen that the calculated electronic entropy contributions remain below 1 meV/atom for temperatures up of 1000 K for all three alloy compositions. The configurational entropy contributions dominate the temperature dependence of the free energy of mixing and serve as the main stabilising force at all temperatures of interest, with vibrational contributions favouring alloy formation to a lesser degree. The critical temperature, i.e. temperatures at which the disordered phase becomes energetically favourable, is seen to decrease for larger tin compositions in a result attributable to the asymmetry observed in fig. 3 for the composition-dependent formation energy $\Delta E_{\mathrm{conf}}$; while this asymmetry is reproduced by both the CE fit and SQS calculations, the formation energy in the BW model reduces to $\Delta E_{\mathrm{conf}}=V[\mathrm{x(1-x)}]$ which is symmetric about $\mathrm{x}=0.5$.

To further investigate the origin of the asymmetry, the formation energy is partitioned as

[TABLE]

where the three contributions correspond to the volume deformation energy, the chemical or spin-flip energy, and the internal relaxation energy. The volume deformation energy $\Delta E_{VD}$ is defined as the energy required to hydrostatically strain Ge and Sn cells to the equilibrium lattice parameter corresponding to a specific alloy composition. The resulting strain energies are weighted by their respective fractional compositions in the alloy

[TABLE]

where $E_{i}(a_{j})$ is the energy for $i$= {Ge,Sn} at the equilibrium lattice constant $a_{j}$ for $j$={Ge,Sn, Ge1-xSnx}. The chemical or spin-flip energy $\delta E^{\mathrm{chem}}_{UR}$ corresponds to the energy gained when germanium and tin atoms occupying (unrelaxed–$UR$) diamond lattice sites form heteronuclear bonds

[TABLE]

The internal relaxation energy $\delta E^{int}$ is the energy gained when the atomic positions in an alloy are allowed to relax by displacing them from the ideal diamond lattice sites:

[TABLE]

where $R$ and $UR$ indicate whether the atomic positions are relaxed or constrained to the ideal diamond lattice sites, respectively. Figure 5 shows the dependence of each contribution as computed for 64-atom SQS across the composition range, and table 2 lists their magnitude at compositions x=0.25, x=0.50, and x=0.75. The large difference between the Ge and Sn equilibrium lattice parameter results in the volume deformation energy dominating and destabilising the alloys, and the marked asymmetry observed for this quantity can be ascribed to the significant difference between components’ bulk moduli; the volume expansion required to accommodate Sn into Ge-rich alloys increases formation energy significantly more than the corresponding volume contraction associated with incorporation of Ge into Sn-rich alloys. Contributions arising from chemical interactions between different atomic species and the relaxation of internal coordinates are of a similar magnitude and act to stabilise the alloys, partially offsetting the effects of the volume deformations. Although the magnitude of stabilising contributions correlates with the magnitude of volume deformation energies, the asymmetry present in the latter dominates formation energies and thus results in larger destabilisation for germanium-rich alloys compared to tin-rich alloys.

The impact of the asymmetry on the system’s critical temperature is revealed by the phase boundary between the disordered phase (i.e. random alloy) and decomposition into the elemental phases. Lattice model Monte Carlo simulations employing the CE fit enable larger supercells than tractable with DFT to be studied, resulting in a much improved treatment for ensemble averages over large numbers of alloy configurations. Temperature effects are included by combining ECIs describing the internal configurational energies with configurational and vibrational free energy contributions. The resulting phase diagram for $\mathrm{Ge_{1-x}Sn_{x}}$ is shown in fig. 6, where two horizontal lines indicating the temperatures for the $\alpha\to\beta$ phase transition and melting point of tin have been included for reference. The higher transition temperatures found with this method when compared to those shown in fig. 4 are attributed to a more accurate treatment of configurational entropy contributions, as the MC simulation accounts for energy variations across alloy configurations. Germanium-rich alloys with compositions up to $x=0.005$ are predicted to be stable below the eutectic temperature, in agreement with previous reports;Olesinski and Abbaschian (1984) for tin-rich alloys stability is only observed for compositions $x>0.995$ at temperatures where $\alpha$-Sn remains stable. Predicted critical temperatures at intermediate compositions are well above the melting temperature for tin. These results clearly indicate that growth of bulk equilibrium $\mathrm{Ge_{1-x}Sn_{x}}$ alloys at technologically relevant compositions is not possible at standard pressure. It is worth noting that inclusion of vibrational free energy contributions reduces predicted critical temperatures, with larger reductions for germanium-rich compositions. This can be attributed to a stabilising effect associated from softening of phonon modes by addition of tin, partially counteracting the destabilising effects of volume deformations. Gassenq et al. (2017); Xu et al. (2019); Bouthillier et al. (2020)

III.2 Epitaxial alloys

The primary factor leading to the high critical temperatures has been shown to arise from the large difference in equilibrium volumes for germanium and tin. Although these findings can be anticipated based on the large difference between the atomic radii of germanium and tin, the purpose of this work is to quantify the magnitude of the energies required to be overcome to stabilise thin films. By coherent epitaxial growth, the lattice spacing parallel to a substrate can to a large extent be constrained within a thin film and result in stabilisation.Wood and Zunger (1988); Zunger and Wood (1989); Ozoliņš et al. (1998a) The extent to which coherent growth can stabilise germanium tin alloys is examined by constraining lattice parameters for the SQS along a selected plane to mimic the effect of growth on substrates for which epitaxially stable films of $\alpha$-Sn and/or germanium-tin alloys have already been demonstrated, namely: $\mathrm{Ge}$ (0.564 nm), $\mathrm{ZnTe/GaSb}$ (0.610 nm), and $\mathrm{CdTe/InSb}$ (0.648 nm). Farrow et al. (1981); Hoffman et al. (1989); John et al. (1989); Bowman et al. (1990); Calavita et al. (1990); Piao et al. (1990) To explore the most promising orientations for stabilising the alloys by coherent growth, an epitaxial softening function is defined as

[TABLE]

which is the ratio between the increase in energy due to biaxial strain as a result of epitaxial growth along a direction $\hat{G}$ on a substrate with lattice constant $a_{S}$, and the increase in energy due to hydrostatically straining to the substrate lattice constant $a_{S}$. The ratio between the two energies quantifies the degree to which out-of-plane relaxation can stabilise an epitaxially grown thin film, with lower values of $q(a_{S},\hat{G})$ signalling increased stability. Lower values of the softening function in general imply a reduction in the formation of dislocations and other strain-induced film or surface defects.Ozoliņš et al. (1998b) Figure 7 shows the epitaxial softening functions of germanium and tin as computed with DFT; it is seen that $\langle 100\rangle$ is the softest direction for germanium and tin for values of $a_{S}$ in the vicinity of their equilibrium lattice constants. This conclusion is retained at all compositions if the alloy epitaxial softening function is estimated as a weighted sum of the elemental softening functions. This suggests that $\langle 100\rangle$-oriented films will tend to minimise internal strain energy and thus structural defects for epitaxially coherent $\mathrm{Ge_{1-x}Sn_{x}}$ films, and is consistent with at least one experimental analysis. Asano et al. (2014)

Epitaxial growth on (001)-oriented substrates is evaluated by constraining the lattice parameters of 64-atom SQS along $[100]$ and $[010]$ to match the lattice spacing to that of the three substrates under consideration. The third lattice parameter for each supercell is allowed to relax. This scheme neglects the role of surface and interfacial effects, but it does provides a baseline for the energetics of the alloy compositions independent of these effects. Figure 8 shows the computed epitaxial formation energy relative to decomposition into the coherently strained elemental constituents. In this case, the substrates suppress the tendency to segregate due to the coherent strain assumed for the segregated elemental components leading to a stabilisation of the alloys even at zero temperature. Note geometry optimizations for SQS biaxially strained to the Ge$(001)$ lattice constant with compositions above $\mathrm{x=0.56}$ resulted in amorphous structures and thus are not included in the analysis of Figure 8(a). Formation energies obtained from DFT calculations have been used to fit cluster expansion ECIs in order to investigate the case of random alloys and eliminate restrictions on the cluster ensemble averages arising from the use of the SQS approach. The CE predict all alloys to be energetically favourable with respect to coherent decomposition with a trend whereby stability decreases along the (ZnTe/GaSb)-(CdTe/InSb)-Ge sequence, indicating growth on germanium to be the least favourable choice for alloy growth if coherent decomposition is assumed. Additionally, fig. 8 includes the magnitude of the out-of-plane lattice parameter $c_{\mathrm{eq}}$ for each substrate as computed with DFT for SQS and a continuum elasticity theory (CET), Zunger and Wood (1989) where good agreement between both methods is found across most of the composition range in Ge and ZnTe/GaSb substrates. However, significant deviations between the two approaches are seen at low tin content alloys coherent to CdTe/InSb substrates due to anharmonic effects at these large strains not taken into account within CET. From analysis of the vibrational contributions to the free energy in a manner similar to fig. 4, it is found that the aforementioned sequence in pseudomorphic alloy stability is maintained at finite temperatures, although we note that the bond stiffness vs. bond length model and harmonic approximation employed in this work are inadequate to describe the vibrational properties of pseudomorphic germanium and tin with up to 15% biaxial strain, as evident from $c_{\mathrm{eq}}$ calculations employing CET.

Decompositions of epitaxial formation energies akin to those defined for bulk alloys in eq. 8 are presented in table 3 for alloys at a composition $\mathrm{x=0.5}$. The effects of growth on lattice-matched substrates can be observed in the formation energy decomposition corresponding to the ZnTe/GaSb substrates. Whereas the computed epitaxial spin-flip $\delta E^{\mathrm{chem},\mathrm{epi}}$ and internal relaxation $\delta E^{\mathrm{int,epi}}$ energies remain similar to the values found for bulk alloys, the volume deformation energy relative to epitaxially coherent decomposition $\Delta E_{VD}^{\mathrm{epi}}$ is reduced by $37\%$ resulting in the observed enhancement to stability. By comparing results across the different substrates, it is observed that destabilisation of component segregation results in an overall reduction of $\Delta E_{VD}^{\mathrm{epi}}$ with increasing substrate lattice parameter due to germanium’s larger bulk modulus. However, the magnitudes of stabilising contributions $\delta E_{UR}^{\mathrm{chem},\mathrm{epi}}$ and $\delta E^{\mathrm{int,epi}}$ are observed to also decrease with increasing bond lengths associated with alloys grown on substrates with larger lattice spacings. This leads to the ZnTe/GaSb substrate exhibiting the lowest value of $\Delta E_{\mathrm{conf}}^{\mathrm{epi}}$ across the three substrates.

As a final point, the relative stability of the alloys grown on different substrates in fig. 9 are compared by computing the formation energy for epitaxial alloys with respect to non-coherent decomposition into their (unstrained) bulk elemental components. The substrate providing the lowest formation energy, and thus highest stability with respect to segregation, is observed to depend on alloy composition. Within the set of substrates considered, Ge is preferred for compositions below $\mathrm{x=0.25}$, CdTe/InSb for compositions above $\mathrm{x=0.76}$, and ZnTe/GaSb for $0.25<\mathrm{x}<0.76$, indicating the latter to be energetically favourable across much of the composition range where the alloy exhibits are expected to display semimetallic behaviour. The asymmetry discussed for bulk alloys whereby tin-rich compositions have relatively lower formation energies is also present for alloys grown on ZnTe/GaSb, as seen in fig. 9.

IV Conclusion

A study for the thermodynamic stability of GeSn alloys has been undertaken with models for both relaxed and pseudomorphic thick GeSn layers. A CE fit is performed and demonstrated to be capable of reproducing energies from DFT calculations for SQS across the full alloy composition range. From the CE, an asymmetry in the formation energies of random alloys as a function of alloy composition is identified whereby introduction of Sn into Ge-rich alloys results in larger formation energies relative to the case of introducing Ge into Sn-rich alloys. A partitioning of the formation energy into a volume deformation, chemical and relaxation terms reveals the volume deformation energy dominates the zero-temperature formation energies and drives immiscibility for these alloys with the chemical and relaxation energies serving to reduce the formation energies. Evaluation of free energy at finite temperatures reveals the configurational entropy dominates stabilising contributions with vibrational terms contributing to a lesser extent, and the electronic contributions have negligible impact in the temperature range of interest. CE fits describing configurational and temperature-dependent vibrational properties are employed to calculate the system’s phase diagram using lattice Monte Carlo simulations. It is shown that the critical temperature for stability of Ge1-xSnx alloys generally decreases at higher Sn composition as a result of the aforementioned asymmetry in the formation energy versus composition. As well, vibrational terms introduce an opposing asymmetry in the phase diagram whereby the stability of Ge-rich alloys is enhanced at finite temperatures due to the softening of phonon modes associated with incorporation of Sn into the lattice. The phase diagram identifies critical temperatures between the range of the melting temperatures for Sn and Ge for essentially the entire composition range, in good agreement with known experimental properties for bulk GeSn alloys. It should be noted that relaxation in thin films or nanowires can be enhanced by strain relaxation normal to the surfaces and can serve to stabilise nanostructures.

The influence of pseudomorphic growth on stability was considered by studying alloys biaxially strained to the lattice spacing of substrates previously used for coherent epitaxial growth of Ge1-xSnx which provide lattice matching at compositions x=0, and approximately for x=0.5 and x=1. The preferred growth orientation is analysed with the epitaxial softening function which indicates the (001)-orientation is preferable at small to moderate strains. Analysis of the formation energies reveals pseudomorphic growth enhances alloy stability by impeding coherent segregation resulting in alloy stability even at zero temperature. Growth on ZnTe/GaSb substrates is predicted to provide the greatest epitaxial stabilisation across most of the composition range. Comparing pseudomorphic formation energies across substrates with respect to incoherent segregation reveals ZnTe/GaSb to be energetically favourable across approximately half of the composition range ($\mathrm{0.25<x<0.76}$), with Ge (CdTe/InSb) substrates providing greater stability for $\mathrm{x<0.25}$ ($\mathrm{x>0.76}$). This result correlates with a recent study where similar crossing points were determined when computing thin film critical thicknesses using continuum elasticity theory. O’Donnell et al. (2020) It is worth noting that even though, for example, the results for the formation energies for an x=0.25 alloy grown on either Ge or ZnTe/GaSb substrates are predicted to be comparable, their resulting electronic properties can be significantly different. Recent calculations report the compressive strain built into Ge0.75Sn0.25/Ge results in a zero band gap semimetal band structure, whereas the tensile strain built into Ge0.75Sn0.25/ZnTe leads to a band overlap (negative band gap) of approximately 600 meV. O’Donnell et al. (2020)

Results presented in this work emphasise built-in strain to be the driver of immiscibility in Ge1-xSnx alloys. Strategies for stabilisation must thus be directed at tackling this issue by balancing several complementary effects such as relaxation normal to surfaces in nanostructures, by targeting Sn-rich compositions, and by choice of substrates that induce epitaxial stabilisation of alloys while delivering suitable critical thicknesses and electronic properties. In the context of materials engineering, selection of alloy composition, strain, film thickness (quantum confinement), and surface chemistry enables a wide range of options for the tailoring of the thermomechanical and electronic properties of GeSn alloys for technological applications.

Acknowledgements.

This work was funded by Science Foundation Ireland through a Principal Investigator award Grant No. 13/IA/1956, by the Nottingham Ningbo New Materials Institute, and the National Science Foundation of China Project Code 61974079.