A theoretical study of thermal stability and electronic properties of wurtzite and zincblende ZnOxS1−x

A theoretical study on the thermodynamic stability and electronic properties of wurtzite (WZ) and zincblende (ZB) ZnOxS1−x as been carried out with first-principles methods. Special quasi-random structures (SQSs) are employed to model the random alloys and from the calculated phase diagram, it is clear that the solubility limits in the mother lattices (sulfur-doped WZ-ZnO and oxygen-doped ZB-ZnS) are very small for both WZ and ZB phases. Detailed electronic and structural properties of quasi-ordered structures were investigated by examining alloy configurations in different supercells. In particular, we found that the concentration dependence of the band gap is highly nonlinear and the large bowing coefficients are also sensitive to the composition. The intriguing difference in band gaps and bowing coefficient between WZ- and ZB-ZnO1−xSx is evident. This result from first-principles calculations predicts that the band gap bowing with the WZ structure can be quite different from that with the ZB structure in some semiconductor alloys and deserves further experimental investigation for verification.


Introduction
ZnO is a promising candidate for optoelectronic and spintronic devices due to its unique physical properties, such as a large exciton binding energy (60 meV) [1]- [3]. Two important issues toward the development of ZnO-based ultraviolet optoelectronic technologies are (1) to engineer the band gap [4,5] and (2) to obtain the stable p-type conductivity [6]. For band gap engineering, the band structure can be modified by the substitution of isoelectronic anions or cations in ZnO. While a lot of experimental and theoretical work on cations-substituted alloys, such as Mg x Zn 1−x O, Cd x Zn 1−x O and Be x Zn 1−x O, has been carried out [7]- [12], relatively less attention has been paid to anion doping. One of the problems in the process of alloying ZnO x S 1−x originates from the significant difference between the stability of ZnS and ZnO at film growth temperature [13]- [15]. On the other hand, the large difference in the electronegativity between O and S is likely to result in a larger bowing parameter in ZnO x S 1−x [15]. Thus it would be possible to achieve a wide band gap modulation with a low fraction of S in ZnO.
In spite of the ongoing efforts, obtaining a stable p-type ZnO is still a challenging task [6,16,17]. Ashrafi et al [18] suggest that a zincblende (ZB) ZnO phase is expected to solve this long-standing problem. However, few reports have been found on ZB-ZnO epitaxy to date due to the lack of suitable substrates [19]- [22]. It was proposed that ZB-ZnOS/ZnS thin alloy layers can be used as interlayers and/or nucleation layers for the ZB-ZnO film epitaxy [19,20,22]. Furthermore, the ZnOS alloy layer can be used to improve ZnO device performance. For example, it has been shown that the use of the ZnOS buffer layer can effectively enhance the conversion efficiency of Cu(In, Ga)Se 2 -based thin film solar cells [23]. Pt contacts on S-treated n-type ZnO results in Schottky behavior in contrast to the Ohmic interface of untreated Pt/ZnO [24]. Recently, the valence-band offset bowing of ZnO 1−x S x was reported to enhance p-type nitrogen doping of ZnO-like alloys [25]. Previous theoretical calculations of ZnSO alloy mainly focused on ZB-ZnO 1−x S x alloy [15,26,27]. For example, Ito et al studied the phase stability of ZB-ZnO 1−x S x alloy with the valence force field (VFF) method [26]. Moon et al. studied the electronic properties of ZB-ZnO 1−x S x alloy by combining special quasi-random structures (SQSs) using first-principles calculation [15]. To our knowledge, there is no theoretical study about the ZB/WZ phase stability of ZnOS alloy by following the change of compositions (oxygen and sulfur).

3
The special role of ZnO 1−x S x alloy in further developments of ZnO-based optoelectronics provides essential motivation for the present work. We engaged systematic first-principles calculations to compare the thermal stability between wurtzite (WZ) and ZB structures. In particular, the composition-temperature (x-T) phase diagram is obtained to shed new lights on the growth process of ZnO x S 1−x . Moreover, we have also investigated the structural and electronic properties of the alloys with different phases. These fundamental physical properties will be important for designing and optimizing the growth of the related quantum wells and superlattices.

