High-Throughput Investigation of the Electron Transport Properties in Si₁-ₓGeₓ Alloys

Si<sub>1–<italic>x</italic></sub>Ge<sub><italic>x</italic></sub> alloys are among the most used materials for power electronics and quantum technology. In most engineering models the parameters used to simulate the material and its electronic transport properties are derived from experimental results using simple semiempirical approaches. In this paper, we present a high-throughput study of the electron transport properties in Si<sub>1–<italic>x</italic></sub>Ge<sub><italic>x</italic></sub> alloys, based on the combination of atomistic first principles calculations and statistical analysis. Our results clarify the effects of the Ge concentration and of disorder on the properties of the Si<sub>1–<italic>x</italic></sub>Ge<sub><italic>x</italic></sub> alloy. We discuss the results in comparison with existing semiempirical methods and we provide a Ge-dependent set of transport parameters that can be used in device modeling.

The versatility of Si 1−x Ge x relies on the tunable electronic structure and transport properties, by modulating the Ge concentration x. On the theoretical side, the most of existing works adopted semiempirical approaches for the electronic structure characterization (such as empiric pseudopotentials [25]- [27] tight-binding [28], [29] or k·p [30]- [32]) and Montecarlo methods [33], [34] for the electron transport characteristics. In those cases, the effect of the Ge inclusion within the Si host is either parametrized from experimental data or deduced using the virtual crystal approximation (VCA) [35], [36]. However, both approaches have severe limitations: the former lacks of predictability when experimental data are not available for a specific alloy concentration, the latter neglects the effects of local disorder and local deformation potentials. First-principles approaches based on density functional theory (DFT) can help overcoming these limitations, provided that a reliable atomic structure is available. In addition, most of existing DFT studies focused on the low-doping regime, where Ge is assumed as a diluted (i.e. isolated) impurity in the ideal Si crystalline host and whose effects can be described as a perturbation of the pristine system. However, relevant technological applications tend to use Si 1−x Ge x alloys with high Ge content (ca. 20-80%). For VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ this class of systems, the diluted model does not hold and new theoretical approached are necessary to take explicitly into account the effects of the disorder and impurity-impurity interaction. Albeit isoelectronic of Si, the larger size of Ge affects the lattice constant of the host, thus originating possible sources of internal strain or at the interface between the active channel and the rest of the device (e.g. source and drain in transistors). The increase of the structural disorder due to the Ge-alloying could either degrade the crystallinity of the system through the formation of Ge-Ge aggregates (i.e. nanoclustering), or favoring the localization of electronic states in the band gap region (i.e. electron traps) that may affect the transport properties of the system. Unfortunately, disorder complicates the theoretical analysis because it is not possible to represent large extended systems, while average periodic models are not sufficiently realistic.
To unravel this problem, we adopted a high-throughput approach that combines atomistic first principles calculations and statistical analysis. This allows us to gain an unbiased characterization of the microscopic effects of Ge on the structural, electronic and transport properties of Si 1−x Ge x alloys. We distinguish between the effects of the compositional concentration and the effects of disorder. Concentration influences the mean properties of the alloys (e.g. lattice parameter), while disorder imparts local features on the electronic structure (e.g. band splitting at critical high symmetry points) that alter the carrier mobility of the system. Our results support the validity of the semiempirical models used in the past and provide a Ge-dependent set of transport parameters (such conductivity, effective mass, Luttinger parameters) that can be used in device modeling.

