Structural and electronic properties of the random alloy ZnSe_xS_1-x

Next, we investigate the Zn-S and Zn-Se nearest-neighbor bond lengths for the first coordination shell. The variation in the average values of the nearest neighbor Zn-S (circles) and Zn-Se distances (squares) as a function of increasing Se concentration is shown in Fig. 2(a). The presence of two distinct distributions of bond lengths represent a strong deviation from simplistic expectations based on Vegard’s law. The latter should basically correspond to results obtained in VCA, which are reported in Fig. 2(b) as a black dashed curve. They show a perfect agreement with those obtained by converting the lattice parameters in effective average bond lengths (up triangles). This is done by imposing the optimized lattice parameters from the alloy systems in the ZnS wurtzite unit cell and then measuring the Zn-S bond lengths. Interestingly, this behavior would also be recovered by taking a full average of the Zn-$X$ ($X$ = S or Se) bond lengths, which will be analyzed below. Overall, Fig. 2(a) and Fig. 2(b) demonstrate that the individual Zn-S/Se bonds try to retain values close to those obtained in the pure end-members, without following Vegard’s law, which would instead dictate a much larger change.

More details on the distribution of the nearest neighbor bond lengths are obtained by examining the corresponding histograms, reported in Fig. 2(c). This picture shows clearly the presence of two distinct distributions, each with a spread of $\sim$ 0.05 Å around the average Zn-S/Se values. This spread is not small with respect to the bond length but is still much smaller than the difference between the average Zn-S/Se values. For reference, we provide a more compact overview of the data shown in these histograms, as well as in the following ones, in a series of Tables included in the SM SM .

The analysis of the previous data show that the average value of all the Zn-S and Zn-Se bonds present in the system together closely follows Vegard’s law (or equivalently what would have been predicted by the VCA). In addition, from Fig. 2(a), we can see that, though the average values of the Zn-S and Zn-Se bond lengths separately do not follow Vegard’s law, both of them slightly increase as a function of increasing Se concentration. Finally, the histograms in Fig. 2(c) show that there is a range of values about the average Zn-S and Zn-Se values. All these observations can be understood in terms of two extreme situations where we substitute a single Se atom in an S-rich environment and vice versa. The presence of a larger Se atom in an S-rich region creates a tensile strain locally due to a local chemical pressure effect. This tensile strain increases the distance between the surrounding atoms and as a result the Zn-S bonds about the substituted Se atom increase from their natural values. This effect increases the average value of the Zn-S bond distances in the system slightly. So we expect the average Zn-S bond length to increase gradually with an increasing Se concentration as shown in Fig. 2(a). Exactly the opposite mechanism takes place when we substitute a smaller S atom in a Se-rich region and hence, the average Zn-Se bond length gradually decreases with increasing S concentration. Now, in a realistic alloy system for any composition, there will not be just one type of environment but a fluctuation in the local composition. This fluctuation in the local composition will give rise to the dispersion in the extents of the volume expansion or contraction, leading to a distribution of Zn-S and Zn-Se bond distances about the average values as we see in Fig. 2(c). These features are very similar to those found for the cubic zinc blende phase as shown by the blue lines in Fig. 2(c). The average values of the Zn-S/Se bonds in the two phases match perfectly, which is consistent with our previous interpretation in terms of local environment only. A minor difference is due to the spread of the distributions, which is $\sim$ 0.04 Å in the cubic phase, i.e. about 0.01 Å smaller than in the wurtzite phase. This is expected, as the wurtzite phase possesses a lower symmetry than the cubic phase, which makes one of the four bonds in the tetrahedral cage slightly bigger than the other three. This is clearly visible in the histograms for pristine ZnS and ZnSe and results into a wider distribution of the bond lengths at intermediate compositions. This diminished symmetry also introduces an ambiguity in how the atoms are included in different shells of neighbors. For reference, our precise partitioning is reported in Table A1, for both wurtzite and zinc blende phases.