Computational methods and theoretical model
The importance of the local chemical environment such as charge transfer and local structural relaxation in determining the thermodynamic and electronic properties of alloys has received considerable attention over the past few decades. It has been shown that the effects of disorder and local correlation of the atomic arrangements on the stability and structure of alloys can be modeled with a moderate-sized supercell, because the most important physical properties can be analyzed correctly by considering the correlation of nearby neighbors [28]. The SQS method, developed by Zunger et al [29] in which the short and intermediate range correlations are considered, can be used to efficiently study the random alloy systems. In this work, an SQS with 32 atoms is employed for modeling the thermodynamics properties of WZ-and ZB-ZnO 1−x S x alloys. The SQS of this size was found to be effective when dealing with other alloy systems such as In x Ga 1−x N [30].
Modern epitaxial growth techniques (P-MBE, LMBE, RP-MOCVD) usually operate at non-equilibrium conditions. Therefore, samples grown are not able to reach thermal equilibrium and are likely to result in short range ordering. To study these alloy structures with quasiordering, we begin by constructing all configurations in three different 16-atom supercells. For the ZB structure, a 1 × 1 × 2 supercell (shown in figure 1(a)) is used. For the WZ structure, we have considered 2 × 2 × 1 and 2 × 1 × 2 supercells (shown in figures 1(b) and (c)). The 2 8 = 256 alloy configurations in each of the three supercells have been simplified by eliminating the symmetrically equivalent ones [10]. We found there are 27, 22 and 34 symmetrically distinct configurations in ZB-1 × 1 × 2, WZ-2 × 2 × 1 and WZ-2 × 1 × 2 supercells, respectively. The averaged electronic properties of these quasi-ordered structures can be analyzed by including the contribution of all the alloy configurations with the same concentration and using the following formula: where P i is the property of the ith configuration, g i is the degeneracy factor arising from the symmetry of the ith configuration, and i max is the number of symmetrically distinct configurations at the concentration x. Noted that the summation is limited to configurations at concentration x.
In the present work, the structural and electronic properties of all the alloy configurations described above are examined by the density functional theory (DFT) under the local density approximation (LDA) [30]- [32] along with the projector augmented wave (PAW) method [33]- [35] as implemented in the VASP program package [36,37]. The PAW potentials used for Zn treat 4p electrons in the valence. The k-space integral and plane wave basis, as detailed below, are tested to ensure the total energy that is converged at the 1 meV cation −1 level. The kinetic energy cutoff of 600 eV for the plane wave expansion is found to be sufficient, while different Monkhorst-Pack-type [38] k-point meshes are chosen for the different supercells. For the ZB-1 × 1 × 2 supercell, a 5 × 5 × 3 Monkhorst-Pack-type mesh is chosen in the first Brillouin zone. With regard to the WZ structures, the 3 × 3 × 4 and 3 × 5 × 3 Monkhorst-Packtype mesh are chosen for the WZ-2 × 2 × 1 and WZ-2 × 1 × 2 supercells. In order to enable the alloy configurations to converge to their most stable conformation, all atomic positions and lattice constants are fully relaxed. The total energies are obtained after the energy minimization is finished. The band gap at -point estimated after a self-consistent calculation that can generate high-quality charge density is finished with the optimized geometrical structure.

