- Email: [email protected]
First-principles calculation of sulfur–selenium segregation in ZnSe1−xSx: The role of lattice vibration
Materials Science in Semiconductor Processing 39 (2015) 96–102
Contents lists available at ScienceDirect
Materials Science in Semiconductor Processing journal homepage: www.elsevier.com/locate/mssp
First-principles calculation of sulfur–selenium segregation in ZnSe1 xSx: The role of lattice vibration Xiao-Kang Li, Hong-Tao Xue n, Fu-Ling Tang n,1, Wen-Jiang Lu State Key Laboratory of Advanced Processing and Recycling of Non-ferrous Metals, Department of Materials Science and Engineering, Lanzhou University of Technology, Lanzhou 730050, China
a r t i c l e in f o
PACS: 31.15.A 64.75. g 64.75.Qr 63.70.þh. Keywords: First-principles calculation Miscibility gap Phase separation Lattice vibrations
abstract First-principles phase diagram calculations based on density functional theory within the generalized gradient approximation in combination with Monte Carlo (MC) techniques and cluster expansion were performed for the ZnSe1 xSx alloys. Formation energies were used as a basis for fitting cluster expansion Hamiltonians. All formation energies of ZnSe1 xSx alloys are positive, showing that ZnSe1 xSx alloys is a miscibility gap system and has a tendency to phase separation. For ZnSe1 xSx alloys the phase diagram computed with conventional cluster expansion shows a miscibility gap with consolute temperature Tc ¼ 327 K. The contributions of lattice vibration reduce Tc to 281 K (about 15%). We presented a MC study of the spatial distribution of S and Se in ZnSe0.5S0.5 alloys. We found that, the system becomes more homogeneous including lattice vibration at lower temperature. It is consistence with the calculation of phase diagram of ZnSe1 xSx alloys. & 2015 Elsevier Ltd. All rights reserved.
1. Introduction Up to now, there has been a renewed interest in the study of the wide band gap II–VI compound semiconductors. Especially, zinc sulfide (ZnS) and zinc selenide (ZnSe) have received lots of attentions because of the possibility of having a vanishingly small conduction band offset [1] and their potential applications in photo detectors and devices operating in the visible-light range and an UV photo detectors [2]. ZnS and ZnSe have a band gap of 3.66 eV and 2.67 eV [3] respectively, which means that ZnS and ZnSe are potentially good materials for optoelectronic devices [4]. Alloying of binary semiconductors is an important and also the most effective way to change the optical and electronic properties of the alloys artificially.
n
Corresponding authors. E-mail addresses: [email protected] (H.-T. Xue), [email protected] (F.-L. Tang). 1 Tel.: þ 86 931 2973941. http://dx.doi.org/10.1016/j.mssp.2015.03.024 1369-8001/& 2015 Elsevier Ltd. All rights reserved.
In this view, it is necessary to study ZnSe, ZnS and ZnSe1 xSx ternary alloys. Both ZnSe and ZnS have zinc blend structure (space group F43m) at normal conditions [5]. There are many methods such as metal-organic chemical vapor deposition [6], atomic layer epitaxy [7], molecular beam epitaxy [8] and so on can be employed to form the ZnSe1 xSx solid solution. Many studies which have calculated the lattice constant, band gap, optical properties and electronic properties of ZnSe1 xSx compounds can be found in earlier literatures [9–11]. In terms of thermal properties, Benkabou et al. used molecular simulation coupled with the empirical interatomic potentials to calculate the structural and thermodynamic properties of the ternary Alloy ZnSe 1 xSx [12]. Saitta et al. calculated the structural and electronic properties of ZnxMg1 xSySe1 y by using selfconsistence density-functional theory (DFT) within the local-density approximation (LDA) [13]. In their study the phase diagram of the quaternary alloy was determined, and the critical temperature of ZnSeyS1 y (Tc ¼254 K) was in the range of typical MBE growth temperature.
X.-K. Li et al. / Materials Science in Semiconductor Processing 39 (2015) 96–102
The aim of this paper is to present a comparative study of the spatial distribution of S and Se atoms in ZnSe1 xSx alloys with and without the lattice vibrational effects and then to investigate the thermodynamic stability of ZnSe1 xSx system. First-principles phase diagram computational approach allows a direct assessment of the effects of lattice vibrations on phase stability by comparing the phase diagrams that were calculated with and without the inclusion of vibrational effects. The paper can be summarized as follows. In Section 2, we describe the computational method. The results and discussions are given in Section 3. And the paper is concluded in Section 4. 2. Computational methodology 2.1. Cluster expansion The computation of phase equilibrium and ground-state structure of the ZnSe1 xSx alloy requires the calculation of total energy for all possible configurations of placing S and Se atoms on N sites of the Bravais lattice. The cluster expansion (CE) method can be considered as a compact representation of the configurational total energy of an alloy. The accuracy of alloy can be improved by adding the number of clusters in the expansion and the number of energies used in the fit when necessary. Detailed illustration for this method can be found in many literatures [14–18].The main description of the CE method is given. The cluster expansion (CE) method constructs an Ising-like Hamiltonian for the energies of the different atomic configurations. Generally speaking for a binary alloy, the Isinglike model consists of designating occupation variables σi to each site-i of the parent lattice. Due to the type of atom occupying the site, we can take the value þ1 or 1. A specific arrangement of a spin of the parent lattice can be represented by σ containing the value of the occupation variable for each site in the parent lattice. In our work of the pseudobinary ZnSe–ZnS system, the alloy configuration is described by a set of “spin” occupation variables σi, which take values σi ¼ þ1 when site-i is occupied by Se atom and σi ¼ 1 when site-i is occupied by S atom. The CE parameterizes the configurational energy of every exchangeable anion, as a polynomial in “spin” occupation variables: X Eðσ Þ ¼ mα J α ∏i A α0 σ i ; ð1Þ α
where α is a cluster (a set of lattice sites i). The sum is taken over all clusters α that are not equivalent by a symmetry operation of the space group of the parent lattice, while the average is taken over all clusters α0 that are equivalent to α by symmetry. The coefficient Jα in this polynomial embodies the information regarding the energetics of the alloy and is called the effective cluster interaction (ECI). The multiplicity mα indicates the number of clusters that are equivalent to α by symmetry. It is seen that when all clusters α are considered in the sum, the cluster expansion is able to represent the energy E (σ) of any configuration σ by an appropriate selection of the values of Jα. The unknown parameters of the cluster
97
expansion (the ECIs) are then determined by fitting them to the energies of a relatively small number of configurations obtained by first-principles calculations [19]. The CE method was performed with the Alloy Theoretic Automated Toolkit (ATAT) which can automate most of the tasks associated with the construction of the CE Hamiltonian and the calculation of thermodynamic properties [19,20,21]. ATAT proceeds by gradually increasing the number of clusters included in the cluster expansion until the desired accuracy is achieved. Using this method, we studied the phase diagram of CuIn(Se1 xSx)2 and found that CuIn(Se1 xSx)2 alloy is a miscibility gap system and has a tendency to phase separation [22]. In order to evaluate the predictive power of the cluster expansion, the cross-validation (CV) score is introduced as a criterion. The CV score is defined as: CV ¼
!1=2 N 4 2 1X Ei Ei ; Ni¼1
ð2Þ 4
where Ei is the calculated energy of structure i, while Ei is the predicted value of the energy of structure i obtained from a least-squares fit to the (N 1) other structural energies [19]. The CV score is an improvement relative to the standard mean square error criterion that only minimizes the error for structures included in the fit and is not monotonically decreasing. Along with the increase of the number of cluster to be fitted, the CV score first decrease. This is because an increasing number of degrees of freedom are used to explain the fluctuation in energy. Then, before increasing the CV score go through a minimum, because of the decrease in predictive power caused by the increase in the noise in the fitted ECI [19]. A cluster expansion can be considered satisfactory when the optimal CV score can be found, typically less than 0.025 eV.
2.2. Total energy calculations The electronic-structure total-energy calculations of ordered structures, which are required for the construction of the cluster expansion Hamiltonian, were performed using the Vienna Ab initio Simulation Package (VASP) [23,24]. The interaction between the core and the valence electrons was included by the frozen-core projector augmented wave (PAW) method [25,26]. Valence electron configurations for the pseudopotentials are Zn¼3d10 4s2, S¼3s2 3p4, and Se¼4s2 4p4. Both cell constant and ionic positions were fully relaxed in all supercell computations and a plane-wave energy cuff of 380 eV was used. The generalized gradient approximation (GGA) functional of Perdew–Burke–Ernzerhof (PBE) was introduced to describe the exchange-correlation energy [27]. ATAT is used as a script interface to VASP. This script defines a parameter called KPPRA, which automatically sets up the k-point mesh for similar systems. In our work, KPPRA was set to 2000, which, for the primitive cell of face centered cubic structure, translates to a 5 5 5 k-points grid. The convergence tests showed that the increase of energy cutoff and k-points change the total energy by less than 0.01 meV.
98
X.-K. Li et al. / Materials Science in Semiconductor Processing 39 (2015) 96–102
2.3. Lattice vibrations The contribution of the lattice vibration to the free energies can be described by using the bond-lengthdependent transferable force constants [28]. This calculation method has the advantage of reducing the computational burden of calculating phonon densities of states for each supercell configuration that is used to fit the cluster expansion. That is to say, the introduction of this scheme made computational efficiency higher. As this method has been used successfully in several studies [29–32], we used it to calculate the contribution of lattice vibrations to the free energies of ZnSe1 xSx alloys. In brief, this calculation method proceeds as follows: first, the reaction forces from various imposed atomic displacements away from their equilibrium positions are calculated, which is done only for a few high-symmetry ordered supercells for a range of lattice parameters. Secondly, the parameters defining the bond-length-dependent transferable force constants relation can be obtained from a polynomial fit of the calculated forces as a function of bond length. Finally, once the bond-length dependence of bond stiffness is known, one can predict the nearest-neighbor interatomic force constants. The configuration dependence of the vibrational free energy is parameterized with a cluster expansion with temperature-dependent effective cluster interactions [30].
Fig. 1. Cluster expansion and direct enumeration ground-state search results of formation energies of ZnSe1 xSx. “Known str” denotes structures which energies have been calculated from first principle. “Predicted” denotes structures which energy has been predicted from the cluster expansion for which the ab initio energy is known.
2.4. Monte Carlo simulations In order to obtain the phase diagram of the ZnSe1 xSx alloys, we use the conventional Monte Carlo (MC) simulations. The energy of the system is specified by the cluster expansion Hamiltonian. The MC simulations were performed using the phb program of ATAT. This program can compute the two-phase equilibrium. When the “end” of the two-phase equilibrium has been reached, the phb can terminate automatically. We used a 25 25 25 supercell for the direct MC simulations. And the averaging times for the given precision on the average concentration of the alloy was set to 0.1% in the current simulations. Detailed description of the algorithm underlying this program can be found in Ref. [20]. In addition, we also study the spatial distribution of S and Se atoms in ZnSe1 xSx alloys by MC simulation, which was performed by using the emc2 program of ATAT. To make the simulation feasible, we keep S and Se atoms fixed at their respective lattice sites and neglect other defects. The MC moves consist of exchanging the position of the S and Se atoms. The model box contains 18 18 18 primitive cells (5832 active atoms). The simulation temperature ranges from 100 K to 500 K. 3. Results and discussion 3.1. Ground-state structures and effective cluster interactions In Fig. 1, we obtained the static formation energies ΔEf (σ) (0 K) by first-principle computations, and they will be used to fit the cluster expansion Hamiltonians. The
Å
Fig. 2. Effective cluster interactions (ECIs) of ZnSe1 xSx alloys.
formation energy
ΔEf (σ) was calculated as follows:
ΔEf ðσ Þ ¼ Eðσ Þ ½ð1 xÞEðZnSeÞ þ xEðZnSÞ:
ð3Þ
where E(σ) is the total energy (per anion) of a given configuration of ZnSe1 xSx, E(ZnSe) is the energy (per anion) of ZnSe, and E(ZnS) is the energy (per anion) of ZnS. From Fig. 1, we can also see that all formation energies of all ordered structures are positive, which indicates a miscibility gap system. The cluster expansions were fit to moderate number of input structures to obtain a good cross validation score. The small CV value indicates that the predictive power of the cluster expansion is excellent in this work. The effective cluster interactions (ECIs) of ZnSe1 xSx alloys obtained from cluster expansion are shown in Fig. 2. The cluster expansion of ZnSe1 xSx includes interactions up to sixth neighbor pairs and three three-body terms to obtain the desired CV score. In Fig. 2, the horizontal axis represents the cluster diameter. The first is for pair clusters beginning with the label “pair”, and the next is for triplet clusters. Cluster diameter was defined as the length of the longest pair
X.-K. Li et al. / Materials Science in Semiconductor Processing 39 (2015) 96–102
contained in the cluster. It is evident that the magnitude of the ECI generally decreases as the distance and the cluster size increases, and the convergence is obtained. The presence of three-body terms (Table 1) in the cluster expansion Hamiltonians for implies asymmetric immiscibility for ZnSe1 xSx. The dependence of the stretch and bending force constant on bond length are shown in Fig. 3. A linear relationship was found to provide a reliable description of the nearest-neighbor force constants in ZnSe1 xSx system. We can see that the stretching constants decrease monotonically with increasing bond length, by contrast the bending terms are not sensitive. Once the bond-length
Table 1 Characteristics of the calculated cluster expansions. ZnSe1 xSx
Number of structures Number of clusters Cross-validation scores (meV/atom)
33 2þ6 þ3 0.174
Å
Characteristics
dependence of bond stiffness is known, one can predict the nearest-neighbor interatomic force constants for any supercell configuration in the cluster expansion fit from the relaxed bond lengths which are obtained from the VASP structure energy minimizations. The standard lattice dynamics based on a nearest-neighbor bond Born–von Kármán model then provides the phonon density of states and, consequently, any thermodynamic property of interest, such as the contribution of lattice vibrations to the free energy.
3.2. Phase diagram Subsolidus phase diagrams of ZnSe1 xSx alloys computed with and without lattice vibration effects are shown in Fig. 4, where the miscibility gaps (binodal curves) correspond to the coexistence of two stable phases. Inside the miscibility gaps the phase separation occurs, and the mixture of two phases of ZnSe and ZnS is stable. Outside the miscibility gaps, single-phase homogenous ZnSe1 xSx is stable. As shown in Fig. 4, the consolute temperature Tc predicted without vibrational contribution is 327 K (xc ¼ 0.571). The ZnSe1 xSx system exhibits slight asymmetric, with maximum closer to the end member ZnS. And our result also shows that the ZnSe1 xSx alloys can exhibit the phase separation at low temperature, the tendency to phase separation can cause alloy inhomogeneity to a certain extent. From Fig. 4, we also can see that the contribution of lattice vibrations reduces the consolute temperature to 281 K (xc E(co). And the miscibility gap is also asymmetric. The xc with the vibrational effect (0.543) is smaller than that without the vibrational effect (0.571). The miscibility gap with the effect of vibration exhibits subtle increase symmetry than that without the effect of vibration. And the consolute temperature of 281 K is only 27 K higher than the result of 254 K obtained by Saitta [13]. Such results about the critical temperatures are in agreement with the few available experimental findings on the coexistence curves of such alloys and indicate the complete miscibility of the ZnSe1 xSx system. By contrast, Tc is reduced by about 15%, by considering the vibrational effect. The reduction in Tc, ΔTc, was
Å
Å
Å
Fig. 3. Nearest-neighbor bond stiffness (against stretching and bending) as a function of bond length for ZnSe1 xSx alloys: (a) Se–Zn bond, (b) S–Zn bond. Points indicate ab initio data and lines are linear fits used in the calculations of the vibrational free energy.
99
Fig. 4. The calculated phase diagram of ZnSe1 xSx alloys.
100
X.-K. Li et al. / Materials Science in Semiconductor Processing 39 (2015) 96–102
calculated according to the formula:
ΔT C ¼
2 T C T VC ib T C þ T VC ib
;
ð4Þ
where T Vib and Tc are the values for Tc that were calculated c with and without lattice vibrations, respectively. The result (ΔTc ¼15%) obtained by us for the ZnSe1 xSx alloys is consistent with the result that including the effect of
Fig. 5. Snapshot of the ZnSe0.5S0.5 with (on the right) and without (on the left) lattice vibration effects at the temperature of 250 K, 275 K, 300 K, and 350 K. Yellow and blue spheres represent S and Se atoms, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
X.-K. Li et al. / Materials Science in Semiconductor Processing 39 (2015) 96–102
lattice vibrations typically leads to a 5–15% reduction in the calculated value for Tc [33]. We think the atomic size mismatch should be the most probable cause to explain the vibrational effect and the asymmetry in the ZnSe1 xSx system. The smaller S atoms (atomic radius 1.04 Å) [34] are under tension and the larger Se atoms (atomic radius 1.17 Å) [34] are under compression. When larger Se atoms sit on neighboring crystal lattices, extra space will be occupied and the amplitude of lattice vibrations is reduced. The ZnSe1 xSx alloy systems tend to be stiffer. On the contrary, when smaller S atoms sit on neighboring crystal lattices, there will be the extra available room which results in a locally softer region in the ZnSe1 xSx alloy system and thus have large vibrational entropy. Therefore, it should be reasonable to use the size mismatch between S and Se atoms to explain the large vibrational effect in ZnSe1 xSx alloys. 3.3. Distribution of S and Se atoms in ZnSe1 xSx alloys We have calculated the spatial distribution of S and Se atoms in ZnSe0.5S0.5 alloys at temperatures of 100 K, 150 K, 200 K, 250 K, 275 K, 300 K, 350 K, 400 K, 450 K, and 500 K by MC simulation, respectively. As the results are similar with each other when the temperature below 250 K and over 350 K, in order to save space, we just displayed several representative results at temperatures of 250 K, 275 K, 300 K, and 350 K, as shown in Fig. 5. As shown in the left part of the Fig. 5, at low temperature, for example, at 250 K obvious phase separation between ZnS and ZnSe can be observed in ZnSe0.5S0.5 alloy, which means that there are two phases in ZnSe0.5S0.5 alloy: a rich-S and a poor-Se ZnSe1 xSx phase. As temperature increase, ZnSe0.5S0.5 becomes more and more homogeneous especially over 300 K. Similar observation can also be found in the right part of the Fig. 5 which shows the spatial distribution of S and Se atoms in ZnSe0.5S0.5 alloy with lattice vibration effects. In order to further analyze the inhomogeneity of ZnSe0.5S0.5 alloy at various temperatures quantitatively, we divided the simulation box into 729 cubic segments and each segment has 8 active atoms (S and Se atoms). The number of S (Se) atoms in each segment was counted and the standard deviation σ was used to describe the inhomogeneity. The dependence of the standard deviation on temperature is shown in Fig. 6. We can see that the standard deviation typically decreases with increasing temperature for both cases that consider the lattice vibration effects or not. It means that the homogeneity of ZnSe0.5S0.5 alloy increases with the increase of the temperature. In Fig. 6, we also can see that σ with lattice vibration effects is lower than that without lattice vibration effects at temperature below 300 K, which indicates that the alloy considering the lattice vibration effects is more homogeneous than that without considering the lattice vibration effects at low temperature. Hence, the effect of the lattice vibration on the homogeneity of ZnSe1 xSx alloys should not be neglected, especially at low temperature. It may also have some influences on the temperature of phase separation or homogeneity for the ZnSe1 xSx alloys. In addition, the inflection points on the
101
Fig. 6. Standard deviation σ of ZnSe0.5S0.5 as a function of temperature with (red) and without (black) lattice vibration effects. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
curves in Fig. 6 likely indicate that the system will undergo a phase transition near them. 4. Conclusions We have studied the phase equilibrium of ZnSe1 xSx alloys system based on density functional theory within the generalized gradient approximation combination Monte Carlo and cluster expansion techniques. For this pseudobinary alloy system, all formation energies are positive, indicating a phase separation tendency. For ZnSe1 xSx alloys the phase diagram computed with conventional cluster expansion shows a miscibility gap with consolute temperature Tc ¼327 K. The contributions of the lattice vibration reduce Tc to 281 K (about 15%). The miscibility gaps for ZnSe1 xSx system are predicted to be asymmetric with maximum closer to the end member with the smaller ion (S). The miscibility gap with the effect of vibration exhibits subtle increase symmetry than that without the effect of vibration. We studied the spatial distribution of S and Se atoms in ZnSe0.5S0.5 alloy, and found that, the system becomes more homogeneous including lattice vibration at low temperature.
Acknowledgments This work was financially supported by the National Natural Science Foundation of China (11164014 and 11364025) and Gansu Science and Technology Pillar Program (1204GKCA057). This work was performed in Gansu Supercomputer Centre. References [1] A. Qteish, R. Said, N. Meskini, A. Nazzal, Phys. Rev. B 52 (1995) 1830–1838. [2] J. Jeon, J. Ding, A.V. Nurmikko, W. Xie, D.C. Griio, M. Kobayashi, L. Gunshor, G.C. Hau, N. Otsuka, Appl. Phys. Lett. 60 (1992) 892–894. [3] C. Persson, A. Zunger, Phys. Rev. B 68 (2003) 073205. -073205.
102
X.-K. Li et al. / Materials Science in Semiconductor Processing 39 (2015) 96–102
[4] J.A. Zapien, Y. Jiang, X.M. Meng, W. Chen, F.C.K. Au, Y. Lifshitz S.T. Lee, Appl. Phys. Lett. 84 (2004) 1189–1191. [5] J.L. Bredin, Phys. Today 47 (1994) 5. [6] P.J. Wright, B. Cockayne, A.F. Cattle, P.J. Dean, A.D. Pitt, J. Cryst. Growth 79 (1986) 357–362. [7] C.T. Hsu, Mater. Chem. Phys. 58 (1999) 6–11. [8] N. Matsumura, M. Tsubokura, J. Saraie, Y. Yo-dogawa, J. Cryst. Growth 86 (1988) 311–317. [9] E. Engel, S.H. Vosko, Phys. Rev. B 50 (1994) 10498–10505. [10] F. El Haj Hassana, B. Amranib, F. Bahsoun, Physica B 391 (2007) 363–370. [11] Z. Nourbakhsh, Physica B 405 (2010) 4173–4187. [12] F.Z. Benkabou, Appl. Phys. Res. 3 (2011) 211–223. [13] A.M. Saitta, S. de G ironcoli, S. Baroni, Phys. Rev. Lett. 80 (1998) 4939–4942. [14] J.M. Sanchez, F. Ducastelle, D. Gratias, Physica A 128 (1984) 334–350. [15] Z.W. Lu, S.H. Wei, A. Zunger, S. Frota-Pessoa, L.G. Ferreira, Phys. Rev. B 44 (1991) 512–544. [16] A. Zunger, L.G. Wang, G.L.W. Hart, M. Sanati, Model. Simul. Mater. Sci. Eng. 10 (2002) 685–706. [17] J.M. Sanchez, Phys. Rev. B 81 (2010) 224202-1–224202-13. [18] C. Ravi, B.K. Panigrahi, M.C. Valsakumar, A. van de Walle, Phys. Rev. B 81 (2010) 104111-1–104111-13.
[19] A. van de Walle, M. Asta, G. Ceder, Calphad 26 (2002) 539–553. [20] A. van de Walle, M. Asta, Model. Simul. Mater. Sci. Eng. 10 (2002) 521–538. [21] A. van de Walle, G. Ceder, J. Phase Equilibria 23 (2002) 348–359. [22] H.T. Xue, F.L. Tang, X.K. Li, F.C. Wan, W.J. Lu, Z.Y. Rui, Y.D. Feng, Mater. Sci. Semicond. Process. 25 (2014) 251–257. [23] G. Kresse, J. Hafner, Phys. Rev. B 47 (1993) 558–561. [24] G. Kresse, J. Furthmuller, Phys. Rev. B 54 (1996) 11169–11186. [25] P.E. Blöchl, Phys. Rev. B 50 (1994) 17953–17979. [26] G. Kresse, D. Joubert, Phys. Rev. B 59 (1999) 1758–1775. [27] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865–3868. [28] A. van de Walle, G. Ceder, Rev. Mod. Phys. 74 (2002) 11–45. [29] A. van de Walle, G. Ceder, Phys. Rev. B 61 (2000) 5972–5978. [30] B.P. Burton, A. van de Walle, Chem. Geol. 225 (2006) 222–229. [31] O. Adjaoud, G.S. Neumann, B.P. Burton, A. van de Walle, Phys. Rev. B 80 (2009) 134112-1–134112-7. [32] P. Dalach, D.E. Ellis, A. van de Walle, Phys. Rev. B 82 (2010) 144117-1–144117-11. [33] B.P. Burton, A. van de Walle, U. Kattner, J. Appl. Phys. 100 (2006) 113528-1–113528-6. [34] M.K. Agarwal, L.T. Talele, Solid State Commun. 59 (1986) 549–551.
Contents lists available at ScienceDirect
Materials Science in Semiconductor Processing journal homepage: www.elsevier.com/locate/mssp
First-principles calculation of sulfur–selenium segregation in ZnSe1 xSx: The role of lattice vibration Xiao-Kang Li, Hong-Tao Xue n, Fu-Ling Tang n,1, Wen-Jiang Lu State Key Laboratory of Advanced Processing and Recycling of Non-ferrous Metals, Department of Materials Science and Engineering, Lanzhou University of Technology, Lanzhou 730050, China
a r t i c l e in f o
PACS: 31.15.A 64.75. g 64.75.Qr 63.70.þh. Keywords: First-principles calculation Miscibility gap Phase separation Lattice vibrations
abstract First-principles phase diagram calculations based on density functional theory within the generalized gradient approximation in combination with Monte Carlo (MC) techniques and cluster expansion were performed for the ZnSe1 xSx alloys. Formation energies were used as a basis for fitting cluster expansion Hamiltonians. All formation energies of ZnSe1 xSx alloys are positive, showing that ZnSe1 xSx alloys is a miscibility gap system and has a tendency to phase separation. For ZnSe1 xSx alloys the phase diagram computed with conventional cluster expansion shows a miscibility gap with consolute temperature Tc ¼ 327 K. The contributions of lattice vibration reduce Tc to 281 K (about 15%). We presented a MC study of the spatial distribution of S and Se in ZnSe0.5S0.5 alloys. We found that, the system becomes more homogeneous including lattice vibration at lower temperature. It is consistence with the calculation of phase diagram of ZnSe1 xSx alloys. & 2015 Elsevier Ltd. All rights reserved.
1. Introduction Up to now, there has been a renewed interest in the study of the wide band gap II–VI compound semiconductors. Especially, zinc sulfide (ZnS) and zinc selenide (ZnSe) have received lots of attentions because of the possibility of having a vanishingly small conduction band offset [1] and their potential applications in photo detectors and devices operating in the visible-light range and an UV photo detectors [2]. ZnS and ZnSe have a band gap of 3.66 eV and 2.67 eV [3] respectively, which means that ZnS and ZnSe are potentially good materials for optoelectronic devices [4]. Alloying of binary semiconductors is an important and also the most effective way to change the optical and electronic properties of the alloys artificially.
n
Corresponding authors. E-mail addresses: [email protected] (H.-T. Xue), [email protected] (F.-L. Tang). 1 Tel.: þ 86 931 2973941. http://dx.doi.org/10.1016/j.mssp.2015.03.024 1369-8001/& 2015 Elsevier Ltd. All rights reserved.
In this view, it is necessary to study ZnSe, ZnS and ZnSe1 xSx ternary alloys. Both ZnSe and ZnS have zinc blend structure (space group F43m) at normal conditions [5]. There are many methods such as metal-organic chemical vapor deposition [6], atomic layer epitaxy [7], molecular beam epitaxy [8] and so on can be employed to form the ZnSe1 xSx solid solution. Many studies which have calculated the lattice constant, band gap, optical properties and electronic properties of ZnSe1 xSx compounds can be found in earlier literatures [9–11]. In terms of thermal properties, Benkabou et al. used molecular simulation coupled with the empirical interatomic potentials to calculate the structural and thermodynamic properties of the ternary Alloy ZnSe 1 xSx [12]. Saitta et al. calculated the structural and electronic properties of ZnxMg1 xSySe1 y by using selfconsistence density-functional theory (DFT) within the local-density approximation (LDA) [13]. In their study the phase diagram of the quaternary alloy was determined, and the critical temperature of ZnSeyS1 y (Tc ¼254 K) was in the range of typical MBE growth temperature.
X.-K. Li et al. / Materials Science in Semiconductor Processing 39 (2015) 96–102
The aim of this paper is to present a comparative study of the spatial distribution of S and Se atoms in ZnSe1 xSx alloys with and without the lattice vibrational effects and then to investigate the thermodynamic stability of ZnSe1 xSx system. First-principles phase diagram computational approach allows a direct assessment of the effects of lattice vibrations on phase stability by comparing the phase diagrams that were calculated with and without the inclusion of vibrational effects. The paper can be summarized as follows. In Section 2, we describe the computational method. The results and discussions are given in Section 3. And the paper is concluded in Section 4. 2. Computational methodology 2.1. Cluster expansion The computation of phase equilibrium and ground-state structure of the ZnSe1 xSx alloy requires the calculation of total energy for all possible configurations of placing S and Se atoms on N sites of the Bravais lattice. The cluster expansion (CE) method can be considered as a compact representation of the configurational total energy of an alloy. The accuracy of alloy can be improved by adding the number of clusters in the expansion and the number of energies used in the fit when necessary. Detailed illustration for this method can be found in many literatures [14–18].The main description of the CE method is given. The cluster expansion (CE) method constructs an Ising-like Hamiltonian for the energies of the different atomic configurations. Generally speaking for a binary alloy, the Isinglike model consists of designating occupation variables σi to each site-i of the parent lattice. Due to the type of atom occupying the site, we can take the value þ1 or 1. A specific arrangement of a spin of the parent lattice can be represented by σ containing the value of the occupation variable for each site in the parent lattice. In our work of the pseudobinary ZnSe–ZnS system, the alloy configuration is described by a set of “spin” occupation variables σi, which take values σi ¼ þ1 when site-i is occupied by Se atom and σi ¼ 1 when site-i is occupied by S atom. The CE parameterizes the configurational energy of every exchangeable anion, as a polynomial in “spin” occupation variables: X Eðσ Þ ¼ mα J α ∏i A α0 σ i ; ð1Þ α
where α is a cluster (a set of lattice sites i). The sum is taken over all clusters α that are not equivalent by a symmetry operation of the space group of the parent lattice, while the average is taken over all clusters α0 that are equivalent to α by symmetry. The coefficient Jα in this polynomial embodies the information regarding the energetics of the alloy and is called the effective cluster interaction (ECI). The multiplicity mα indicates the number of clusters that are equivalent to α by symmetry. It is seen that when all clusters α are considered in the sum, the cluster expansion is able to represent the energy E (σ) of any configuration σ by an appropriate selection of the values of Jα. The unknown parameters of the cluster
97
expansion (the ECIs) are then determined by fitting them to the energies of a relatively small number of configurations obtained by first-principles calculations [19]. The CE method was performed with the Alloy Theoretic Automated Toolkit (ATAT) which can automate most of the tasks associated with the construction of the CE Hamiltonian and the calculation of thermodynamic properties [19,20,21]. ATAT proceeds by gradually increasing the number of clusters included in the cluster expansion until the desired accuracy is achieved. Using this method, we studied the phase diagram of CuIn(Se1 xSx)2 and found that CuIn(Se1 xSx)2 alloy is a miscibility gap system and has a tendency to phase separation [22]. In order to evaluate the predictive power of the cluster expansion, the cross-validation (CV) score is introduced as a criterion. The CV score is defined as: CV ¼
!1=2 N 4 2 1X Ei Ei ; Ni¼1
ð2Þ 4
where Ei is the calculated energy of structure i, while Ei is the predicted value of the energy of structure i obtained from a least-squares fit to the (N 1) other structural energies [19]. The CV score is an improvement relative to the standard mean square error criterion that only minimizes the error for structures included in the fit and is not monotonically decreasing. Along with the increase of the number of cluster to be fitted, the CV score first decrease. This is because an increasing number of degrees of freedom are used to explain the fluctuation in energy. Then, before increasing the CV score go through a minimum, because of the decrease in predictive power caused by the increase in the noise in the fitted ECI [19]. A cluster expansion can be considered satisfactory when the optimal CV score can be found, typically less than 0.025 eV.
2.2. Total energy calculations The electronic-structure total-energy calculations of ordered structures, which are required for the construction of the cluster expansion Hamiltonian, were performed using the Vienna Ab initio Simulation Package (VASP) [23,24]. The interaction between the core and the valence electrons was included by the frozen-core projector augmented wave (PAW) method [25,26]. Valence electron configurations for the pseudopotentials are Zn¼3d10 4s2, S¼3s2 3p4, and Se¼4s2 4p4. Both cell constant and ionic positions were fully relaxed in all supercell computations and a plane-wave energy cuff of 380 eV was used. The generalized gradient approximation (GGA) functional of Perdew–Burke–Ernzerhof (PBE) was introduced to describe the exchange-correlation energy [27]. ATAT is used as a script interface to VASP. This script defines a parameter called KPPRA, which automatically sets up the k-point mesh for similar systems. In our work, KPPRA was set to 2000, which, for the primitive cell of face centered cubic structure, translates to a 5 5 5 k-points grid. The convergence tests showed that the increase of energy cutoff and k-points change the total energy by less than 0.01 meV.
98
X.-K. Li et al. / Materials Science in Semiconductor Processing 39 (2015) 96–102
2.3. Lattice vibrations The contribution of the lattice vibration to the free energies can be described by using the bond-lengthdependent transferable force constants [28]. This calculation method has the advantage of reducing the computational burden of calculating phonon densities of states for each supercell configuration that is used to fit the cluster expansion. That is to say, the introduction of this scheme made computational efficiency higher. As this method has been used successfully in several studies [29–32], we used it to calculate the contribution of lattice vibrations to the free energies of ZnSe1 xSx alloys. In brief, this calculation method proceeds as follows: first, the reaction forces from various imposed atomic displacements away from their equilibrium positions are calculated, which is done only for a few high-symmetry ordered supercells for a range of lattice parameters. Secondly, the parameters defining the bond-length-dependent transferable force constants relation can be obtained from a polynomial fit of the calculated forces as a function of bond length. Finally, once the bond-length dependence of bond stiffness is known, one can predict the nearest-neighbor interatomic force constants. The configuration dependence of the vibrational free energy is parameterized with a cluster expansion with temperature-dependent effective cluster interactions [30].
Fig. 1. Cluster expansion and direct enumeration ground-state search results of formation energies of ZnSe1 xSx. “Known str” denotes structures which energies have been calculated from first principle. “Predicted” denotes structures which energy has been predicted from the cluster expansion for which the ab initio energy is known.
2.4. Monte Carlo simulations In order to obtain the phase diagram of the ZnSe1 xSx alloys, we use the conventional Monte Carlo (MC) simulations. The energy of the system is specified by the cluster expansion Hamiltonian. The MC simulations were performed using the phb program of ATAT. This program can compute the two-phase equilibrium. When the “end” of the two-phase equilibrium has been reached, the phb can terminate automatically. We used a 25 25 25 supercell for the direct MC simulations. And the averaging times for the given precision on the average concentration of the alloy was set to 0.1% in the current simulations. Detailed description of the algorithm underlying this program can be found in Ref. [20]. In addition, we also study the spatial distribution of S and Se atoms in ZnSe1 xSx alloys by MC simulation, which was performed by using the emc2 program of ATAT. To make the simulation feasible, we keep S and Se atoms fixed at their respective lattice sites and neglect other defects. The MC moves consist of exchanging the position of the S and Se atoms. The model box contains 18 18 18 primitive cells (5832 active atoms). The simulation temperature ranges from 100 K to 500 K. 3. Results and discussion 3.1. Ground-state structures and effective cluster interactions In Fig. 1, we obtained the static formation energies ΔEf (σ) (0 K) by first-principle computations, and they will be used to fit the cluster expansion Hamiltonians. The
Å
Fig. 2. Effective cluster interactions (ECIs) of ZnSe1 xSx alloys.
formation energy
ΔEf (σ) was calculated as follows:
ΔEf ðσ Þ ¼ Eðσ Þ ½ð1 xÞEðZnSeÞ þ xEðZnSÞ:
ð3Þ
where E(σ) is the total energy (per anion) of a given configuration of ZnSe1 xSx, E(ZnSe) is the energy (per anion) of ZnSe, and E(ZnS) is the energy (per anion) of ZnS. From Fig. 1, we can also see that all formation energies of all ordered structures are positive, which indicates a miscibility gap system. The cluster expansions were fit to moderate number of input structures to obtain a good cross validation score. The small CV value indicates that the predictive power of the cluster expansion is excellent in this work. The effective cluster interactions (ECIs) of ZnSe1 xSx alloys obtained from cluster expansion are shown in Fig. 2. The cluster expansion of ZnSe1 xSx includes interactions up to sixth neighbor pairs and three three-body terms to obtain the desired CV score. In Fig. 2, the horizontal axis represents the cluster diameter. The first is for pair clusters beginning with the label “pair”, and the next is for triplet clusters. Cluster diameter was defined as the length of the longest pair
X.-K. Li et al. / Materials Science in Semiconductor Processing 39 (2015) 96–102
contained in the cluster. It is evident that the magnitude of the ECI generally decreases as the distance and the cluster size increases, and the convergence is obtained. The presence of three-body terms (Table 1) in the cluster expansion Hamiltonians for implies asymmetric immiscibility for ZnSe1 xSx. The dependence of the stretch and bending force constant on bond length are shown in Fig. 3. A linear relationship was found to provide a reliable description of the nearest-neighbor force constants in ZnSe1 xSx system. We can see that the stretching constants decrease monotonically with increasing bond length, by contrast the bending terms are not sensitive. Once the bond-length
Table 1 Characteristics of the calculated cluster expansions. ZnSe1 xSx
Number of structures Number of clusters Cross-validation scores (meV/atom)
33 2þ6 þ3 0.174
Å
Characteristics
dependence of bond stiffness is known, one can predict the nearest-neighbor interatomic force constants for any supercell configuration in the cluster expansion fit from the relaxed bond lengths which are obtained from the VASP structure energy minimizations. The standard lattice dynamics based on a nearest-neighbor bond Born–von Kármán model then provides the phonon density of states and, consequently, any thermodynamic property of interest, such as the contribution of lattice vibrations to the free energy.
3.2. Phase diagram Subsolidus phase diagrams of ZnSe1 xSx alloys computed with and without lattice vibration effects are shown in Fig. 4, where the miscibility gaps (binodal curves) correspond to the coexistence of two stable phases. Inside the miscibility gaps the phase separation occurs, and the mixture of two phases of ZnSe and ZnS is stable. Outside the miscibility gaps, single-phase homogenous ZnSe1 xSx is stable. As shown in Fig. 4, the consolute temperature Tc predicted without vibrational contribution is 327 K (xc ¼ 0.571). The ZnSe1 xSx system exhibits slight asymmetric, with maximum closer to the end member ZnS. And our result also shows that the ZnSe1 xSx alloys can exhibit the phase separation at low temperature, the tendency to phase separation can cause alloy inhomogeneity to a certain extent. From Fig. 4, we also can see that the contribution of lattice vibrations reduces the consolute temperature to 281 K (xc E(co). And the miscibility gap is also asymmetric. The xc with the vibrational effect (0.543) is smaller than that without the vibrational effect (0.571). The miscibility gap with the effect of vibration exhibits subtle increase symmetry than that without the effect of vibration. And the consolute temperature of 281 K is only 27 K higher than the result of 254 K obtained by Saitta [13]. Such results about the critical temperatures are in agreement with the few available experimental findings on the coexistence curves of such alloys and indicate the complete miscibility of the ZnSe1 xSx system. By contrast, Tc is reduced by about 15%, by considering the vibrational effect. The reduction in Tc, ΔTc, was
Å
Å
Å
Fig. 3. Nearest-neighbor bond stiffness (against stretching and bending) as a function of bond length for ZnSe1 xSx alloys: (a) Se–Zn bond, (b) S–Zn bond. Points indicate ab initio data and lines are linear fits used in the calculations of the vibrational free energy.
99
Fig. 4. The calculated phase diagram of ZnSe1 xSx alloys.
100
X.-K. Li et al. / Materials Science in Semiconductor Processing 39 (2015) 96–102
calculated according to the formula:
ΔT C ¼
2 T C T VC ib T C þ T VC ib
;
ð4Þ
where T Vib and Tc are the values for Tc that were calculated c with and without lattice vibrations, respectively. The result (ΔTc ¼15%) obtained by us for the ZnSe1 xSx alloys is consistent with the result that including the effect of
Fig. 5. Snapshot of the ZnSe0.5S0.5 with (on the right) and without (on the left) lattice vibration effects at the temperature of 250 K, 275 K, 300 K, and 350 K. Yellow and blue spheres represent S and Se atoms, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
X.-K. Li et al. / Materials Science in Semiconductor Processing 39 (2015) 96–102
lattice vibrations typically leads to a 5–15% reduction in the calculated value for Tc [33]. We think the atomic size mismatch should be the most probable cause to explain the vibrational effect and the asymmetry in the ZnSe1 xSx system. The smaller S atoms (atomic radius 1.04 Å) [34] are under tension and the larger Se atoms (atomic radius 1.17 Å) [34] are under compression. When larger Se atoms sit on neighboring crystal lattices, extra space will be occupied and the amplitude of lattice vibrations is reduced. The ZnSe1 xSx alloy systems tend to be stiffer. On the contrary, when smaller S atoms sit on neighboring crystal lattices, there will be the extra available room which results in a locally softer region in the ZnSe1 xSx alloy system and thus have large vibrational entropy. Therefore, it should be reasonable to use the size mismatch between S and Se atoms to explain the large vibrational effect in ZnSe1 xSx alloys. 3.3. Distribution of S and Se atoms in ZnSe1 xSx alloys We have calculated the spatial distribution of S and Se atoms in ZnSe0.5S0.5 alloys at temperatures of 100 K, 150 K, 200 K, 250 K, 275 K, 300 K, 350 K, 400 K, 450 K, and 500 K by MC simulation, respectively. As the results are similar with each other when the temperature below 250 K and over 350 K, in order to save space, we just displayed several representative results at temperatures of 250 K, 275 K, 300 K, and 350 K, as shown in Fig. 5. As shown in the left part of the Fig. 5, at low temperature, for example, at 250 K obvious phase separation between ZnS and ZnSe can be observed in ZnSe0.5S0.5 alloy, which means that there are two phases in ZnSe0.5S0.5 alloy: a rich-S and a poor-Se ZnSe1 xSx phase. As temperature increase, ZnSe0.5S0.5 becomes more and more homogeneous especially over 300 K. Similar observation can also be found in the right part of the Fig. 5 which shows the spatial distribution of S and Se atoms in ZnSe0.5S0.5 alloy with lattice vibration effects. In order to further analyze the inhomogeneity of ZnSe0.5S0.5 alloy at various temperatures quantitatively, we divided the simulation box into 729 cubic segments and each segment has 8 active atoms (S and Se atoms). The number of S (Se) atoms in each segment was counted and the standard deviation σ was used to describe the inhomogeneity. The dependence of the standard deviation on temperature is shown in Fig. 6. We can see that the standard deviation typically decreases with increasing temperature for both cases that consider the lattice vibration effects or not. It means that the homogeneity of ZnSe0.5S0.5 alloy increases with the increase of the temperature. In Fig. 6, we also can see that σ with lattice vibration effects is lower than that without lattice vibration effects at temperature below 300 K, which indicates that the alloy considering the lattice vibration effects is more homogeneous than that without considering the lattice vibration effects at low temperature. Hence, the effect of the lattice vibration on the homogeneity of ZnSe1 xSx alloys should not be neglected, especially at low temperature. It may also have some influences on the temperature of phase separation or homogeneity for the ZnSe1 xSx alloys. In addition, the inflection points on the
101
Fig. 6. Standard deviation σ of ZnSe0.5S0.5 as a function of temperature with (red) and without (black) lattice vibration effects. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
curves in Fig. 6 likely indicate that the system will undergo a phase transition near them. 4. Conclusions We have studied the phase equilibrium of ZnSe1 xSx alloys system based on density functional theory within the generalized gradient approximation combination Monte Carlo and cluster expansion techniques. For this pseudobinary alloy system, all formation energies are positive, indicating a phase separation tendency. For ZnSe1 xSx alloys the phase diagram computed with conventional cluster expansion shows a miscibility gap with consolute temperature Tc ¼327 K. The contributions of the lattice vibration reduce Tc to 281 K (about 15%). The miscibility gaps for ZnSe1 xSx system are predicted to be asymmetric with maximum closer to the end member with the smaller ion (S). The miscibility gap with the effect of vibration exhibits subtle increase symmetry than that without the effect of vibration. We studied the spatial distribution of S and Se atoms in ZnSe0.5S0.5 alloy, and found that, the system becomes more homogeneous including lattice vibration at low temperature.
Acknowledgments This work was financially supported by the National Natural Science Foundation of China (11164014 and 11364025) and Gansu Science and Technology Pillar Program (1204GKCA057). This work was performed in Gansu Supercomputer Centre. References [1] A. Qteish, R. Said, N. Meskini, A. Nazzal, Phys. Rev. B 52 (1995) 1830–1838. [2] J. Jeon, J. Ding, A.V. Nurmikko, W. Xie, D.C. Griio, M. Kobayashi, L. Gunshor, G.C. Hau, N. Otsuka, Appl. Phys. Lett. 60 (1992) 892–894. [3] C. Persson, A. Zunger, Phys. Rev. B 68 (2003) 073205. -073205.
102
X.-K. Li et al. / Materials Science in Semiconductor Processing 39 (2015) 96–102
[4] J.A. Zapien, Y. Jiang, X.M. Meng, W. Chen, F.C.K. Au, Y. Lifshitz S.T. Lee, Appl. Phys. Lett. 84 (2004) 1189–1191. [5] J.L. Bredin, Phys. Today 47 (1994) 5. [6] P.J. Wright, B. Cockayne, A.F. Cattle, P.J. Dean, A.D. Pitt, J. Cryst. Growth 79 (1986) 357–362. [7] C.T. Hsu, Mater. Chem. Phys. 58 (1999) 6–11. [8] N. Matsumura, M. Tsubokura, J. Saraie, Y. Yo-dogawa, J. Cryst. Growth 86 (1988) 311–317. [9] E. Engel, S.H. Vosko, Phys. Rev. B 50 (1994) 10498–10505. [10] F. El Haj Hassana, B. Amranib, F. Bahsoun, Physica B 391 (2007) 363–370. [11] Z. Nourbakhsh, Physica B 405 (2010) 4173–4187. [12] F.Z. Benkabou, Appl. Phys. Res. 3 (2011) 211–223. [13] A.M. Saitta, S. de G ironcoli, S. Baroni, Phys. Rev. Lett. 80 (1998) 4939–4942. [14] J.M. Sanchez, F. Ducastelle, D. Gratias, Physica A 128 (1984) 334–350. [15] Z.W. Lu, S.H. Wei, A. Zunger, S. Frota-Pessoa, L.G. Ferreira, Phys. Rev. B 44 (1991) 512–544. [16] A. Zunger, L.G. Wang, G.L.W. Hart, M. Sanati, Model. Simul. Mater. Sci. Eng. 10 (2002) 685–706. [17] J.M. Sanchez, Phys. Rev. B 81 (2010) 224202-1–224202-13. [18] C. Ravi, B.K. Panigrahi, M.C. Valsakumar, A. van de Walle, Phys. Rev. B 81 (2010) 104111-1–104111-13.
[19] A. van de Walle, M. Asta, G. Ceder, Calphad 26 (2002) 539–553. [20] A. van de Walle, M. Asta, Model. Simul. Mater. Sci. Eng. 10 (2002) 521–538. [21] A. van de Walle, G. Ceder, J. Phase Equilibria 23 (2002) 348–359. [22] H.T. Xue, F.L. Tang, X.K. Li, F.C. Wan, W.J. Lu, Z.Y. Rui, Y.D. Feng, Mater. Sci. Semicond. Process. 25 (2014) 251–257. [23] G. Kresse, J. Hafner, Phys. Rev. B 47 (1993) 558–561. [24] G. Kresse, J. Furthmuller, Phys. Rev. B 54 (1996) 11169–11186. [25] P.E. Blöchl, Phys. Rev. B 50 (1994) 17953–17979. [26] G. Kresse, D. Joubert, Phys. Rev. B 59 (1999) 1758–1775. [27] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865–3868. [28] A. van de Walle, G. Ceder, Rev. Mod. Phys. 74 (2002) 11–45. [29] A. van de Walle, G. Ceder, Phys. Rev. B 61 (2000) 5972–5978. [30] B.P. Burton, A. van de Walle, Chem. Geol. 225 (2006) 222–229. [31] O. Adjaoud, G.S. Neumann, B.P. Burton, A. van de Walle, Phys. Rev. B 80 (2009) 134112-1–134112-7. [32] P. Dalach, D.E. Ellis, A. van de Walle, Phys. Rev. B 82 (2010) 144117-1–144117-11. [33] B.P. Burton, A. van de Walle, U. Kattner, J. Appl. Phys. 100 (2006) 113528-1–113528-6. [34] M.K. Agarwal, L.T. Talele, Solid State Commun. 59 (1986) 549–551.