In the ideal wurtzite structure, the second coordination shell of any cation (anion) consists of 12 cations (anions). The only difference between the cationic and anionic sublattices in the case of the alloy systems under study is that, for the cationic sublattice, 12 atoms around each Zn atom are Zn atoms, but for the anionic sublattice both S and Se atoms are present around any S or Se atom. The distribution of the calculated Zn-Zn second neighbor bonds in the cationic sublattice is illustrated in Fig. 3(c). For all the intermediate concentrations, the distributions are clearly bimodal, albeit the peaks are less marked than those associated to the two separate nearest-neighbor distributions plotted in Fig. 2(c). This bimodality indicates the presence of a group of shorter and larger Zn-Zn bonds coming from the SZn₄ and SeZn₄ tetrahedral units present in the system. Furthermore, it allows us to calculate separate average values for the Zn-Zn distance coming from SZn₄ and SeZn₄ tetrahedral units, shown as respectively circles and squares in Fig. 3(a). In the same plot, we also report the composition weighted average of both type of Zn-Zn bonds (diamond), as a function of increasing Se concentration. The black dotted line represents the variation of a single hypothetical Zn-Zn bond that would have been predicted by VCA. We can see that the shorter and larger Zn-Zn bonds separately do not follow Vegard’s law, which suggests that the Zn-Zn distances around S/Se atoms in the alloy systems retain values close to their values in the pure systems. However, as for the first shell, the average value of all the second neighbor Zn-Zn bond lengths lies very close to the black dotted line corresponding to the VCA results.

We can then proceed to the analysis of the S-S, Se-Se, and S-Se second neighbor distances around the Zn atoms. The corresponding histograms, reported in Figs. 4(a), 4(b) and 4(c), are very different from the other plots we showed so far. For all the intermediate concentrations, the distributions of the bond lengths show a significant overlap, which is in sharp contrast with the two well marked peaks visible in Fig. 2(c). Nevertheless, since these distributions describe different pairs of atoms, it is still possible to calculate meaningful average values. Fig. 3(b) shows the average values of the S-S (up triangle), Se-Se (down triangle), and S-Se (right triangle) second neighbor bonds as a function of increasing Se concentration. The black dotted line again represents the variation of a single hypothetical average $X-X$ ($X$ = S or Se) bond that would have been predicted by VCA. From Fig. 3(b) we can see that the second neighbor bond lengths from the anionic sublattice follow Vegard’s law more closely compared to the second neighbor Zn-Zn - see Fig. 3(a) - or nearest neighbor Zn-S/Se bond lengths - see Fig. 2(a). This same trend can also be seen in the cubic case, where the average values as well as the range of the quantities described above closely match with the values obtained for the wurtzite system. For example, the distributions of the Zn-Zn bonds from the cubic phase are shown alongside that of the wurtzite phase in Fig. 3(c), with unfilled blue lines. We see a very good qualitative match between these two phases. The only minor difference between wurtzite and cubic structures is that the second shell of wurtzite ZnS and ZnSe consist of two different sets of bond with a difference of $\sim$ 0.01 Å, whereas in the cubic structure the 12 bonds are all of the same length. This is also visible from the histograms of the two end members in Fig. 2(c). As discussed above, these features are a result of the lower symmetry of the wurtzite phase, when compared to the zinc blende phase.

These observations in the solid solution $\text{ZnSe}_{x}\text{S}_{1-x}$ can only be understood by taking into account the variations in the cation centered $X$-Zn-$X$ and anion centered Zn-$X$-Zn ($X$ = S or Se) bond angles from their ideal value, obtained in pure ZnS or ZnSe. The ideal value here refers to the value of 109.47^o that is obtained for a regular tetrahedron, where all four faces are equilateral triangles. We first consider the cation-centered features. The variations of the average values of the S-Zn-S (circle), Se-Zn-Se (square), and S-Zn-Se (diamond) bond angles as a function of increasing Se concentration are shown in Fig. 5(a). For completeness, the histograms of their distributions are reported in Appendix B, for all compositions. The black dotted line in Fig. 5(a) corresponds to the ideal value of 109.47^o across the series. From this figure, we can see that the S-Zn-S angle increases gradually with increasing Se concentration from its value in ZnS. On the contrary, the Se-Zn-Se angle decreases with increasing S concentration. The S-Zn-Se bond angle values lie in the region between the S-Zn-S and Se-Zn-Se values. Next, we proceed to analyze the two anion-centered bond angles, whose distributions are also detailed in Appendix B. In Fig. 5(b), we report the variations of the average values of the Zn-S-Zn and Zn-Se-Zn bond angles as a function of increasing Se concentration. We can see that the average values of both angles remain almost unchanged with respect to a change of composition. These trends are also similar to those reported for the cubic phase and, as for the bond lengths, the only difference concerns the spread of the distribution. The cation-centered angles in the cubic system show a larger spread of $\sim$ 1^o compared to the wurtzite system. For the anion-centered angles, it is the opposite. Here the spread in the wurtzite system is $\sim$ 1^o larger compared to the cubic system.