Basic properties of the parent materials
The LDA-calculated structural parameters of WZ and ZB (ZnO and ZnS) are listed in table 1.
All calculated values agree with the experimental ones within 2%. Hence, it is reasonable to expect that the lattice parameters of the alloys can be described with similar accuracy in our calculations. Usually, the larger ionicity will induce the smaller c/a ratio in the WZ  [48]. b Decremps et al [49]. c Madelung [44]. d Hellwege and Madelung [50]. e Kim et al [21]. f Ashrafi and Jagadish [18].
structures [18]. From the calculated c/a ratio, the Zn O bond has more ionic characteristics than the Zn S bond. It is well known that compounds with stronger polar bonds have a tendency to crystallize in structures with larger Madelung constants. Since the Madelung constant (MC WZ = 1.641) of the WZ structure is greater than that of the ZB structure (MC ZB = 1.638), it is reasonable that ZnO prefers to crystallize in the WZ structure and ZnS prefers to crystallize in the ZB structure. We found that WZ-ZnO is more stable than ZB-ZnO by 14.87 meV cation −1 and the reverse situation occurs for ZnS (ZB-ZnS is more stable than WZ-ZnS by 11.16 meV cation −1 ). The small differences in the total energy between WZ and ZB forms of ZnO and ZnS are due to the subtle differences between WZ and ZB structures. This is reasonable as ideal WZ and ZB structures are quite similar and only differ from the third nearest neighbor onwards.

Composition-temperature diagram by SQS
For alloy materials with non-isostructural parent compounds, the formation enthalpies of these SQS structures with different phases can be expressed as: where γ can be either WZ or ZB. The formation enthalpies of the SQSs are shown in figure 2.
We fit the data of the ZB structures to the form H = (α + βx) · x · (1 − x) + 14.87x and obtain α = 0.689 eV cation −1 and β = 0.007 eV cation −1 . For WZ structures, we found that . From the crossover of the two curves (shown in figure 2), we can roughly estimate that the phase transition point should be around x = 0.25. In addition, it can be found from the figure that the formation enthalpy curve of the ZB-ZnOS alloy is more symmetric than that of its WZ counterpart. At finite temperature, the thermal stability of the solid state system is determined by the Helmholtz free energy [39], where H is the mixing energy (formation enthalpy) and S is the mixing entropy. Within the mean-field approximation, the entropy is given by where k B is the Boltzmann constant. The binodal curve, which indicates the equilibrium solubility limits as a function of temperature, can be calculated by the common tangent approach from the Helmholtz free energy. As shown in figure 3, the critical temperatures, above which complete 7 miscibility is possible for some concentrations, are 4018 and 3831 K for ZB and WZ phases, respectively. These high critical temperatures reveal that it is difficult to mix randomly ZnO and ZnS irrespective of whether the alloy is synthesized as a ZB or WZ structure. For example, the solubility limit is less than 2% at growth temperature of ∼1000 K.
With the spinodal curve, we can draw out the region of metastablity. On the O-rich side, the ranges of the metastable region are relatively larger for B4 (WZ) phase than for B3 (ZB) phase. On the S-rich side, the reverse situation occurs. The concentration is more than 92% in the O-rich region for B4 phase and the limit of phase separation is 7% in the S-rich region for B3 phase at the growth temperature of 1000 K. The values are consistent with recent experimental results where the author achieved the ZB structure ZnO x S 1−x for 0 x 0.4 and the WZ structure ZnO x S 1−x for 95% x 1 in the bulk crystal [40]. So, combing the SQS structure with mean-field approximations, we can obtain properly the concentration-temperature phase diagram of ZnOS with B3 or B4 phase in the condition of the thermodynamic equilibrium. Here, we should point out that the contribution of phonons to the free energy was not included in the calculations due to large computational costs. According to a previous study on the phonon entropy in InGaN [30], the critical temperature should be reduced to ∼20-∼30% after the lattice vibration is taken into account. As a result of this, the experimental results fall in the vicinity of the stable region from the calculations and make the theoretical calculation more consistent with the experimental data. Additionally, it is found that the critical temperature still remains high for ZnO x S 1−x alloys. This means that the random mixing of ZnO and ZnS for all composition ranges is actually difficult under the thermal equilibrium conditions, whatever the WZ and ZB structures.