II. COMPUTATIONAL DETAILS
The structural and electronic properties of Si 1−x Ge x compounds are calculated by using a first principles totalenergy-and-forces approach based on DFT, as implemented in QUANTUM-ESPRESSO (QE) [37]. The Perdew, Burke, and Ernzerhof (PBE) [38] generalized gradient approximation is adopted for the parametrization of the exchangecorrelation (XC) functional. Ultrasoft pseudopotentials of the Vanderbilt's type [39] are used to described the atomic potentials of Si and Ge species. Single particle orbitals (charge) are expanded in plane waves up to a kinetic energy cutoff of 40 Ry (400 Ry), respectively. Uniform (4 × 4 × 4) and (6 × 6×6) k-point grids are used for summations over the Brillouin Zone (BZ) of the (2 × 2 × 2) cubic supercells during lattice optimization and electronic structure, respectively. Spin-orbit coupling (SOC) is included by using the non-collinear fullyrelativistic approach implemented in QE [40].
Electron transport properties are evaluated by solving the Boltzmann semiclassical equation, within the scattering time approximation [41]. Electron conductivity σ ij is calculated using the PAOFLOW package [42], which relies on a minimal-space tight-binding (TB) hamiltonian representation, resulting from the projection of the original DFT wavefunctions onto a pseudoatomic orbital (PAO) basis set [43].
Effective mass tensors are calculated through the evaluation of the curvature of the band energy dispersion E n (k) at selected critical points [44], [45]. This is done by considering the full mathematical complexity of the overall 3D band structure, and not simply along the 2D band structure plots [45]. This approach adopts a double layer highthroughput engine for (i) the generation of the TB-PAO Hamiltonian used to interpolate the DFT band structure on high-dense k-point meshes, and (ii) for extraction of the band-to-band effective mass tensor at critical points. Highthroughput ground state and transport calculations have been run by using the automatic workflows implemented in the AFLOWπ infrastructure [46].