Before providing a microscopic explanation of these results, it is useful to make a comparison of the present anionic solid solution with the cationic solid solution $\text{Zn}_{x}\text{Cd}_{1-x}\text{S}$ of Ref. SMukherjee_PRB_2014 . While the nearest neighbor Zn-S/Cd-S bond behave similarly to Zn-S/Zn-Se bonds when concentration is changed, the behavior of the second shell environment changes with respect to the different doping process, affecting the cationic sublattice instead of the anionic one. For $\text{Zn}_{x}\text{Cd}_{1-x}\text{S}$, the second coordination shell around any cationic site consists of both Zn and Cd atoms whereas any anionic site has only S atoms in its second coordination shell. This is just the opposite in the anion substituted $\text{ZnSe}_{x}\text{S}_{1-x}$ alloy system. As a result, in Ref. SMukherjee_PRB_2014 we see that the second neighbor Zn-Zn, Cd-Cd, and Zn-Cd bonds in the cationic sublattice follows Vegard’s law closely, similarly to what we report here for the S-S, Se-Se, and S-Se bonds in the anionic sublattice in $\text{ZnSe}_{x}\text{S}_{1-x}$. Concerning the bond angles, instead, we see that the average values of the cation centered S-Zn-S and S-Cd-S angles remain almost unchanged with respect to compositional change. It is the anion-centered bond angles, namely Zn-S-Zn, Cd-S-Cd, and Zn-S-Cd, that vary. With increasing Cd concentration, the average Zn-S-Zn angle gradually increases from its ideal value, whereas the average Cd-S-Cd angle decreases from its ideal value with increasing Zn concentration. The Zn-S-Cd angle lies in between Zn-S-Zn and Cd-S-Cd and behaves similar to the S-Zn-Se angle in the anion substituted system. This difference in the behavior of the second neighbor bonds between cation and anion substituted systems indicates that the structural changes at the atomic scale occurring in these two types of systems are very different, although the nearest neighbor cation-anion bonds follow the same trends.

The variation of the average second neighbor $X-X$ ($X$ = S or Se) bond lengths reported in Fig. 3(b) shows a behavior resembling what one would expect from Vegard’s law. Nevertheless, some deviations from it are still visible and therefore one may ask to what extent this behavior persists into the next shells of neighbors. In the case of pure ZnS/ZnSe in the wurtzite structure, the third shell around each Zn atom consists of 10 S/Se atoms, divided in three different types by symmetry (see Table A1). Two of those types have roughly the same bond lengths, which means that we expect to see only two different bond lengths. For example, the Zn-S bond in pristine ZnS has 1 short value of about 3.95 Å and 9 larger values of about 4.51 Å. Variations in the average value of short and long Zn-S (circle and square) and Zn-Se (down and up triangles ) bonds between Zn and S/Se atoms in the third coordination shell are illustrated in Fig. 6(a). The black dashed line indicates the variation of a single Zn-V bond between Zn and a composition averaged virtual atom (V) replacing both S and Se atoms, which would be predicted by VCA. Looking at short and long bonds separately, the three curves are on top of each other, which clearly shows that the third neighbor bond lengths follow Vegard’s law more accurately compared to those of the first and second neighbors. We stress again that the origin of these different bond lengths is not a different local environment, as in Fig. 2(a) or Fig. 3(a), but simply the result of the crystal symmetry of the wurtzite phase. This feature is, in fact, absent in the data for the cubic phase, where 12 equidistant anions are present in the third coordination shell around each Zn atom (see Table A1). The histograms showing the distributions of the Zn-S and Zn-Se bond-lengths are reported in Figs. 7(a) and 7(b). The two separate regions corresponding to the two different bond lengths in the third coordination shell are evident. For the end members, we can also observe the split of the long bonds into two types, with 3 and 6 elements. Overall, no significant variation is observed between the bonds Zn-$X$, depending on choice of $X$. The distributions have now a substantial width, spreading over about 0.3 Å. A similar spread is also observed in the zinc blend structure, as reported in the SM SM . It is somehow expected that the differences between wurtzite and cubic structures become less clear for pairs of atoms that are farther apart.