Quasi-ordered structures
In modern optoelectronic devices, the active media usually have the form of thin film grown by various epitaxial techniques. The non-equilibrium environments offered by these techniques will result in a wider range of the alloy concentration. As Yoo et al [13] reported, the range of concentration is 0.85 < x < 1 for the ZnO x S 1−x films by pulsed laser deposition. By RF reactive sputtering, Meyer et al [41] even extend the composition range to 0 x 1 in the temperature ∼620 K. The increase of the composition range can be attributed to the existence of the quasiordered structure induced by the non-equilibrium conditions.
In figure 4(a), we show the formation energy of all quasi-ordered configurations in the three supercells shown in figure 1. For ZB-ZnO x S 1−x , the metastable ground states are at the concentration of 25% with the Famatinite structure and at concentration of 50% with the Chalcopyrite structure. Both of these structures are layer structures oriented along the (201) crystallographic direction. This finding is consistent with the report by Liu et al [42] in which they have shown that isovalent ZB alloys can exhibit ordered coherent, tetrahedral networks characterized by alternate atomic bilayers along the (201) direction. For WZ-ZnO x S 1−x , we find that the alloy configurations in the WZ-2 × 2 × 1 supercell are in general more stable than that in the WZ-2 × 1 × 2 supercell and this appears to be a general trend in cation-substituted II-VI alloys [10]. The metastable ground states are at concentrations of 25, 50 and 75%, respectively.
Comparing these quasi-ordered structures with the random SQS structures, it is clear that some of quasi-ordered ground states have relatively lower formation energy, and therefore they are more likely to exist in samples prepared by epitaxial growth. At the concentration of 50%, the energy of the ZB Chalcopyrite structure is lower than that of the WZ ground state. This  means that in the range of 0 x 0.5, the alloy film will tend to form in the ZB structure under highly non-equilibrium epitaxial conditions. Our finding here is inline with the idea of using ZB-ZnO x S 1−x alloy as the nucleation layers for the ZB-ZnO film growth, since it is easy to grow ZB-ZnO x S 1−x in low temperature. In addition, the lattice constants of all quasi-ordered configurations are calculated and shown in figure 4(b). The lattice constants of ZnO x S 1−x alloy with WZ or ZB vary linearly with the change of the oxygen concentration, as indicated by the XRD experimental results [41]. This linear behavior governed by Vergard's rule shows that both phases will be structurally stable in Table 2. Band gaps of ZnO and ZnS with WZ and ZB structures. The first column shows the band gaps at -point the point calculated by the LDA method. Experimental data [44] are listed in the second column for comparison. The 'predicted' band gaps (shown in the third column) for those meta-stable forms that are obtained by shifting the LDA-calculated band gaps by the difference between the LDA-calculated and experimental values of the oxides in their stable phases. all ranges of concentration and therefore further supports the fact that using ZB-ZnO x S 1−x alloy as the nucleation layer will be beneficial to the growth of the ZB-ZnO film.