A. SYSTEM GENERATION AND STRUCTURAL ANALYSIS
A (2 × 2 × 2) cubic supercell with 64 Si atoms in the diamond structure was assumed as the initial structure for modeling the Si 1−x Ge x alloys (Fig. 1a). The cell parameter of the reference Si crystal is a 0 = 2 × a Si 0 = 10.94 Å, where a Si 0 = 5.47 Å is the lattice constant of Si bulk, obtained through the geometry optimization of the primitive fcc cell. High symmetry directions of cubic lattice are aligned to cartesian axes, as shown in Fig. 1a. Different Ge concentrations were considered for the Si 1−x Ge x alloys, with x ∈ [0, 1], where x = 0 (x = 1) corresponds to pure Si (Ge). The initial structures were prepared substituting n = 3, 8, 16, 19, 32, 45, and 48 Si atoms per cell with an equivalent number of Ge atoms (Fig. 1a). This corresponds to x = 0.04, 0.12, 0.25, 0.30, 0.50, 0.70, and 0.75, i.e. Ge concentration equal to 4, 12, 25, 30, 50, 70, and 75%.
As the number of Ge atoms increases, a single structure is not sufficient to model the Si 1−x Ge x alloy. Since Ge atoms may arrange in several nonequivalent configurations within the host, it is necessary to consider all the possible positions of the impurities within the cell and their mutual interactions. The total number of configurations to simulate rapidly diverges to an unaffordable number of systems. In order to have a statistically meaningful but computationally accessible sampling of the possible atomic distribution of Si 1−x Ge x , for each Ge concentration x we prepared 100 initial random configurations. The initial position of the n Ge-substituted sites in the supercell is defined by a string of n randomly generated numbers, ranging from 1 to 64. In principle, this random choice distribution does not span all possible configurations nor does it assure that those 100 configurations are inequivalent (because of the periodic boundary conditions), but it provides an unbiased way to generate statistically reliable structures.
For all the configurations at each concentration (i.e. 100 × 7 systems), we performed cell optimization through Murnaghan equation of state allowing for the complete relaxation of the atomic positions within the cell of all the structures. Results are summarized in Fig. 1 and Table 1. For all Ge concentrations, the distribution of the optimized lattice parameters has a Gaussian-like character centered around the mean valuesā 0 , as shown in Figure 1b for the case x = 12%. In the picture a Gaussian fit (red line) is superimposed to DFT data (blue line), for comparison. Although spatially limited to a few hundredths of Å, the full width at half maximum (FWHM) of the lattice distributions increases with the Ge amount, that is an indication of the disorder induced by the dopants. Nonetheless, these sharp and symmetric distributions justify the assumption that the centralā 0 value is representative of the lattice parameter for the corresponding Ge concentration.
As the Ge concentration increases, the mean lattice parameterā 0 also increases, as shown in Fig. 1c and in agreement with the reported experimental [47]- [50] and theoretical [51] data. Our theoretical results are systematically larger than experimental ones by ∼1.5%. The difference is expected and ascribable to the choice of the PBE XC functional. With respect to the empirical Vegard's law, which predicts TABLE 1. The calculated mean values of the optimized lattice parameters a 0 and lattice mismatch a 0 with respect to Si crystal. Energy bandgap E g and split-off gap so evaluated for representative structures at a 0 . a linear lattice expansion as a function of x, our results exhibit a deviation from perfect linearity especially in the intermediate doping range 25% ≤ x ≤ 75%. Both positive and negative discrepancies from the Vegard's law have been extensively discussed and several non-linear corrections have been proposed to best fit the experimental data [31], [47], [52]. In the present case, the main deviations are associated to the symmetry constrain (cubic simulation cell) during the lattice optimization process. The relaxation of this constrain, i.e. slight cell deformation, would help to optimize the excess volume of the mixing for the intermediate concentrations.
The overall increase of the lattice constant implies an increase of lattice mismatch a 0 with respect to the ideal Si. Within a device this may lead to a surface strain at the interface between the Si 1−x Ge x layers and the remaining parts of the device, such as a doped Si substrate. Another relevant case is represented by fully depleted silicon on insulator (FDSOI) p-MOSFET where both the source/drain and the channel are often made of Si 1−x Ge x with different Ge concentrations [9]. Interface strain is also relevant for systems like Si 1−x Ge x /Ge/Si 1−x Ge x quantum wells, where the biaxial strain induced by the interface mismatch is exploited to manipulated the electronic levels and reduce the interband scattering [16], [32].
While the volume expansion of the SiGe can be correctly predicted by VCA [53], [54], this approximation does not account for the internal distortions or the modifications of the atomic coordination induced by the Ge inclusion. Fig. 1d shows the calculated radial distribution function (RDF or g(r)) for the Si-Si and Ge-Ge distances in all the optimized geometries at each concentration. This analysis takes into account the internal atomic relaxation and gives a statistical insight into the local deviations from the crystalline symmetry of the pristine Si host. Dashed black (red) vertical lines mark the lowest Si-Si (Ge-Ge) distances in the ideal Si (Ge) fcc crystals, respectively. For all the concentrations, the main features of the crystalline system are clearly recognizable, proving that all compounds maintain an overall zincblende structure and neither phase change nor amorphization take place.
The inclusion of Ge in the Si host progressively shifts the RDF features typical of the Si crystal closer to those of the Ge one. The Ge-Ge distances (2.37-2.39 Å) from  [49]. The average elongation of both Si-Si and Ge-Ge bonds and the absence of new peaks in the RDF indicate that all structures conserved the crystalline phase features with local disorder around the ideal lattice site to accommodate Ge in the host. This is proved also by the slight but progressive broadening of the RDF peaks as the Ge content increases. These features, left out by VCA approaches, are responsible for so-called alloy scattering that affect the carrier mobility in SiGe alloys [25]. Our results rule out the formation of Ge nanoclusters within the Si matrix, since this would be characterized by the coexistence, even at low Ge concentrations, of both bulk like Si-Si and Ge-Ge distances. The coordination number of both Si and Ge results to be independent from the composition, which means that Si and Ge atoms are randomly distributed on the sites of the cubic lattice and are fully miscible at all compositions, in agreement with the experimental findings [55].
We can conclude that main effect on the structural properties is to be ascribed to the amount of Ge concentration x, while the local spatial distributions of the Ge atoms in the cell (i.e. structural disorder) has a minor impact on the overall structure of the SiGe alloy. This supports, a posteriori, the experimental practice of dealing with SiGe as a homogeneous system, even without the atomistic spatial control of the Ge species during the growth, and justifies the validity of the VCA approaches in predicting the structure of Si 1−x Ge x alloys.