At last, we also investigate the variation in the atomic bond lengths arising from the fourth shell of neighbors. This shell is particularly interesting for our analysis, since no experiment has been able to probe these longer distances. The analysis of the trends reported for the first three shells suggests that a behavior in accordance with Vegard’s law should be expected. The data for the average Zn-Zn bond lengths as a function of increasing Se concentration, shown in Fig. 6(b), confirm this hypothesis. This establishes firmly that any technique probing over distances above the third shell of neighbors will provide results showing a linear trend with respect to the composition, at least for $\text{ZnSe}_{x}\text{S}_{1-x}$. This is the situation typical of diffraction experiments, which are used to affirm Vegard’s law. As expected, a single mode characterizes the distributions of the fourth neighbor bond lengths, shown in the histograms of Fig. 7(c). However, a substantial width of the distribution is still visible, which raises the interesting question on what coordination shell will show a narrow distribution that is indistinguishable from the behavior of the lattice constant. We compared the histograms of the distributions of farther and farther neighbors and observed a slow decrease of the width (data not shown). However, our analysis is probably limited by the size of our supercells, which provide an exact match of the pair correlation function of a random alloy only up to the $7^{\text{th}}$ coordination shell. Investigating larger supercells should provide a more detailed answer to this question, but this analysis is out of the scope of the present manuscript.

If we consider the behavior of the first nearest neighbor Zn-S and Zn-Se bond lengths, as well as the behavior of the second neighbor bond lengths and their bond angles, we can construct a model for the structural changes at the atomic scale taking place in the anion substituted $\text{ZnSe}_{x}\text{S}_{1-x}$ alloy system in either wurtzite/cubic phase SMukherjee_thesis_2018 ; TDan_thesis_2020 . In these alloy systems, where the substitution is done at the anionic sublattice, the SZn₄ and SeZn₄ tetrahedral units act as the basic building blocks and are connected via a corner shared Zn atom. This is a consequence of the almost invariant nature of the local Zn-S and Zn-Se bonds and anion-centered bond angles. Now, as a result of their different sizes, the SZn₄ and SeZn₄ units rotate/tilt with respect to each other, to accommodate and reduce the local strain. This results in a change of the $X$-Zn-$X$ bond angles, which we know, from Fig. 5(a), to vary as a function of composition. Thus, this leads to a much rapid increase or decrease in the $X-X$ second neighbor bond lengths, whose behavior becomes closer to what a microscopic analogue of Vegard’s law would predict.