Electronic properties
As an optoelectronic material, the fundamental energy gaps of ZnO x S 1−x alloys deserve a more detailed investigation. The LDA-calculated band gaps of the binary constituents (ZnS and ZnO) with WZ and ZB structures are tabulated and compared with the experimental values in table 2.
The systematic underestimation of the LDA-calculated band gap is a well-known problem of the LDA [43,44], which fails quantitatively in dealing with excited-state properties. In order to compare with the experimental report, the LDA band gaps are 'corrected' by linear interpolation to fit the bulk values of WZ-ZnO and ZB-ZnS to yield the 'predicted' values, as suggested by Janotti et al [45]. The 'predicted' values of ZB-ZnO and WZ-ZnS by the above procedure are comparable with the values obtained in experiments. Since both ZnS and ZnO have fundamental energy gaps with the direct c (min) → v (max) transition, the LDA-calculated band gaps of each configuration are extracted by analyzing the -point energy levels. The original LDA calculated values are shown in figure 7(a). The band gaps of the ZB configurations are consistently larger than that of WZ configurations. This may be a result of the crystal-field splitting in the WZ structure due to its relatively lower symmetry. The band gaps of different alloy configurations with the same concentration can be quite different, for example, the band gaps of ZnO 0.5 S 0.5 can differ by more than 1 eV. This diversity in the band gap of alloy configurations signifies the importance of taking the alloy statistics into account when analyzing the electronic properties. Here, by considering the degeneracy factor of each calculated configuration in a supercell (shown in table 3) with formula (1), the averaged band gap is obtained as a function of concentration, as shown in figure 7(a).
The band gaps of the parent materials are used to correct the LDA-calculated band gaps of the alloy configurations. Quasi-particle corrections for the band gap have been shown to vary almost linearly with the composition in ternary nitride alloys [46]. It is therefore reasonable to expect that the same correction relation is valid for ternary oxide alloys. In order to match the experimental values for the band gap of the parent materials, we have therefore applied a linear Table 3. All the symmetrically distinct alloy configurations in ZB-1 × 1 × 2, WZ-2 × 2 × 1 and WZ-2 × 1 × 2 supercells. #(O) indicates the number of O atoms in the configurations and g j is the degeneracy factor defined in formula (1).  4  8  22  5  16  22  5  8  23  6  16  23  5  16  24  6  4  24  5  8  25  7  8  25  5  8  26  8  1  26  5  16  27  6  4  28  6  4  29  6  8  30  6  4  31  6  8  32  7  8  33  8  1 shift in the band gap for alloy ZnO x S 1−x . The corrected averaged band gaps and the experimental values are illustrated in figure 7(b). The corrected calculated band gaps of WZ-ZnO x S 1−x are consistent with the experimental values [41], except at very low oxygen concentrations. In this diluted O-doped limit, the bowing of the band gap is very large, as seen from the experimental data. This may be due to the large size and chemical mismatch of oxygen and sulfur. It should be noted that similar huge bowing of the band gaps have been observed in the As-rich impurity limit of GaN x As 1−x alloy as well.
Because of the importance of the fundamental band gap for device design, its nonlinear dependence on the concentration deserves further qualitative understanding and quantitative analysis. The linear behavior of the averaged lattice shown constant in figure 4(b) is due to the relaxation effect of the localized Zn O and Zn S bonds. At the same time, this relaxation effect will result in the positional disorder of the constituent atoms (Zn, S and O) due to the shift of atomic positions away from the sites of the parent lattice (WZ or ZB). Both volume deformation and positional disorder due to the relaxation effect will change the electronic structure. Furthermore, the compositional disorder will also result in the change of the electronic structure. If we consider the effects of relaxation and compositional disorder as a weak potential fluctuation, the change of the electronic structure can be evaluated by solving the self-energy equation associated with the scattering potential and it can be found that the change of the fundamental band gap will include a quadratic term bx(1 − x) with a positional bowing coefficient b [39]. So a quadratic function of the concentration is typically used to describe the downward shift of the band gap of isovalent semiconductor alloys with a constant bowing coefficient b as the fomula [47]: Both the calculated values and the experimental band gaps of alloy ZnO x S 1−x , shown in figure 7(b), indicate a very significant downward bowing of the band gaps. A constant bowing parameter is clearly not sufficient for describing the band gap. The nonlinearity of the bowing coefficient has, therefore, been fitted to the corrected averaged band gaps by b(x) = P 0 + P 1 x + P 2 x 2 + P 3 x 3 and shown in the inset of figure 7(b). For both WZ and ZB structures, the bowing parameters are relatively low in the center of the composition range and are largest for both phases in the S-rich limit. Furthermore, we found that the bowing parameter of WZ structures is systematically larger than that of ZB structures. While the electronic properties of ZB-ZnO x S 1−x have not been determined by experiments, it is clear from figure 7(b) that the electronic properties of ZB-ZnO x S 1−x are different from that of the WZ-ZnO x S 1−x . This result from first-principles calculation predicts that the band gap bowing with the WZ structure can be quite different from that with the ZB structure in some semiconductor alloys.
To shed some light on these puzzling results, we consider the structural properties of the ZnO x S 1−x alloy based on the configurations in ZB-1 × 1 × 2, WZ-2 × 2 × 1 and WZ-2 × 1 × 2 supercells. As shown in figure 5, the average internal parameter u just has a small change, following the change in the concentration of sulfur. This confirms that the WZ phase can be structurally stable in all ranges of concentrations, as the demonstration of the lattice constants. Nevertheless, the lattice parameter u from the Zn O bond is very different from that of the Zn S bond for each concentration. This means the diversity in the sizes of oxygen and sulfur ions can make a very big difference between the localized chemical environment of the Zn O bond and that of the Zn S bond and can result in a strong potential fluctuation, although the local strain due to the difference in size has been eliminated by charge transfer and structural relaxation. This may be the reason that the bowing coefficient b will change by following the variety of composition, whatever be the ZB and WZ structures. Lattice parameters u of all the symmetrically distinct alloy configurations in the WZ-2 × 2 × 1 supercell (represented by ∇) and the WZ-2 × 1 × 2 supercell (represented by ). Note that the blue triangles stand for the parameter u from Zn S bonds, the black triangles represent that from Zn O bonds, the red triangles depict the averaged parameter u for each configuration. The solid lines are the statistical averages obtained from equation (1). alloys, such as CdZnO. It is also found that the difference of the Zn O (or Zn S) bond length in both WZ and ZB crystal lattices is very small and it is impossible to induce the obvious difference of the bowing coefficient b of ZnO x S 1−x alloy with the WZ structure and that with the ZB structure. The Zn O and Zn S bond lengths along the 0001 direction and that in the (0001) plane with the WZ structure are shown in figure 6(b). It is found that the difference of the lengths of Zn O (or Zn S) bonds along the 0001 direction and that in the (0001) plane in ZnO x S 1−x alloy is more than that in pure crystals ZnO (or ZnS). This means the anisotropic crystal-filed splitting ( 0001 direction and (0001) plane) is larger in ZnO x S 1−x alloy than in pure crystals (ZnO and ZnS) with the WZ structure. So the difference of band gap bowing with WZ and that with the ZB structure may be attributed to the large anisotropic crystal-filed splitting in the WZ structure with the large chemical mismatch of oxygen and sulfur in ZnO x S 1−x alloy.