B. ELECTRONIC PROPERTIES
For each Ge concentration and each structure configuration, we calculated the electronic ground state and the corresponding band structure. At fixed concentration x, the internal distribution of the Ge atoms in the supercell does not significantly modify the overall electronic properties of the alloy. This is consistent with the restricted variability of both the lattice parameter and the interatomic bond length. Thus, on the basis of the structural analysis of the previous section, for each concentration we assumed a Si 1−x Ge x structural model corresponding to the mean value of the optimized lattice parameters a 0 as representative of the SiGe alloy at fixed Ge content. Fig. 2(a) shows the band structure of Si 1−x Ge x for the cases x = 12% (left panel), and x = 75% (right panel), taken as examples of medium and high Ge concentration, respectively. The other systems have similar trends. We restricted the plot around , along the high symmetry directions (R--X) of the (2 × 2 × 2) cubic cell, taking into account the effect of the spin-orbit coupling (SOC). With respect to the well-known band structure of the Si crystal [56], we can detect three main aspects related to the increasing inclusion of Ge atoms in the Si host. The first effect is associated to the disorder: the presence of Ge in the supercell locally breaks the translation invariance of the original diamond structure. This imparts a splitting in the band degeneracy at the edge of the Brillouin zone, especially at the R point for the highest valence bands and at the point for the lowest conduction bands. The details of the band splitting depend on the specific Ge distribution of the single configurations simulated at each x concentration. It is important to note that, even though these differences are small in absolute value (1-10 meV) and may be irrelevant for transport in standard devices (e.g. MOSFET), they may be relevant for properties related to the topological character of bands or for systems (e.g. strained quantum dots [32]) that exploits band slitting to modulate the particle population within the system.
The second relevant aspect is the orbital composition of the bands. The inclusion of Ge does not insert any defect state in the band structure. Rather Ge atoms contribute to both valence and conduction bands forming bonding and antibonding sp 3 states similar to silicon, and whose ratio with respect to the total DOS depends on the actual concentration x, as shown in Fig. 2(b), for the case x = 12%. Increasing x the band structure of Si evolves into the Ge one. This involves an energy redistribution of the lowest conduction bands and the well-known differences in the band gaps of Si and Ge bulk crystals, between 25 − X 1 and + 8 − L + 6 high symmetry points of the fcc BZ, respectively [57]. As a consequence, we remark a progressive reduction of the band gap (E g ) as the Ge concentration increases. The numerical values of E g for the systems under evaluation are reported in Table 1 and shown in panel 2(c). The present E g values are systematically underestimated with respect to the experimental data. This is due to the generalized gradient approximation used to treat the exchange correlation in the present work. This can be easily corrected, e.g., by using hybrid functionals or many-body approaches [31]. However, this does not change the presented trends or the analysis presented below.
The third relevant aspect is a genuine quantum mechanical effect due to the different spin-orbit coupling of Si and Ge species. SOC is responsible for the splitting of the three double degenerate bands at the top of valence band at into the heavy-hole (hh), light-hole (lh) and split-off (so) bands (Fig. 2). The effect of spin-orbit coupling is higher for Ge due to its heavier core. Thus, the overall band split increases as the Ge concentration increases. This is particularly evident for the split-off gap so between the hh/lh and the so bands at (Fig. 2 and Table 1). Our values agree well with experimental data on pure Si ( so = 44 meV) and Ge ( so = 289 meV) [58]. Interestingly, the symmetry reduction due to the disorder within the cell imparts a further split of the hh/lh bands at point, which increases with the Ge concentration up to ∼18 meV in the case x = 75%. We remark that this band splitting and its modulation cannot be described by semiempirical approaches, unless applying distortion to the lattice (e.g. axial strain) [31]. The combination of SOC and disorder causes a significant change in the parabolic-like shape of the three main bands (especially lh) near . Such distortions are not isotropic in the k-space, being more pronounced, for instance, along the R− than the −X direction (Fig. 2). These are responsible for the anisotropy in the hole conduction measured along different direction in Ge quantum wells.