The previous analysis allows us to model the details of the two extreme situations, where we substitute a single Se atom in an S-rich region and vice versa. Substituting a Se atom in an S-rich region is schematized in Fig. 8(a). The larger volume of the Se atom (green sphere) leads to a larger SeZn₄ tetrahedral unit, as indicated by the red triangle. The nearest neighbor Zn atoms (gray spheres) around the Se atom are displaced from their ideal positions, causing the associated S-Zn-S angles to increase. This has been shown by the small red arrows. Such movements of the Zn atoms further lead to a rotation of the surrounding SZn₄ tetrahedral units (gray triangles) as a result of the chemical pressure, depicted through black curved arrows. These cooperative rotations allow the S atoms (yellow spheres) inside the tetrahedral units to come closer to the Se atom, causing the S-Zn-Se bond angles to decrease. The initial state of the tetrahedral units before the substitution of the Se atom at the S site is indicated by the black dotted lines. An analogous scheme can be drawn for the situation where an S atom is substituted in a Se-rich region, shown in Fig. 8(b). A smaller SZn₄ tetrahedral unit (blue triangle) leads to a decreased Se-Zn-Se bond angles (indicated by the small blue arrows), forcing the surrounding SeZn₄ tetrahedral units to rotate in such a manner that the Se-Zn-S bond angles increase. So, from these two extreme pictures we understand that for an increasing Se concentration the S-Zn-S/S-Zn-Se bond angles shall gradually increase/decrease - see Fig. 5(a) - which leads to a more rapid increase in the second neighbor S-S distance with respect to the first neighbors, in agreement with Fig. 3(b). Similarly, for an increasing S concentration, the Se-Zn-Se/Se-Zn-S bond angles shall gradually decrease/increase - see Fig. 5(a) - ensuring a rapid decrease in the Se-Se second neighbor distance with increasing S concentration, in agreement with Fig. 3(b). From these structural considerations at the atomic scale, we can also observe an interesting competition in the structural relaxation. An increase in the S-Zn-S bond angle favors a decrease of the S-Zn-Se bond angle, whereas a decrease in the Se-Zn-Se bond angle favors an increase of the S-Zn-Se bond angle (we consider the S-Zn-Se and Se-Zn-S angles to be coincident here). These competing mechanisms lead to average values of the S-Zn-Se angle that lie in the region between the S-Zn-S and Se-Zn-Se values, exactly as depicted in Fig. 5(a). Furthermore, as the rotations of the tetrahedral units take place around the S/Se atoms, we expect the Zn-$X$-Zn (X = S or Se) bond angles to remain invariant when the the composition is changed, as visible in Fig. 5(b). This reveals unambiguously that the SZn₄ and SeZn₄ tetrahedral units try to constitute the basic and probably rigid building blocks of the alloy systems. Similar considerations based on the bond lengths and bond angles can be made for the cation substituted $\text{Zn}_{x}\text{Cd}_{1-x}\text{S}$ alloy system SMukherjee_PRB_2014 in the wurtzite phase, suggesting that the ZnS₄ and CdS₄ tetrahedral units can be identified as the basic building blocks. Hence, these observations from two very different alloy systems suggest that, in general, one cannot consider either the cation or anion-centered tetrahedral units as the basic building blocks. In the pristine systems, both units are equivalent. The difference between them occurs only through the alloying process that promotes one over the other as a basic unit, based on the nature of the substitution (cationic/anionic).

After having identified cation or anion-centered tetrahedral units as the basic building blocks in the alloy systems, we can proceed to analyze their effective rigidity. Here, by rigidity we mean that, in the case of the anion substituted system, the SZn₄ and SeZn₄ units can maintain their tetrahedral shapes in the alloy systems/solid solutions. This rigidity has been suggested in previous literature SMukherjee_PRB_2014 ; DD2 , on the basis of the variation of the average values of bond lengths and bond angles. However, the aforementioned analysis ignores the total picture that we get from their corresponding distributions and therefore may be misleading. As we described above, the variation of the average values can be understood from the picture of substituting a S/Se atom in the ZnSe/ZnS host in the dilute limit, but things can be very different for a realistic doping concentration, where disorder plays an important role. Fig. 3(c) illustrates clearly that, for a 50% concentration of Se, the Zn-Zn second neighbor bonds have a very broad distribution that cannot really be interpreted as bimodal. A similar absence of bimodality can be inferred from the analysis of the Zn-S-Zn or Zn-Se-Zn bond angles, as visible from Fig. B2 in Appendix B. This observation is important, as the short and long Zn-Zn bonds construct the tetrahedral units around the S and Se atoms respectively. Hence, the absence of bimodality in these distributions suggest that the SZn₄ and SeZn₄ units are not able to preserve their tetrahedral shapes in the alloy systems and cannot be considered as rigid. So, if they are not rigid enough to preserve their tetrahedral shapes, the next question is whether the SZn₄ and SeZn₄ units preserve their volumes or not. To answer this question, we calculated the volumes associated with all the SZn₄ and SeZn₄ units present in the system. Their distribution is illustrated in Fig. 9(a) for the three alloy systems we considered here. Since a clear bimodal distribution is observed, the SZn₄ and SeZn₄ units can be described as not so rigid, but volume preserving. The presence of volume-preserving distortions also explains the invariance of the average anion-centered bond angles, reported in Fig. 5(b). A very different scenario is observed for the cation-centered ZnX₄ units, which not only are non-rigid but do not even preserve their volumes. This is clearly visible in the histogram of the calculated volumes of the ZnX₄ units, shown in Fig. 9(b). The broad distribution observed in this case is most likely due to the fact that ZnX₄ units have five possible variants, ZnS₄, ZnS₃Se₁, ZnS₂Se₂, ZnS₁Se₃, and ZnSe₄, with all of them having different volumes. In conclusion, in this kind of alloy systems, the nearest neighbor cation-anion bonds always try to retain their value close to the values observed in the pure end members. The local chemical pressure resulting from the disordered substitution of a different atom at a given lattice site is compensated by a change in the bond angles leading to a rotation and tilt of the tetrahedral units that constitutes the basic building block of the system, which are SZn₄ and SeZn₄ units in our case. In the case of alloying at the cationic sublattice, cation-centered units act as the basic building blocks.

