Composition and temperature dependence of self-diffusion in Si1−xGex alloys

The knowledge of diffusion processes in semiconducting alloys is very important both technologically and from a theoretical point of view. Here we show that, self-diffusion in Si1−xGex alloys as a function of temperature and Ge concentration can be described by the cBΩ thermodynamic model. This model connects the activation Gibbs free energy of point defects formation and migration with the elastic and expansion properties of the bulk material. The approach allows the systematic investigation of point defect thermodynamic parameters such as activation enthalpy, activation entropy and activation volume, based on the thermo-elastic properties (bulk modulus and its derivatives, mean atomic volume and thermal expansion coefficient) of the two end-members of the Si1−xGex alloy. Considerable deviations from Vegard’s law are observed, due to the diversification of the bulk properties of Si and Ge, in complete agreement with the available experimental data.

Over the past years in microelectronics there was the technological drive to replace silicon (Si) with higher mobility substrates such as silicon germanium (Si 1−x Ge x ) alloys or germanium (Ge). These materials have some common features to Si, however, their defect processes differ and were not as well established as in Si [1][2][3][4][5][6][7][8][9] . For example, although Ge is isostructural to Si its defect processes are very different and this consitutes the formation of n-type doped Ge regions problematic (high n-type dopant diffusion) 8 . Considering Si 1−x Ge x it can be described as a group IV semiconductor random alloy as effectively there is one lattice site but two atomic species that can occupy it. Therefore, in Si 1−x Ge x there is a range of local environments including Si-rich and Ge-rich regions that can influence defect processes such as the formation of dopant-defect clusters and self-diffusion 10 . From an experimental viewpoint self-diffusion in Si 1−x Ge x has been studied for numerous years [11][12][13] , whereas the increasing computational resources and the use of density functional theory (DFT) over the past years have facilitated the application of theoretical approaches to study Si 1−x Ge x 14,15 . Although DFT can provide insights into the diffusion properties of ordered materials it is more difficult to implement when considering random alloys as this will require numerous large cells. Even the use of methods such as special quasirandom structures, which can constitute most random alloy issues computational tractable will require extensive resources when considering the complete composition and temperature range for random alloys, given that ab initio molecular dynamic calculations will be required to study diffusion.
Thermodynamic approaches may bridge this gap as they can be employed in synergy to experiment or advanced computational modeling. Such a well-established thermodynamic model is the so-called cBΩ model, which has been employed during the last four decades (e.g., see ref. 16) to describe the point defect thermodynamic parameters in numerous materials, including metals, oxides, semiconductors alkali and silver halides, diamond and minerals of geophysical interest, as well as in materials that exhibit superionic conductivity at high temperatures [17][18][19][20][21][22][23][24][25][26][27] . The model is based on the theoretical justification that the activation Gibbs free energy of the formation (or migration, or activation) of a point defect is proportional to the bulk modulus of the solid material, B and its mean atomic volume, Ω, i.e., g i = c i BΩ, where i refers to the formation, migration or activation process [28][29][30] .
The significance of the cBΩ model to describe successfully the point defect thermodynamic parameters in different categories of solids has recently emerged with its implementation in semiconductors (Si, Ge, GaAs) and nuclear fuels [31][32][33][34][35][36][37] . The potential of the model can be also extended to describe self-or hetero-diffusion in alloys but the examples are rather limited and are restricted only to alkali-or silver-halides mixed crystals and mixed-oxide nuclear fuels [37][38][39][40] . In the present study, self-diffusion in Si 1−x Ge x is investigated by employing the cBΩ thermodynamic model, in conjunction with recent experimental results. Based on the bulk properties of the two end-members, various point defect thermodynamic parameters, such as activation enthalpy, activation entropy and activation volume have been calculated as a function of temperature and Ge concentration.