C. ELECTRON TRANSPORT PROPERTIES
The electrical conductivity tensor σ ij is evaluated by solving the Boltzmann equation in the scattering-time approximation, as [41]: where τ is the constant relaxation time, v i n (k) is the i th component of the electron velocityv calculated for the n th band for each k point in the BZ, f 0 (T ) is the equilibrium distribution function at the temperature T , and ε is the electron energy. The evaluation of (1) requires an accurate integration over a fine grid of k points in the BZ. Here, we used a tight-binding like representation obtained from a PAO projection procedure of the DFT electronic wavefunctions [42]. The group velocity is calculated from the expectation value of the momentum operatorp [43]: where m 0 is the free electron mass and |ψ n (k) = exp(−ik · r)|u n (k) are the Bloch functions resulting from DFT calculations. The gradient of the Hamiltonian operator that enters in (2) is easily evaluated in terms of the ab initio TB hamiltonian H (r ) obtained projecting the Bloch wavefunctions onto a PAO basis set [59]: Here, we focused on the hole transport in the energy region close to the top of the valence band. For each Ge concentration, we calculated the conductivity for the representative structures analyzed above and corresponding to the mean value of the optimized lattice parameters a 0 . Fig. 3 displays the diagonal terms (σ ii ) of the conductivity tensor, with i = {x, y, z} and calculated at T = 300K. The relaxation time τ for each concentration has been evaluated through a linear fit of the limiting cases τ = 1 × 10 −12 s and τ = 4 × 10 −11 s for Si and Ge, respectively, at T = 300K, from Ref. [60], which include several scattering interactions (e.g. acoustic intravalley, optical intravalley, intervalley scattering). The three diagonal components of the conductivity tensor are almost identical, while the off-diagonal terms (σ ij ) are negligible for all the Ge concentrations, being four orders of magnitude lower than the diagonal ones. This confirms that the bulk alloys maintain the main symmetry properties of cubic systems, despite the internal structural disorder. Fig. 3 shows the evolutions of transport properties from heavy p-doped (µ = −0.5 eV) to intrinsic semiconductor systems. Conductivity for µ > 0 stems from the thermal broadening. This contribution is larger for high Ge content, as the band gap is smaller (Fig. 2c). The increase of Ge concentration also imparts a systematic increase of the conductivity within the overall range of µ [with µ(T = 0K ) = E F ], in agreement with the general statement that Ge or Ge-rich alloys have lower resistivity (i.e. higher conductivity) than Si or Si-rich compounds [33]. A direct comparison with experiments is not trivial: available data for electrical conductivity exhibit a huge variability related to specific the doping level, presence of defects, temperature, and experimental setups. Being σ ij an integrated quantity over the BZ, the differences due to the k-directionality, disorder, and SOC discussed for the band structures, become less relevant in the case of the conductivity. Indeed, all curves in Fig. 3 have the same functional trend, even close to the top of valence band (µ = 0.0 eV) where bands are mainly different.
To gain more insight on the details of the electronic structure, we considered the analysis of the effective mass m * . For each band n, the mass tensor M * is associated to the curvature of the energy dispersion E n (k), through the relation: where k 0 is a generic critical (e.g. minimum, maximum, saddle) point in the entire Brillouin zone, i.e. a k-point that satisfies the condition In principle, E n (k) is a mathematical function with domain in three dimensions. Very often, the complexity of E n (k) is overlooked and the effective mass is calculated as simple parabolic fit of the 2D band-plots (as those in Fig. 2a) along a few high symmetry directions. This approach would disregard the tensorial nature of M * [44].
We compute the Hessian matrix using the second-order Fourier derivatives [61] at the valence band maxima: The final expression (4) is obtained by inverting the Hessian matrix and including the appropriate constant prefactor. The validity of the Fourier approach holds only for nonwarped bands, for which the E n (k) is second-order differentiable and the Hessian matrix is symmetric. In the presence of warped bands, as in the case of ideal fcc bulk semiconductors (e.g. C, Si, Ge [62]), other approaches have to be considered [44]. However, since warping is strictly related to analytical band degeneracies, the presence of Ge structural disorder, which relaxes the fcc symmetry constrains in the cubic supercell, prevents the appearance of non-analytic points in the Si 1−x Ge x alloys. Thus, the method described above (Eqs. 4-6) can be used for all the systems considered in this work.
The effective mass concept is largely used in several electronic structure (e.g. k·p) and transport models for two/three terminal devices, quantum-well and superlattices. Depending on the specific application, the effective mass is often expressed in different forms [63]. For systems where directionality is not controlled or for polycrystalline systems the relevant quantity is the density of state mass m * DOS . In the case of a single valley, m * DOS gives a geometric mean value of the effective mass components. Starting from the effective mass tensor M * , for each band the m * DOS can be obtained as: where m 1 , m 2 , m 3 are the eigenvalues of the effective mass tensor.
In other cases, such as nano-sized channels in transistors [64], quantum wells or strained multilayers [65], it is useful to discriminate between longitudinal (in-plane) and transverse (out-of-plane) component of the effective mass, namely m * L and m * T [33]. Directional effective masses result from matrix-vector product: where i = T , L and The angle θ is defined with respect to the cartesian axis x, oriented along the [100] direction of the cubic crystal lattice, as in Fig. 1a. The final longitudinal mass is evaluated as the average over the angle θ with θ ∈ [0 − 2π].
Finally, the hole effective masses are related to the Luttinger parameters γ 1 , γ 2 , γ 3 within the k·p formalism through the relations [30]: We focused on the hole effective masses close to the top of the valence band. Thus, for each configuration and each Ge concentration, we calculated the effective mass tensors for the hh, lh, and so bands at the point. The different formulation of the effective masses and the Luttinger parameters were first calculated using Eqs. (7)-(11) for any configuration and then averaged over the statistical sample. The results, summarized in Table 2   results. This confirms that the SiGe alloys act on average as an effective cubic system [66]. We expect that larger differences between longitudinal and traverse directions can be detected for electron effective masses, due to the multi-valley structure of the conduction band [57], or for biaxial strained systems. (iv) The Luttinger parameters follow the opposite trend, by increasing along with the Ge increment. Calculated values well agree with previous results extrapolated from cyclotron resonance measurements for Si [68] and Ge [69], as well as with theoretical data from k·p Si 1−x Ge x models [70]. We note, however, that the γ 3 values calculated here are slightly lower than the ones calculated from semiempirical approaches [26], [31]. In the k·p model, γ 3 is closely related to the hh/lh mixing at the point. In the previous works the structural disorder is not considered and the hh, lh bands are degenerate at , unless (bi)axial strain is imposed to the system. Here, the effect of disorder is intrinsically taken into account and the band splitting is a natural consequence of the local symmetry break, as shown in Fig. 2. Thus, the lower value of γ 3 reflects the lower mixing of the hh, lh bands induced by the inclusion of Ge.
In the k·p approach, the calculation of the split-off hole mass (m * so ) −1 would be particularly critical: since it requires the evaluation of the additional parameter E p ∝ | S|p x |P | 2 that takes into account the effect of remote bands through the momentum matrix element between the s-like conduction bands (S) and p-like valence bands (P) [30].
In the case of SiGe alloys E p is hard to extrapolate from experimental data and the role of remote bands can not be estimated on the basis of simple symmetry considerations (as it happens for ideal fcc crystals), because of the band shifting and splitting imparted by the Ge atoms (Fig. 2a). On the contrary, in our approach the masses for split-off band are evaluated at the same level of accuracy of the other bands, since they derive from the complete DFT electronic structure of the system.