While Vegard’s law predicts a linear behavior for the evolution of the lattice parameters in a solid solution alloy, experimental measurements show a quadratic behavior for the band gap. The deviation from linearity is associated to a bowing parameter $b$, whose size increases when the difference in band gap between the end members increases JEBernard_PRB_1987 ; JWu_SSC_2003 . During the years, there have been a few attempts to provide a formula connecting the bowing parameter to relevant microscopic quantities, such as the mismatch in size, band gap or electronegativity, but only with limited success JEBernard_PRB_1987 ; JWu_SSC_2003 ; Mader_PRB_1995 ; Wei_PRB_1997 ; Shan_PRL_1999 ; Kent_PRB_2001 ; Ferhat_PRB_2002 ; Wu_PRB_2003 . In this section, we will perform a systematic analysis of various effects contributing to the bowing parameter, in the attempt of identifying what mechanisms play the most crucial role in $\text{ZnSe}_{x}\text{S}_{1-x}$.

We start by focusing on the atom-projected density of states for pristine ZnS and ZnSe in their DFT wurtzite optimized structure, reported in Fig. C1. The valence band is mainly formed by S/Se-$p$ states, plus smaller contributions due to Zn-$s$ and Zn-$p$ states. The conduction band is, instead, constituted in equal measure of Zn-$s$, Zn-$p$ and S/Se-$p$ states. Overall, Fig. C1 is in excellent agreement with previous reports ATorabi_JPCL_2015 ; Hussain2019MT ; SSapra_PRB_2002 ; RViswanatha_PRB_2005 ; SDatta_JPCM_2008 . Since the band gap is fully specified by Zn-$s$ and S/Se-$p$ states, we can focus on those for the analysis of the evolution of the density of states in the $\text{ZnSe}_{x}\text{S}_{1-x}$ alloy, illustrated in Fig. C2. We can see that doping does not induce any new states, but the band edges move continuously across the various compositions. This behavior is consistent with the typical experimental observations on isovalent alloys, although there are a few exceptions to this rule, as e.g. $\text{GaP}_{x}\text{N}_{1-x}$ and $\text{GaP}_{x}\text{Bi}_{1-x}$ JEBernard_PRB_1987 . To further verify that this picture holds in $\text{ZnSe}_{x}\text{S}_{1-x}$, we performed additional calculations for a single Se (S) impurity in a 128 atom supercell of ZnS (ZnSe). These results (data not shown) confirm that no detached impurity level forms, even in the dilute limit.