Results and Discussion
Determining the Thermo-elastic properties of the Si 1−x Ge x alloys. To describe self-diffusion in a binary alloy A 1−x B x (such as Si 1−x Ge x ), the cBΩ thermodynamic model can be applied by considering that, each atom of the component B that is added to the homogeneous crystal of the pure component A can be treated as a point defect 30 . The various point defect thermodynamic parameters, such as activation Gibbs free energy g act , activation enthalpy h act , activation entropy s act and activation volume υ act , can be expressed through the same relations, as in the case of a unary solid (refer to Eqs 7 and 9-11 of Methods), where the thermo-elastic properties, i.e., the bulk modulus, B, its derivatives (∂B/∂T| P and ∂B/∂P| T ) and the volume thermal expansion coefficient, β refer to the alloy and are functions of temperature and the molar concentration, x. To estimate these bulk properties for the Si 1−x Ge x alloy, from the corresponding properties of the two constituents (Si and Ge), we proceed as follows.
The composition and temperature dependence of the molar volume of the Si 1−x Ge x alloy may be expressed to a first approximation according to the following relation where x is the molar concentration of Ge in the alloy and V Si , V Ge denote the molar volumes of Si and Ge respectively. The derivation of the above equation is based on the assumption that the volume change of the alloy due to the replacement of one atom of the constituent A with an atom of constituent B is independent of the composition, x 30 . In Eq. 1 we may substitute the molar volumes V i with the mean atomic volumes Ω i of the alloy and its end-members. Considering that the lattice constant of the Si 1−x Ge x alloy deviates slightly from Vegard's law 41,42 , Eq. 1 may be used as a first approximation to estimate the mean atomic volume of the alloy, as a function of concentration and temperature (for more details refer to the Supplementary Information).
The composition (and temperature) dependence of the bulk modulus, is estimated from the pressure derivative of the molar volume of the alloy, as given by Eq. 1. Recalling that B = −V(∂P/∂V| T ): where we have replaced the ratio of the molar volumes of the two end members (V Ge /V Si ) with the ratio of their mean atomic volumes (Ω Ge /Ω Si ). This non-linear equation in x provides a direct estimation of the bulk modulus of the Si 1−x Ge x alloy, at any desired concentration and temperature, which is based solely on the bulk properties of the end members 30 . Regarding Ge, the available data of B Ge (T) and Ω Ge (T) suggest linear relations with respect to temperature, i.e., where the subscript refers to the corresponding properties at room temperature 32 . These values [43][44][45][46] are summarized in Table 1. For Si, the mean atomic volume, Ω Si (T) was estimated from the lattice parameter of the Si crystal structure, Derived from the analytical expression of β Derived from the analytical expression of β Based on the previous considerations, the 3D surface plot of the bulk modulus, as a function of temperature and Ge concentration is depicted in Fig. 1. The temperature range is restricted to 963 K-1543 K, where diffusion measurements of Si and Ge in Si 1−x Ge x alloys have been reported and are used in the present study 2,9 .
The isobaric temperature derivative of the bulk modulus of the alloy, Similarly, the isothermal pressure derivative of the bulk modulus of the alloy, Si Ge x x 1 may be estimated from the following expression 30 :  a non-linear behavior of activation enthalpy with Ge concentration, i.e., an upward bowing for both, Si and Ge diffusion in the alloy, with increasing x. In the following, we will show that these experimental findings are interpreted within the framework of the cBΩ thermodynamic model. To proceed with the implementation of the cBΩ model to the Si 1−x Ge x alloy, the reported experimental values of Si and Ge self-diffusion coefficients by Kube et al. 2,9 have been plotted as a function of the quantity, BΩ/k B T (see Fig. 2) as proposed in ref. 50. We recall that the bulk modulus B and the mean atomic volume Ω of the alloy have been estimated as a function of temperature and Ge concentration, according to Eqs 1 and 2. We observe that linear relations hold for both, Si and Ge diffusivities, at any concentration x, implying the validity of the cBΩ model (in accordance with Eq. 12 of the Method). The parameters c act have been estimated at each concentration from the slopes of the linear fittings and their values are shown in Fig. 3. A second order polynomial fitting has been applied to the derived values, in order to extrapolate the values of c act to the entire range of Ge concentration (0 ≤ x ≤ 1). Subsequently, these values are necessary to calculate the point defect thermodynamic parameters (according to Eqs 7 and 9-11 of the Method).
The activation Gibbs free energy, g Ge act of Ge self-diffusion in Si 1−x Ge x alloys, as a function of temperature and Ge concentration, is illustrated in Fig. 4a. It is observed that, g Ge act exhibits an upward bowing and a maximum value at x = 0.10-0.15, depending on the temperature. The effect of temperature to g Ge act is also significant at the entire range of Ge concentration, resulting in a variation from 2.74 to 3.03 eV for self-diffusion in Ge (x = 1.0), and from 3.44 to 3.89 eV, in the case of Ge diffusion in Si (x = 0). Obviously, the observed variation of g Ge act should have a considerable contribution to the calculated activation enthalpy.  To estimate the activation entropy, s Ge act and the activation enthalpy, h Ge act (through Eqs 9 and 10 of the Method), the isothermal pressure derivative of the bulk modulus of Si, ∂B Si /∂P| T should be determined, (refer to Eq. 5) which in turn defines the thermal expansion coefficient, β of the alloy. To the best of our knowledge, the temperature dependence of the pressure derivative of Si bulk modulus, B Si has not been determined experimentally, but it can be roughly estimated by using the Rose-Vinet universal equation of state (EoS) which states that 51, 52 where the zero subscript refers to zero (ambient) pressure, T r is a reference temperature, X = [Ω 0 (T)/Ω 0 (T r )] 1/3 and η 0 (T) = (3/2)(∂B/∂P| p=0 −1). In this way, the calculation of s Ge act and h Ge act becomes feasible over the entire temperature range and at any Ge concentration. Thus, the energy term, Ts Ge act , as well as the activation enthalpy h Ge act of Ge diffusion have been plotted with respect to T and x in Fig. 4b,c, respectively. A considerable change of the term Ts Ge act is observed for Ge diffusion in Si (x = 0), i.e., from 0.02 to 1.15 eV, while this change is less pronounced (0.46-0.75 eV) for Ge self-diffusion (x = 1). This in turn, causes a monotonic increase of Ts Ge act with increasing x at low temperatures which, however, is reversed at higher temperatures. The activation enthalpy, h Ge act , with respect to x, exhibits a maximum value (4.10 eV at x = 0.35) at low temperature (953 K) that is, however, shifted at lower concentrations, with increasing temperature. Finally, this maximum disappears as we approach high temperatures. At 1543 K, h Ge act decreases from 4.59 eV to 3.50 eV, with increasing Ge concentration, exhibiting an upward bowing. This non-linear behavior is in good agreement with the experimental values of activation enthalpy of diffusion reported by Kube et al. 2,9 . Specifically, they reported a variation from 4.83 to 3.13 eV for Ge diffusion, which has been described by a quadratic correction term in Vegard's law i.e., Q(x) = (1 − x)Q(0) + xQ(1) + x(1 − x)Θ, with Θ denoting the bowing parameter. This empirical description arises effortlessly within the framework of the cBΩ model and may be attributed to the diversification of the bulk properties of Si and Ge. Specifically, it has been recently reported 33 that the non-linear temperature dependence of activation enthalpy and activation entropy of self-diffusion in Si, which has been experimentally measured by Kube et al. 53 , can be explained in terms of the cBΩ model, by considering the non-linear anharmonic behavior of the isothermal bulk modulus of Si 48 . This peculiar behavior of diffusion in Si is clearly observed in Fig. 4c where h Ge act in Si (x = 0) varies non-linearly with temperature from 3.91 to 4.59 eV. In contrast, at high Ge content (x > 0.6), no temperature dependence of h Ge act is observed, due to the linear variation of the bulk modulus of Ge at this temperature range, in agreement with a previous study of self-diffusion in Ge 32 . Finally, the temperature and concentration dependence of the activation volume, υ Ge act is illustrated in Fig. 4d. We observe that the temperature dependence of υ Ge act is negligible all over the concentration range of Ge, while a considerable variation of υ Ge act occurs with increasing concentration. Specifically, at 963 K, υ Ge act decreases from 29.7 Å 3 (at x = 0) to 14.2 Å 3 for Ge self-diffusion (x = 1). These values correspond to (1.48 ± 0.07)Ω o,Si and (0.62 ± 0.06)Ω o,Ge , respectively, in agreement with previous calculations of self-diffusion in Si and Ge 32,33,54 . The sign and magnitude of the activation volume provides evidence on the diffusion mechanism, i.e., vacancy (V) or self-interstitial (I), through the relation υ Ω υ υ = ± + + ο  V I  act  V I  r  V I  m  , , , , where the positive sign of Ω o refers to V and the negative to I formation, υ V I r , is the relaxation volume around the point defect (V or I) and υ V I m , is the corresponding migration volume 33,55 . In the present study, the positive sign of υ Ge act and the range of the calculated values indicate that the vacancy mechanism is prevalent throughout the temperature range considered. Ignoring the negligible effect of temperature to υ Ge act , we observe that Vegard's law describes in a good approximation the activation volumes of Si 1−x Ge x alloys with Ge concentration, x, i.e., υ Overcoming limitations and prospects of the cBΩ model. In order to apply the cBΩ thermodynamic model to estimate various point defect parameters, the thermo-elastic properties of the two end members of the binary alloy should be known as a function of temperature and/or pressure. However, these bulk properties are not always known from experimental or theoretical studies. Even for "simple" well studied binary systems such the Si 1−x Ge x alloy of the present study, the derivative of the bulk modulus has not been defined experimentally and it has been approximated by using the Rose-Vinet universal EoS 51,52 . Furthermore, the temperature derivative of ∂B/∂P| T (see Eq. 5) may be difficult to find in the literature, however, it can be estimated via the approximation, ∂B/∂P| T ≈ ∂B S /∂P| T + 2T βγ, where γ is the Grüneisen constant and B S is the adiabatic bulk modulus 30 . By disregarding small temperature dependencies of ∂B S /∂P and γ, the aforementioned unknown quantity is finally approximated with the expression 2γ(β + T(∂β/∂T| P ) which contains easily accessible terms 30 .
The application of the cBΩ model in binary alloys will depend upon the availability of parameters that can be calculated either by experiment and/or DFT calculations. As it is discussed above, approximations can also be used to overcome the lack of certain parameters thus the calculation of the point defect parameters becomes feasible.
Summary. In the present study, the self-diffusion of Si and Ge in Si 1−x Ge x alloys has been investigated in the framework of the cBΩ thermodynamic model, which allows the calculation of point defect thermodynamic parameters from the bulk properties of the alloy. The consideration of a wide temperature range and the whole composition range of Si 1−x Ge x in conjunction with the excellent agreement of the calculated values as compared to the available experimental data demonstrates the efficacy of the approach.
Here we demonstrate how the cBΩ model can become applicable and provide valuable information for the self-or hetero-diffusion and point defect thermodynamic parameters in binary alloys. This method in conjunction with experiment and/or advanced modeling techniques can be employed in numerous systems including for example solid solution MAX phases, nuclear materials, ternary semiconductors and disordered ionic conductors [56][57][58] .

Methods
In the context of the cBΩ model [28][29][30] , the activation Gibbs free energy g act due to the formation and migration of a point defect in a solid is related to its elastic and expansion properties via: In equation (7), B is the isothermal bulk modulus, Ω stands for the mean atomic volume, while c act is a dimensionless constant which is independent of temperature and pressure but depends on the diffusion mechanism (i.e., vacancy or interstitial) and the host material. The diffusion coefficients D of a single diffusion mechanism exhibiting an Arrhenius behavior are: where f is the diffusion correlation factor which depends on the diffusion mechanism and the crystal structure, g is a geometrical factor, a o is the lattice parameter, ν is the attempt frequency and k B is Boltzmann's constant. The connection of point defect thermodynamic parameters, such as the activation entropy s act and the activation enthalpy h act , to the elastic and expansion properties of the bulk material, which is in essence the cBΩ model, is expressed through 28 where β is the volume thermal expansion coefficient which depends on temperature and pressure. The above important thermodynamic parameters (i.e., s act and h act ) are determined experimentally and thus equations (9) and (10) can be used to validate the cBΩ model and calculate the activation enthalpy and activation entropy when there is insufficient experimental diffusion data. The activation volume, υ act in terms of the cBΩ model 30 , is expressed as act act T act T For the proper implementation of the cBΩ model according to equations (7)(8)(9)(10)(11), the estimation of the constant c act is necessary. At zero temperature, g act equals to h o act and thus, in principle, the constant c act has the value Ω h B / o act o o , where the subscripts refer to T = 0 K 27 . The most reliable method of the calculation of c act is the mean value method which is applicable when experimental diffusion data are available over a broad temperature or pressure range 22,23,30,32 . By taking the natural logarithm of both sides in equation (8), we obtain: act B 0 2 According to equation (12), a linear dependence of lnD versus the quantity BΩ/k B T indicates the validity of the cBΩ model as concern a single diffusion mechanism, and the constant c act arises directly from the slope of the linear fitting.