Summary
We have performed first-principles density-functional calculations to study the thermodynamic stability and material properties of semiconductor alloy ZnO x S 1−x with both ZB and WZ structures. The disordered alloy was modeled using the SQS approach and the formation enthalpy of the SQSs shows that the structural phase transition occurrs at x ∼0.25 due to the difference of the stable phases of the parent compounds (WZ-ZnO and ZB-ZnS). The large miscibility gap in the calculated x-T phase diagram suggests that at the typical growth temperature (∼1000 K), the solid solubility limits in mother lattices are very small (<2%) from both sides (sulfur-doped ZnO and oxygen-doped ZnS).
The quasi-ordered structures, likely to exist in epitaxially grown samples, are investigated by constructing all the alloy configurations in three distinct 32-atom supercells. We found that  some ordered structures are more stable than the random phase. In particular, ZB-ZnO x S 1−x structures are systematically more stable on the S-rich size, which implies that under epitaxial conditions, ZB-ZnO x S 1−x alloy can be used as good nucleation layers for the epitaxy of the ZB ZnO film. The electronic properties of the alloy ZnO x S 1−x are studied by analyzing the quasi-ordered alloy configurations. In both WZ and ZB forms, large bowing parameters are found to be asymmetric and composition dependent due to the large size and chemical mismatch. It is intriguing to note that both the band gaps and the bowing coefficient of WZ-ZnO x S 1−x are  [15] and the experimental data (magenta squares) [41]. The concentration dependence of the bowing coefficient is shown in the inset of (b).
significantly different from that of ZB-ZnO x S 1−x . Our prediction on the disparity between the electronic properties of the semiconductor alloys with WZ and ZB structures calls for further experiments for verification.