The aforementioned analysis of the spectra can be used to extract the band gap E_g. The values for pristine ZnS and ZnSe in their DFT optimized wurtzite structure amount to 2.08 and 1.18 eV, respectively. The experimentally reported band gap of ZnS in the wurtzite structure is in the range of 3.86 - 3.98 eV, whereas for ZnSe measurements indicate the value of about 2.87 eV MOshikiri_PRB_1999 ; Semi_Springer_1992 . This strong discrepancy between theory and experiment is due to the usage of a semi-local exchange-correlation functional, as the GGA-PBE, as already emphasized in previous works ATorabi_JPCL_2015 ; Hussain2019MT ; SSapra_PRB_2002 ; RViswanatha_PRB_2005 ; SDatta_JPCM_2008 . As shown below, a more accurate evaluation of the band gap can be obtained via more complex functionals, as e.g. meta-GGA or hybrid functionals ATorabi_JPCL_2015 . The relative variation of the band gap as a function of increasing Se content in the $\text{ZnSe}_{x}\text{S}_{1-x}$ solid solution is shown in Fig. 10(a). The deviation from linearity is evident, and fitting the data points to a quadratic equation leads to a bowing parameter b = 0.44 eV. This value is in very good agreement with previously reported theoretical values, for example, with the value of 0.46 eV reported in Ref. DLong_CMS_2018 by using a $\Delta$-Sol correction to DFT in GGA. A good agreement is also observed with respect to reported experimental values, which are in the range of 0.35 - 0.68 eV JEBernard_PRB_1987 ; JLu_Nanoscale_2012 ; AEbina_PRB_1974 ; MAAviles_InoChem_2019 . This large interval reflects variations in sample homogeneity and the phase admixture between wurtzite and zinc blende that may depend on the S/Se concentration MAAviles_InoChem_2019 . To verify if the error on the precise size of the band gap affects the bowing parameter significantly, we performed additional calculations with the mBJ exchange-potential in combination with LDA correlations MBJ_1 ; MBJ_2 . The predicted band gaps for pristine ZnS and ZnSe in the GGA-PBE optimized wurtzite structure amount to 3.70 and 2.68 eV, respectively. As expected, the agreement with the experimental values is significantly improved with respect to PBE. The variation of the band gap for increasing Se concentration is reported in Fig. 10(b). Fitting the data points to a quadratic equation leads to a bowing parameter b = 0.60 eV. This demonstrates that the choice of the exchange-correlation functional does not only affect the values of the band gap across the full range of compositions, but also changes the bowing parameter itself.

A better understanding of the nature of the bowing parameter in $\text{ZnSe}_{x}\text{S}_{1-x}$ can be achieved by identifying the most important contributions to it. For computational convenience, this analysis will be based on the GGA-PBE calculations, but our conclusions are most likely valid for other DFT functionals as well, as discussed at the end of this section. Our current understanding of the evolution of the band gap across different concentrations in binary alloys is based on the theoretical work of Ref. JEBernard_PRB_1987, , where the authors separate the effects due to volume change, electronegativity and structural deformation. The two steric mechanisms are those that we have already analyzed in the previous section, see e.g. Fig. 8. Starting from a ZnS environment, if we gradually replace S with Se, the larger size of the latter will drive a volume change of the tetrahedral building blocks. Moreover, the disorder of the dopants will induce a volume preserving structural distortion of the tetrahedra (see Section III.D), changing the cation-anion-cation or anion-cation-anion bond angles. These geometrical changes will modify the local electronic structure, by reducing the band gap. The third effect is associated directly to the increased electronegativity of Se with respect to S. In Ref. JEBernard_PRB_1987, , all these mechanisms were analyzed by extracting their corresponding contribution from a few limited calculations, at a fixed concentration. However, due to the larger computational resources available today, we are able to calculate directly the evolution of the band gap for the various supercells under well-defined constraints. This analysis is illustrated in Fig. 11(a), where the black circles represent the fully self-consistent data already shown in Fig. 10(a). The green circles represent the results of calculations where the structural deformation beyond the change of volume due to Vegard’s law is not allowed. The linear behavior shows unequivocally that the local structural deformation is the main factor contributing to the bowing parameter, while the remaining contributions, associated to volume deformation and electronegativity, are very small. This comparison is in agreement with the aforementioned studies, but provides also a new insight, i.e. the fact that the deviation from linearity arises mainly in the S-rich region, as shown by the more marked differences between black and green curves there. This feature may appear obvious at first, since we expect an asymmetric behavior due to the fact that adding a building block whose local electronic structure has a smaller band gap with respect to its host (as in the ZnS limit) decreases the system band gap, while the opposite situation (as in the ZnSe limit) does not lead to any increase. However, one may expect this asymmetry to be connected to the ionic nature of the substitution and hence it should not be fully canceled when structural deformations are suppressed. Since this asymmetry is not visible in the green curve, the main mechanism driving the asymmetry of the band gap bowing should have a different origin.