IV. CONCLUSION
We presented a high-troughput first-principles investigation of the structural, electronic and transport properties of VOLUME 9, 2021 Si 1−x Ge x alloys. The joint atomistic and statistical nature of the proposed approach allowed us to discriminate between the overall effects due to the Ge concentration, which is responsible for the mean character of the alloy, and those related to the local disorder and the actual spatial distribution of Ge within the Si host. The lattice parameter and the conductivity are mainly affected by the former, the band structure and the effective masses can change with the latter. Our results provide a microscopic justification to the empirical approaches used so far in the literature to model the electronic and transport properties of SiGe. In addition, our approach allows to calculate and customize the parameters that enter in those models (e.g. Luttinger parameters) on the basis of the specific characteristics of the samples under investigation (concentrations, impurities, strain, etc), without any fitting of the experimental data. The capability of advanced modeling to provide materials properties ''on demand'' (i.e. targeted to a specific application) for complex (e.g. disordered) systems fosters the development of materialdevice codesign for disruptive electronics and quantum technology. He is currently a Professor with the Department of Physics, Central Michigan University. He has expertise in computational materials physics and condensed matter theory, especially the electronic structure problem in complex materials, such as thermo-electrics and piezo-electrics. He is also a Founding Member of the AFLOW consortium and applies and develops ab initio methodologies including and data analysis tools. His research interests include finding hidden relationships between coexisting and/or competing degrees of freedom in materials to unravel the origin of macroscopic functionalities.