A better understanding of these data can be achieved by removing all the effects due to the different electronegativity. To this aim, we evaluate the band gaps of pure ZnS and ZnSe in the perfect wurtzite structure, but using lattice parameters obtained from Vegard’s law across the various concentrations. These data are illustrated by the two red curves in Fig. 11(a). As expected, a perfectly linear behavior is observed. We can now evaluate the role of the structural distortions while ignoring the chemical composition by calculating pure ZnS and ZnSe system in the structure corresponding to the optimized supercells of the solid solutions at each concentration. These data correspond to the two blue curves in Fig. 11(a). A substantial bowing is visible and can be quantified as $b$ equal to 0.24 eV and 0.17 eV for ZnS and ZnSe, respectively. We can see, then, that the same structural distortions in ZnS give rise to a larger bowing in comparison to ZnSe, albeit no asymmetry is noted with respect to the 50% concentration. Hence, the fact that the bowing parameter arises mainly in the ZnS region can be explained as a coexistence of the previous two effects, namely the asymmetric effect of inducing a local electronic structure modification in a host with a given band gap and the larger response of the electronic structure of ZnS to equivalent structural modifications. This also indicates that structural distortion as well as the difference in the atomic properties of S and Se atoms are necessary to account for the actual bowing. Finally, the previous analysis raises the question why the same structural distortion gives rise to a slightly larger bowing in the ZnS system compared to ZnSe. To provide an answer, we calculated the energy vs volume curves for both ZnS and ZnSe, shown in Fig. 11(b). The curve for ZnS is clearly stiffer than the one for ZnSe, indicating that inducing the same amount of local structural distortions costs more energy in the former than in the latter. The nearest neighbor Zn-S bonds and their corresponding bond angles are hence more rigid with respect to nearest neighbor Zn-Se bonds and corresponding angles. A more quantitative comparison of the stiffness of the curves in Fig. 11(b) can be obtained by calculating the bulk modulus $B$. Fitting the data by means of a Birch-Murnaghan equation of state Murnaghan ; Birch leads to values of 68.76 GPa for ZnS and 56.10 GPa for ZnSe. For reference, the experimentally reported values for the zinc blende phase amount to 77.1 GPa and 62.0 GPa, for ZnS and ZnSe respectively THomann_SSS_2006 ; OMadelung_springer_1987 . These data demonstrate that a similar amount of structural distortion is going to result in a larger change in the interaction between the cationic and anionic orbitals in ZnS, if compared to ZnSe. This difference is then reflected in the variations of the electronic structure induced by said distortion. Finally, we may speculate on how the presented analysis will depend on the choice of a different exchange-correlation functional. The bowing parameter is affected by the exchange-correlation functional directly, by modifying the band gap, and indirectly, by modifying the ionic positions. While the former effect may only bring quantitative changes, the latter may potentially change the decomposition illustrated in Fig. 11(a). However, this change is likely to be small, since GGA is very accurate in describing the structural properties, as shown by our analysis of pure ZnS and ZnSe. Previous calculations in LDA, where the largest discrepancy is expected, have corroborated this expectation by reaching the same conclusion as ours, i.e. that the local structural deformation is the main factor contributing to the bowing parameter JEBernard_PRB_1987 . Analogously, the asymmetry of the bowing parameter with respect to the $50\%$ concentration is not supposed to be affected by a different functional. In fact, we demonstrated that this asymmetry is a consequence of having different bulk moduli for ZnS and ZnSe, in spite of their precise value. This will be qualitatively predicted by most functionals.