Density-functional theory simulations of CZTS, CZTSe, and CZTS0.25Se0.75 photovoltaic compounds have been performed to investigate the stability of the CZTS0.25Se0.75 alloy vs. decomposition into CZTS, CZTSe, and other secondary compounds. The Gibbs energy for vibrational contributions was estimated by calculating phonon spectra and thermodynamic properties at finite temperatures. It was demonstrated that the CZTS0.25Se0.75 alloy is stabilized not by enthalpy of formation but primarily by the mixing contributions to the Gibbs energy. The Gibbs energy gains/losses for several decomposition reactions were calculated as a function of temperature with/without intermixing and vibration contributions to the Gibbs energy. A set of phase diagrams was built in the multidimensional space of chemical potentials at 300 K and 900 K temperatures to demonstrate alloy stability and boundary compounds at various chemical conditions. It demonstrated for CZTS0.25Se0.75 that the chemical potentials for stability differ between typical processing temperature (∼900 K) and operating temperature (300 K). This implies that as cooling progresses, the flux/concentration of S should be increased in MBE growth to maintain the CZTS0.25Se0.75 in a thermodynamically stable state to minimize phase decomposition.

While more than 80% of solar energy in the world is generated by crystalline Si-based photovoltaic (PV) devices, thin film solar technology has great promise to gain a larger share in the PV market. Thin film PV devices have the advantage of using polycrystalline absorbers with direct-bandgap and short absorption lengths, which reduces the manufacturing costs and enables monolithic integration of the solar cells.1,2 Chalcogenide thin films, particularly CdTe and Cu(In,Ga)Se2, are well-developed PV technologies that are already commercialized and manufactured as PV modules with efficiencies greater than 10%.3 Scarcity of In and Te along with toxicity of Cd can impose serious limitations on the large-scale production of chalcogenide solar cells and, consequently, may require switching to different thin film materials.4,5 Therefore, in order to meet the global terawatt-scale energy demands, development of thin films with earth-abundant and non-toxic elements is essential.

Cu2ZnSn(S,Se)4 (CZTSSe) is the most promising earth-abundant material to replace thin film chalcogenides for solar cell applications because of its optimal direct, and tunable, bandgap energy (1.0-1.5 eV) and large absorption coefficient (>104 cm−1).6,7 However, the efficiency of CZTSSe champion cells is at least 6%-8% (absolute) lower than that for CdTe and CIGSe, mostly due to its low open-circuit-voltage (Voc).8 Further improvements in the PV performance of CZTSSe devices strongly rely on identification and passivation of the non-radiative carrier recombination sites that limit the Voc. Many of the previous theoretical and experimental investigations attributed the deficit in Voc to either potential fluctuations due to charged defects or to larger-scale inhomogeneities in the phase, crystal structure, or composition of the CZTSSe films.7,9,10

Phase and composition inhomogeneity of CZTSSe films stem from a very small stability region of single-phase kesterite in the equilibrium phase diagrams, which results in the formation of secondary phases.11,12 While the composition of secondary phases strongly depends on the growth conditions including temperature and partial pressure (or chemical potential) of CZTSSe constituents, ZnS and Cu2SnS3 are the two most widely reported secondary compounds for CZTS.13–15 Depending on the composition and electrical properties, many of the secondary phases that grow along with the CZTSSe films could limit the PV performance by acting as a barrier within the film or by forming interfaces with the kesterites with very large defect densities.16–18 Consequently, many theoretical studies have focused on finding the optimal conditions for synthesis of single-phase CZTSSe.11,19–22Ab initio DFT (density functional theory) calculations predicted that CZTS is stable with respect to its major binary compounds at temperatures below 1100 K under modest sulfur pressures.20 Similar ab initio calculations showed that, with the introduction of point defects including vacancies and antisites, formation of single-phase CZTS is unlikely.22 

Another important example of atomic scale compositional inhomogeneity is the intermixing between the Cu and Zn ions in the kesterite structure, forming the [CuZn + ZnCu] defect complex.11,23 Intermixing in combination with the low dielectric constant of CZTS-Se causes electrostatic potential fluctuations and the formation of band tail states that have been identified as one of the major factors limiting the Voc.10 Near-resonant Raman scattering studies of CZTS films found the critical temperature for the order-disorder transitions for Cu and Zn sites to be 533 K.24 In addition, neutron powder diffraction on CZTS films has illustrated a strong dependence of Cu/Zn intermixing disorder on the cooling regimen from the synthesis temperature to the room temperature, where slow-cooled films showed 70% order as opposed to quenched films with complete disorder.25 Therefore, it is critical to design the film synthesis and post-growth treatment processes in a temperature range that does not lead to the formation of significant amounts of secondary phases as well as disorder in the CZTSSe structure. In this study, by employing DFT calculations, the stability of CZTSSe alloys with respect to formation of pure CZTS and CZTSe as well as binary and ternary chalcogenides has been evaluated, considering not only enthalpies of formation, but also Gibbs energy contributions from thermal vibration and intermixing.

The DFT simulations were performed using the Vienna Ab Initio Simulation Package (VASP) plane-wave simulation package using projector augmented-wave (PAW) pseudopotentials (PP).26–31 The valence electrons for Sn included 5s25p2, for S—3s23p4, for Se—4s24p4, for Cu included 3d104s1, and for Zn included 3d104s2 states. The DFT simulations were performed with 500 eV plane-wave cutoff and 10−8 eV energy convergence level. The K-point grids used for all compounds are summarized in Table S1 of the supplementary material.32 The bulk compounds were relaxed at variable volume using the conjugate-gradient algorithm with 0.01 eV/Å force tolerance level, and their relaxed parameters are summarized in Table S2 of the supplementary material.32 The projection operators were evaluated in reciprocal space. For enhanced accuracy, an additional support grid for the evaluation of the augmentation charges was used. The additional support grid had 8× the number of points as the standard grid. The precision of the charge-density grids was set automatically to avoid aliasing errors (PREC = Accurate in VASP).

The relaxed variable volume unit cells were replicated to create perfect supercells which were later used to calculate force constants. The supercells sizes are summarized in Table S1 of the supplementary material.32 The force constants were calculated in VASP using the two-point central difference scheme with the finite-displacement step of 0.01 Å taking into account system symmetry. The obtained force constants were uploaded into the PHONOPY phonon calculation package to obtain phonon spectra and derive the thermodynamics properties (Helmholtz free energy, heat capacity, and thermal entropy) in a temperature range of 0-1000 K.33,34 The phonon spectra and thermodynamic properties were integrated with a 24 × 24 × 24 grid.

To provide an accurate representation of the CZTSxSe1−x semiconductor bandgap, the HSE06 hybrid exchange-correlation functional was used for some simulations.35–37 To make the bandgap representation even more accurate, the alpha (α) mixing parameter of the HSE06 functional that defines the mixing partition between the Perdew, Burke, and Ernzerhof (PBE) functional exchange and Hartree-Fock exchange terms was varied, the bandgaps were calculated and compared to the experimentally measured bandgap of 1.13 eV for CZTS0.30Se0.70 measured by Haight et al.38–40 The mixing value α = 0.23 was found to be the most optimal, providing bandgaps of 0.92 eV for CZTSe (vs. 0.96 eV for CZTSe,41 1.04 eV for CZTS0.03Se0.9742), 1.45 eV for CZTS (vs. CZTS values of 1.4-1.5 eV41,43–45), and 1.06 eV for CZTS0.30Se0.70 (vs. 1.13 eV40). The experimental samples used as a reference for the DFT calculations had high quality (as defined by the optical properties) to minimize possible microscopic defects, which potentially could affect the experimental bandgap. In addition, these experimental samples had been measured at room or even lower temperatures to minimize the effect of thermal expansion on the measured bandgap. Since the HSE06 fitted mixing parameter of 0.23 is very close to the well-tested default HSE06 mixing parameter of 0.25 and was fitted to the high-quality experimentally measured bandgaps, it provides one of the most reliable approaches with available experimental data. The more extended research investigating effect of HSE06 mixing parameter variation is provided elsewhere.46 All HSE06 simulations below were performed with a mixing parameter α = 0.23, which is very close to the default HSE06 mixing parameter of 0.25.

The CZTS0.25Se0.75 alloy was simulated with random placement of S and Se atoms in S sites of the CZTS kesterite unit cell which was relaxed at variable volume and later expanded to a 64-atom 2 × 2 × 1 supercell to calculate force constants. The supercell provides interatomic force constants over a larger real-space range, which increases the accuracy of dynamical matrix calculations and phonon spectra. Originally 5 different S/Se random placement configurations in CZTS0.25Se0.75 unit cell were relaxed choosing the configuration with the minimal total energy. All 5 configurations had close total energies since S/Se intermixing in CZTS0.25Se0.75 has low defect formation energy of 0.03 eV/pair. In addition to the S-Se alloy, pure kesterite CZTS and CZTSe unit cells were relaxed at variable volume and later expanded to 2 × 2 × 1 supercells to calculate force constants. To build phase diagrams, possible secondary phases should be simulated to delineate regions of stability. As secondary phases, supercells of Cu2SnS3, Cu2SnSe3, CuS, CuSe, SnS, SnSe, ZnS, and ZnSe were simulated. It is also possible that alloy secondary phases such as ZnS0.25Se0.75 and Cu2SnS0.25Se0.75 can be formed; however, since the alloys occupy such a small portion of the available phase space, the chances that the decomposition product alloys would have a phase space border with CZTS0.25Se0.75 are remote, and it would be even more improbable that the phase space border to CZTS0.25Se0.75 would persist in a significant temperature range since as shown below the chemical potentials for phase stability change with temperature. To obtain compound formation energies, pure bulk Cu, Zn, Sn, S, and Se elementary phases were also calculated. To minimize cell-size effects on the force constants required for phonon calculations, the supercells of secondary phases were constructed to have a total number of atoms as close as possible to the 64 atoms of the CZTSxSe1−x supercells and shapes close to cubic (Table S132). All these CZTSxSe1−x, secondary, and elementary phase unitcells were initially relaxed at variable volume below the 0.01 eV/Å force tolerance level with the PBE exchange-correlation functional and later were re-relaxed at variable volume with the HSE06 (0.23 mixing parameter) hybrid functional providing two datasets using PBE and HSE06 functionals.

The phase diagram can be built solving the appropriate choice of the following equations for CZTS0.25Se0.75 system stability:

2×μ(Cu)+μ(Zn)+μ(Sn)+4×μ(S)=ΔG(CZTS),
(1)
2×μ(Cu)+μ(Zn)+μ(Sn)+4×μ(Se)=ΔG(CZTSe),
(2)
2×μ(Cu)+μ(Zn)+μ(Sn)+μ(S)+3×μ(Se)=ΔG(CZTS0.25Se0.75),
(3)

where μ is the chemical potential and ΔG is the Gibbs energy of the corresponding compound. To require stability of CZTSxSe1−x compounds over secondary phases, the following relevant conditions for secondary phases should be satisfied:

μ(Cu)+μ(S)<ΔG(CuS),
(4)
μ(Cu)+μ(Se)<ΔG(CuSe),
(5)
μ(Sn)+μ(S)<ΔG(SnS),
(6)
μ(Sn)+μ(Se)<ΔG(SnSe),
(7)
μ(Zn)+μ(S)<ΔG(ZnS),
(8)
μ(Zn)+μ(Se)<ΔG(ZnSe),
(9)
2×μ(Cu)+μ(Sn)+3×μ(S)<ΔG(Cu2SnS3),
(10)
2×μ(Cu)+μ(Sn)+3×μ(Se)<ΔG(Cu2SnSe3).
(11)

In addition, all elemental components should favor formation of compounds instead of pure elements, which provides a 3rd set of conditions: μ(Cu) < 0, μ(Zn) < 0, μ(Sn) < 0, μ(S) < 0, μ(Se) < 0. The values of DFT calculated formation energies for PBE and HSE06 (α = 0.23) are presented in Table I. The differences between PBE- and HSE06-calculated formation energies result from the difference in the treatment of the exchange energy by PBE and HSE06 functionals. The local density approximation (LDA) and generalized gradient approximation (GGA) DFT functionals often underestimate the binding energy and over-estimate bond-length vs. experimental values. As a result HSE06, which is more accurate, provides larger (more negative) formation energies than PBE (Table I). Since formation energy is provided per formula unit, the accumulated difference is larger for big molecules like CZTS-Se, CZTS, CZTSe (all with 8 atoms/f.u.) than for smaller molecules like SnS, SnSe with 2 atoms/f.u.

TABLE I.

DFT calculated formation energies per formula unit with PBE and HSE06 (α = 0.23) functionals in comparison with other simulated values.

CompoundHSE06 (eV)PBE (eV)Other sources (eV)
CZTS −4.596 −3.758 −4.21 (Ref. 7
CZTSe −4.238 −3.268 −3.31 (Ref. 19
CZTS0.25Se0.75 −4.302 −3.368  
CuS −0.500 −0.402 −0.49 (Ref. 7
CuSe −0.427 −0.271 −0.30 (Ref. 19
SnS −0.987 −0.947 −1.01 (Ref. 7
SnSe −0.976 −0.906 −0.90 (Ref. 19
ZnS −1.870 −1.627 −1.75 (Ref. 7
ZnSe −1.674 −1.435 −1.45 (Ref. 19
Cu2SnS3 −2.627 −2.060 −2.36 (Ref. 7
Cu2SnSe3 −2.462 −1.789 −1.80 (Ref. 19
CompoundHSE06 (eV)PBE (eV)Other sources (eV)
CZTS −4.596 −3.758 −4.21 (Ref. 7
CZTSe −4.238 −3.268 −3.31 (Ref. 19
CZTS0.25Se0.75 −4.302 −3.368  
CuS −0.500 −0.402 −0.49 (Ref. 7
CuSe −0.427 −0.271 −0.30 (Ref. 19
SnS −0.987 −0.947 −1.01 (Ref. 7
SnSe −0.976 −0.906 −0.90 (Ref. 19
ZnS −1.870 −1.627 −1.75 (Ref. 7
ZnSe −1.674 −1.435 −1.45 (Ref. 19
Cu2SnS3 −2.627 −2.060 −2.36 (Ref. 7
Cu2SnSe3 −2.462 −1.789 −1.80 (Ref. 19

Fig. 1(a) presents the combined phase diagram for defect-free CZTS and CZTSe compounds at 0 K without any intermixing and thermal vibration contributions to the Gibbs energy. In this case the CZTS and CZTSe stability regions do not overlap, potentially eliminating formation of the CZTS0.25Se0.75 alloy. To double-check this instability without intermixing at 0 K, the formation energy of the CZTS0.25Se0.75 alloy was calculated considering CZTS and CZTSe as initial phases,

ΔHf(CZTS0.25Se0.75) = Etot(CZTS0.25Se0.75) − 1/4 × Etot(CZTS) − 3/4 × Etot(CZTSe), leading to near-zero positive values of +2.56 × 10−2 eV (HSE06) and +2.27 × 10−2 eV (PBE) consistent with ideal CZTS0.25Se0.75 formation being thermodynamically unfavorable at 0 K.

Fig. 1(b) presents the same combined phase diagram taking into account intermixing between Cu and Zn atoms lying in Cu-Zn planes where such switching is possible; as explained below both entropic and internal energy change effects were included. This has been experimentally shown by neutron powder diffraction measurements.19 The intermixing Gibbs energy contribution was calculated as ΔGm = − TΔSm = − T(kB × lnΩ), where Sm = intermixing entropy, T = temperature, kB = Boltzmann constant, and Ω = number of ways to arrange the atoms. The ΔGm per molecule was evaluated using Stirling’s approximation.

Since Cu-Zn intermixing creates [CuZn + ZnCu] switching defects with non-zero defect formation energy, the proper energy correction dU should be added to energy balance equations. The k × [CuZn + ZnCu] (k = 1, 2, 4) defects in CZTS0.25Se0.75 were recently investigated in more detail indicating that the lowest defect formation energy per pair was found for clustering defects.23 Five different configurations with 4 × [CuZn + ZnCu] were simulated in CZTS0.25Se0.75 with defect formation energies of 0.118, 0.133, 0.145, 0.163, and 0.183 eV/pair (HSE06) demonstrating low variation for different defect placement.23 The lowest defect formation energy +0.118 eV/pair (HSE06) was selected as the most probable value. The highest defect concentration is achieved when half of Cu/Zn pairs are switched, and the other half are left intact. Since the CZTS kesterite unit cell has 2 formula units, a dU energy contribution per formula unit is defined as dU = +0.118/2 = +0.059 eV (HSE06). The same defect was recalculated with PBE functional providing dU = +0.046 eV (PBE). Although the dU term was added for general completeness in Figs. 1(b), 2 and 3, it provides no visible change to phase diagrams due to its small value. Due to chemical similarity, the S and Se intermixing in CZTSSe has even lower defect formation energy of 0.03 eV/pair (HSE06). Therefore it was ignored for the dU term but was included into the Gibbs mixing energy term.

The inclusion of the Cu/Zn intermixing contribution to the Gibbs energy expands the CZTS and CZTSe regions of stability, creating a noticeable overlap between them and making possible alloy formation (Fig. 1(b)). For CZTS0.25Se0.75, intermixing between S and Se is also possible (as mentioned previously) and should be included into the intermixing Gibbs energy contribution thereby further expanding the region of alloy stability. Gibbs energy mixing contribution in the CZTS0.25Se0.75 alloy was evaluated at 900 K using Stirling’s approximation for different types of intermixing: (a) Cu/Zn mixing in Cu/Zn planes provides ΔGmix = −0.1075 eV/formula unit, (b) mixing of all S/Se provides ΔGmix = −0.1745 eV/formula unit, and (c) combined S/Se and Cu/Zn (in-plane) mixing provides ΔGmix = − 0.2820 eV/formula unit.

FIG. 1.

Combined phase diagram for CZTS (black region) and CZTSe (red region) (a) at 0 K with no intermixing and no thermal vibration contribution to Gibbs energy; (b) with included Cu/Zn (in-plane) Gibbs intermixing energy at 900 K.

FIG. 1.

Combined phase diagram for CZTS (black region) and CZTSe (red region) (a) at 0 K with no intermixing and no thermal vibration contribution to Gibbs energy; (b) with included Cu/Zn (in-plane) Gibbs intermixing energy at 900 K.

Close modal

At 0 K, there is no overlap between CZTS and CZTSe (Fig. 1(a)). Including the intermixing between Cu/Zn at 900 K (typical processing temperature) (Fig. 1(b)) creates an overlap between CZTS and CZTSe consistent with stability of CZTS-Se at 900 K but not at 0 K.

Besides intermixing, Gibbs energy contributions from ion thermal vibration at finite temperature should be taken into account. This term was estimated calculating phonon spectra and deriving the thermodynamic properties for CZTS0.25Se0.75 alloy and the secondary phases at a temperature range of 0-1000 K. Due to the high computational cost of HSE06 force constant calculations, which in the case of low alloy symmetry can require up to ∼2 × 3 × N (N = number of atoms) total energy calculations, the phonon simulations were performed with the more efficient PBE-GGA exchange-correlation functional.38,39 The obtained phonon spectra for CZTS0.25Se0.75, CZTS, CZTSe, and secondary phases (CuS, CuSe, SnS, SnSe, ZnS, ZnSe, Cu2SnS3, and Cu2SnSe3) are summarized in the supplementary material Figs. S3-S7.32 To verify convergence of phonon spectra vs. plane-wave cutoff, SnS unit cell relaxations and supercell phonon spectra simulations were repeated with 500, 600, and 750 eV plane-wave cutoffs; all produced the same geometry and practically identical phonon spectra as shown in the supplementary material Fig. S8.32 Using calculated phonon spectra, major thermodynamic properties such as Helmholtz free energy, heat capacity, and thermal entropy were obtained in a temperature range of 0-1000 K using standard principles of statistical thermodynamics.33,47 The calculated thermodynamic property curves are presented in the supplementary material Figs. S3-S7.32 Gibbs energy is defined as G = U − TS + PV (U-internal energy, T-temperature, S-entropy, P-pressure, V-volume). However, under typical conditions of photovoltaic growth and usage, the PV term is negligible (around 10−14 eV) and can be discarded, making it possible to use the calculated Helmholtz free energy value A = U − TS. This does not take into account the effect of thermal expansion on the phonon frequencies. The effect of thermal expansion potentially could be included applying quasi-harmonic approximation, however, it is much beyond the scope of this paper.

FIG. 2.

Phase diagram for CZTS0.25Se0.75 (black) (a) Gibbs mixing only, HSE06, T = 900 K; (b) Gibbs mixing only, PBE, T = 900 K; (c) Gibbs mixing and vibrational, PBE, T = 300 K; (d) Gibbs mixing and vibrational, PBE, T = 900 K; (e) the same as (c) but with CZTS and CZTSe added; (f) the same as (d) but with CZTS and CZTSe added.

FIG. 2.

Phase diagram for CZTS0.25Se0.75 (black) (a) Gibbs mixing only, HSE06, T = 900 K; (b) Gibbs mixing only, PBE, T = 900 K; (c) Gibbs mixing and vibrational, PBE, T = 300 K; (d) Gibbs mixing and vibrational, PBE, T = 900 K; (e) the same as (c) but with CZTS and CZTSe added; (f) the same as (d) but with CZTS and CZTSe added.

Close modal

The CZTS0.25Se0.75 phase diagrams calculated according to Equations (3)-(11) with added dS, dU terms, and μ(X) < 0 conditions are presented in Fig. 2. These phase diagrams are plotted as a function of Zn and Sn chemical potentials. The chemical potential of Cu was fixed at −0.2 eV for all 6 panes. The chemical potential of S for different panes had values of −0.5, −0.9, and −0.4 eV. The selected S chemical potentials roughly correspond to the middle point of S chemical potential range producing stable CZTS0.25Se0.75 at specific conditions. The chemical potential of Se can be derived from the other 4 chemical potentials. Figs. 2(a) and 2(b) compare the CZTS0.25Se0.75 phase diagrams with HSE06 and PBE including Gibbs mixing only since Gibbs vibration energy calculations with finite displacements are too computationally expensive with the HSE06 hybrid functional. All calculations with the Gibbs vibrational contribution were calculated with the PBE functional. Numerical analysis demonstrated that it was impossible to choose one value of S chemical potential with the available stable CZTS0.25Se0.75 phase for all 6 panes of Fig. 2, since intervals of S chemical potential for different panes with stable CZTS0.25Se0.75 had no overlapping interval. As shown below in the 3D phase diagrams (Fig. 3), while at higher temperature there is a larger range of S chemical potential over which CZTS0.25Se0.75 is stable, there is no value of S chemical potential which is common for the equilibrium regions of phase stability for CZTS0.25Se0.75 at both 300 K and 900 K. Although CZTS and CZTSe are not observed as secondary phases probably due to likely high kinetic barriers, the panes in Figs. 2(e) and 2(f) include CZTS and CZTSe limiting lines for general comparison.48,49 For the panes in Figs. 2(e) and 2(f), possible Cu/Zn intermixing in CZTS and CZTSe was taken into account. It is noted that these compounds slightly reduce the region of equilibrium phase stability both at 300 K and 900 K.

Fig. 2(a) presents the region of CZTS0.25Se0.75 stability calculated with the HSE06 hybrid-functional taking into account only the Gibbs mixing energy contribution at 900 K (without a thermal vibration contribution). Consistent with the overlap between CZTS and CZTSe in Fig. 1(b), there is a region of chemical potential over which CZTS0.25Se0.75 is stable at 900 K since the intermixing contribution to the Gibbs free energy decreases (becomes more negative and stable) with temperature. Fig. 2(b) presents the same region of stability but calculated with a PBE DFT functional. The transition from Fig. 2(a) to Fig. 2(b) shows that switching functional from HSE06 to PBE changes and shifts the region of stability.

All 6 panes of Fig. 2 have strictly vertical alloy boundaries with Cu2SnSe3 and ZnS secondary phases. In addition, the boundaries for Cu2SnS3 and SnS are horizontal for all 6 panes. The vertical and horizontal nature of these boundaries is easily derivable from chemical potential Equations (3)-(11) and formation energies when chemical potentials of Cu and S are fixed. For example, the vertical boundaries imply that, for some Sn chemical potentials at low Zn concentration, Cu2SnSe3 is formed and at high Zn concentration ZnS is formed at constant Cu and S concentrations.

Fig. 2(d) presents the CZTS0.25Se0.75 phase diagram with intermixing and thermal vibration contributions to the Gibbs energy at 900 K, a typical processing temperature. Comparison with Fig. 2(b) indicates that adding a thermal vibration contribution shifts the region of alloy stability in chemical potential space and changes adjacent secondary phases. As previously, the stability region is limited by Cu2SnSe3 and ZnS as a function of μZn chemical potential. However, the upper border (high μSn) is changed from SnSe to Cu2SnS3 (Figs. 2(b) and 2(d)).

Fig. 2(c) presents the CZTS0.25Se0.75 phase diagram with intermixing and thermal vibration contributions to the Gibbs energy at 300 K, a typical PV operating temperature. Comparison of Fig. 2(c) with Fig. 2(d) demonstrates that the region of CZTS0.25Se0.75 alloy stability is much larger at 900 K than at 300 K.

Analyzing CZTS0.25Se0.75 alloy neighbors in the phase diagrams, for all presented cases, the alloy stability region is limited from the left and right (regarding Zn chemical potential) by Cu2SnSe3 and ZnS, respectively (Fig. 2). This was confirmed by experiments of super-slow melt cooling.16,24,50 For the bottom boundaries, at lower Sn content and chemical potential, the alloy stability region is mainly limited by ZnSe (magenta) and CuSe (brown) secondary compounds, except in Figs. 2(b) and 2(c) where only ZnSe limits alloy formation for lower Sn chemical potential and except for the more hypothetical case in Fig. 2(e) where CZTSe forms a low boundary. For the upper boundaries, at high Sn content and chemical potential, the alloy stability region is limited by various secondary phases for different presented cases: by Cu2SnS3 (green) for the cases in Figs. 2(a) and 2(d), by SnSe (violet) for the cases in Figs. 2(b), 2(c), and 2(e), and by both CZTS (red) and Cu2SnS3 (green) for the not very realistic case in Fig. 2(f) due to likely kinetic limitations.48,49

FIG. 3.

Stability region of CZTS0.25Se0.75 alloy (red): (a) and (b) as a function of Zn, Sn, and S chemical potentials at μ(Cu) = − 0.20 eV and its projections (green and blue). (c) and (d) As a function of Cu, Zn, and Sn chemical potentials at μ(S) = − 0.50 eV and its projections. PBE functional. Mixing, thermal vibration entropy, and dU included.

FIG. 3.

Stability region of CZTS0.25Se0.75 alloy (red): (a) and (b) as a function of Zn, Sn, and S chemical potentials at μ(Cu) = − 0.20 eV and its projections (green and blue). (c) and (d) As a function of Cu, Zn, and Sn chemical potentials at μ(S) = − 0.50 eV and its projections. PBE functional. Mixing, thermal vibration entropy, and dU included.

Close modal

Figs. 3(a) and 3(b) present stability regions of the CZTS0.25Se0.75 alloy in 3D as a function of Zn, Sn, and S chemical potentials at 900 K and 300 K. The chemical potential of Cu is fixed at −0.20 eV, while the Se chemical potential can be derived from the other four and the alloy formation energy. Note that although chemical potentials can change, the correct stoichiometry of CZTS, CZTSe, and CZTS0.25Se0.75 is maintained by the chemical potential multipliers in Eqs. (1)-(11).

Figs. 3(c) and 3(d) present the same CZTS0.25Se0.75 alloy stability region but as a function of Cu, Zn, and Sn chemical potentials at the fixed μ(S) = − 0.50 eV. In this case Cu, Zn, and Sn chemical potentials are independent with the fixed μ(S) and dependent (derived) μ(Se). As previously stated, the stability regions were calculated for 300 K and 900 K temperatures. The 2-D projections provide visual representation of chemical potentials at alloy stability. Stability regions of Fig. 3 correspond to phase diagrams presented in Figs. 2(c) and 2(d).

In general, the phase diagrams presented in Fig. 2 and stability regions in Fig. 3 demonstrate that higher temperature significantly expands the alloy stability region against possible secondary phases. These diagrams also provide a guide to the experimentalists on how to adjust growth and cooling conditions to maximize the probability of forming CZTS0.25Se0.75 at high temperature and cooling without phase decomposition. The phase diagrams in Figs. 2(c) and 2(d) and stability regions in Fig. 3 can be used as a guide to choose the proper chemical environment for CZTS0.25Se0.75 synthesis and subsequent processing. In simple terms, the chemical potential has logarithmic dependence as a function of chemical concentrations and thereby related to flux in MBE with chemical potential of zero corresponding to the bulk concentration of each element. The main challenge is that CZTS0.25Se0.75 stability region in Fig. 3 is small and shifts in the multi-dimensional space of chemical potentials when the temperature is changed. Figure 3 reveals a major experimental challenge: there are no overlapping regions of stability of CZTS0.25Se0.75 with constant S chemical potential. This implies that if the system were in chemical equilibrium while cooling from 900 K to 300 K, there would be some secondary phase formation. However, due to diffusion barriers, it is possible that the single phase CZTS0.25Se0.75 could be kinetically frozen at lower temperatures.

The three dimensional phase diagrams suggest how to improve the probability of maintaining phase stability during cooling. For example, as shown in Figs. 3(a) and 3(b) for Cu = −0.2 eV , there is overlap in the regions of phase stability for points near μ(Zn) = − 1.4 eV and μ(Sn) = − 1.0 eV but there are no overlapping points for the S chemical potential. Instead at 900 K, the maximum stable point is μ(S) = − 0.85 eV while at 300 K the minimum stable point is μ(S) = − 0.45 eV. This suggests that as cooling progresses, the flux/concentration of S should be increased in MBE growth. Figs. 3(c) and 3(d) further illustrate the challenge of constant sulfur chemical potential if the material stays in equilibrium while cooling. Here the chemical potential of sulfur is fixed, and the regions of phase stability at 900 K and 300 K have no overlap in chemical potentials for Sn, Cu, or Zn. This implies that one would need to raise the flux/concentration of all three while cooling to main phase stability. It is noted that kinetic freezing may occur during cooling before the material leaves the region of phase stability; the phase diagram suggests the region which should have the highest probability of preventing phase instability in growth is near the maximum possible sulfur flux/concentration since this is closest to the region of phase stability at low temperature.

Several possible decomposition reactions were investigated by calculating the change of Gibbs energy as a function of temperature. Four possible decomposition reactions were investigated,

CZTS0.25Se0.7514CZTS+34CZTSe,
(12)
CZTS0.25Se0.7514[ZnS+Cu2SnS3]+34[ZnSe+Cu2SnSe3],
(13)
CZTSeZnSe+Cu2SnSe3,
(14)
CZTSZnS+Cu2SnS3.
(15)

The change of Gibbs energy ΔG = Gproducts − Greactants was plotted as a function of temperature taking into account (a) the G mixing term (without thermal vibration), (b) the thermal vibration term (without the mixing term), and (c) both mixing and thermal vibration terms. Note that enthalpy of formation was included in each curve for every compound, which makes these curves non-additive. The obtained curves are presented in Fig. 4; note these contain both the contribution of entropy and dU from intermixing.

FIG. 4.

Change of Gibbs energy for possible decomposition reactions as a function of temperature. The entropy and dU contributions are included.

FIG. 4.

Change of Gibbs energy for possible decomposition reactions as a function of temperature. The entropy and dU contributions are included.

Close modal

The decomposition reaction (12) is energetically unfavorable at full entropy (mixing + thermal vibration) at temperatures above 100 K (Fig. 4(a)). Below ∼100 K, CZTS0.25Se0.75 is thermodynamically favored to decompose by at most −0.02 eV. The low exothermicity and likely kinetic barriers at 100 K will make this decomposition difficult to observe experimentally. The curve for the same reaction with only the thermal vibration entropy contribution (green) shows minor decomposition favorability around −0.02 eV or less at higher temperatures.

The decomposition reaction (13) is energetically unfavorable at all temperatures in the range of 80-1000 K for all 3 combinations of entropy contributions (Fig. 4(b)). It shows limited stability less than 0.02 eV below 80 K with mixing and thermal entropy contributions. The CZTS0.25Se0.75 alloy has much lower symmetry due to random placement of atoms relative to much more ordered secondary phases in reactions (12) and (13). Since intermixing entropy contributions to Gibbs energy stabilize the alloy, it is logical that at higher temperatures CZTS0.25Se0.75 becomes more and more stable relative to more structurally simple and more ordered secondary phases.

The decomposition reaction (14) is energetically unfavorable in the 0-1000 K temperature range for all combinations of components: mixing only, “mixing + vibration” and vibration only entropy contributions (Fig. 4(c)). However at 900 K, the “mixing + vibration” stability is only 0.075 eV which is about 30% of the stability of CZTS0.25Se0.75. Since the stability is less than the thermal energy at typical processing temperatures, formation of secondary phases is almost impossible to avoid at equilibrium.

The decomposition reaction (15) is energetically unfavorable in the 0-1000 K range for all 3 entropy combinations (Fig. 4(d)). However, the thermal stability is less than 0.1 eV at typical processing temperatures which is ∼40% of the thermal stability of CZTS0.25Se0.75 making secondary phase formation far more likely at equilibrium for CZTS than CZTS0.25Se0.75.

The DFT-calculated phase diagrams for the CZTS0.25Se0.75 alloy have been determined with various exchange-correlation functionals as a function of temperature and chemical potentials. Formation of the CZTS0.25Se0.75 alloy vs. pure CZTS and CZTSe is stabilized not by enthalpy of formation, which has almost zero gain, but by entropy contributions to the Gibbs energy primarily from intermixing. In multidimensional chemical potential space, CZTS0.25Se0.75 is stable in a small region, which shifts as a function of temperature. It demonstrated for CZTS0.25Se0.75 that the chemical potentials for stability differ between typical processing temperature (∼900 K) and operating temperature (300 K). This implies that as cooling progresses, the flux/concentration of S should be increased in MBE growth to maintain the CZTS0.25Se0.75 in a thermodynamically stable state to minimize phase decomposition. As the temperature increases, the CZTS0.25Se0.75 stability region increases consistent with the stabilization effects of mixing. Several decomposition reactions of CZTS0.25Se0.75, pure CZTS and CZTSe to secondary phases were investigated, demonstrating that CZTS0.25Se0.75, CZTS, and CZTSe are stable at temperatures above 100 K when taking into account full thermal vibration and mixing entropy contributions; however, at typical processing temperatures CZTS0.25Se0.75 has more than 2.5× the thermodynamic stability of CZTS and more than 3× the thermodynamic stability of CZTSe. The presented DFT simulated data strongly correlate to experimental measurements.

The information, data, or work presented herein was funded in part by the U.S. Department of Energy, Energy Efficiency and Renewable Energy Program, under Award No. DE- EE0006334. The information, data, or work presented herein was funded in part by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, make any warranty, express or implied, or assume any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represent that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

1.
D. B.
Mitzi
,
O.
Gunawan
,
T. K.
Todorov
, and
D. A. R.
Barkhouse
,
Philos. Trans. R. Soc., A
371
,
20110432
(
2013
).
2.
C. A.
Wolden
,
J.
Kurtin
,
J. B.
Baxter
,
I.
Repins
,
S. E.
Shaheen
,
J. T.
Torvik
,
A. A.
Rockett
,
V. M.
Fthenakis
, and
E. S.
Aydil
,
J. Vac. Sci. Technol., A
29
,
030801
(
2011
).
3.
D. B.
Mitzi
,
O.
Gunawan
,
T. K.
Todorov
,
K.
Wang
, and
S.
Guha
,
Sol. Energy Mater. Sol. Cells
95
,
1421
(
2011
).
4.
T. K.
Todorov
,
K. B.
Reuter
, and
D. B.
Mitzi
,
Adv. Mater.
22
,
E156
(
2010
).
5.
C.
Wadia
,
A. P.
Alivisatos
, and
D. M.
Kammen
,
Environ. Sci. Technol.
43
,
2072
(
2009
).
6.
J. B.
Li
,
V.
Chawla
, and
B. M.
Clemens
,
Adv. Mater.
24
,
720
(
2012
).
7.
A.
Walsh
,
S. Y.
Chen
,
S. H.
Wei
, and
X. G.
Gong
,
Adv. Energy. Mater.
2
,
400
(
2012
).
8.
W.
Wang
,
M. T.
Winkler
,
O.
Gunawan
,
T.
Gokmen
,
T. K.
Todorov
,
Y.
Zhu
, and
D. B.
Mitzi
,
Adv. Energy Mater.
4
,
1301465
(
2014
).
9.
K.
Yu
and
E. A.
Carter
,
Chem. Mater.
27
,
2920
(
2015
).
10.
T.
Gokmen
,
O.
Gunawan
,
T. K.
Todorov
, and
D. B.
Mitzi
,
Appl. Phys. Lett.
103
,
103506
(
2013
).
11.
S. Y.
Chen
,
J. H.
Yang
,
X. G.
Gong
,
A.
Walsh
, and
S. H.
Wei
,
Phys. Rev. B
81
,
245204
(
2010
).
12.
I. D.
Oleksseyuk
,
I. V.
Dudchak
, and
L. V.
Piskach
,
J. Alloys Compd.
368
,
135
(
2004
).
13.
H.
Du
,
F.
Yan
,
M.
Young
,
B.
To
,
C. S.
Jiang
,
P.
Dippo
,
D.
Kuciauskas
,
Z. H.
Chi
,
E. A.
Lund
,
C.
Hancock
,
O. O. W. M.
Hlaing
,
M. A.
Scarpulla
, and
G.
Teeter
,
J. Appl. Phys.
115
,
173502
(
2014
).
14.
E. A.
Lund
,
H.
Du
,
W. M. H.
Oo
,
G.
Teeter
, and
M. A.
Scarpulla
,
J. Appl. Phys.
115
,
173503
(
2014
).
15.
J. J.
Scragg
,
T.
Ericson
,
T.
Kubart
,
M.
Edoff
, and
C.
Platzer-Bjorkman
,
Chem. Mater.
23
,
4625
(
2011
).
16.
X.
Fontane
,
L.
Calvo-Barrio
,
V.
Izquierdo-Roca
,
E.
Saucedo
,
A.
Perez-Rodriguez
,
J. R.
Morante
,
D. M.
Berg
,
P. J.
Dale
, and
S.
Siebentritt
,
Appl. Phys. Lett.
98
,
181905
(
2011
).
17.
J.
Just
,
D.
Luzenkirchen-Hecht
,
R.
Frahm
,
S.
Schorr
, and
T.
Unold
,
Appl. Phys. Lett.
99
,
262105
(
2011
).
18.
K.
Wang
,
B.
Shin
,
K. B.
Reuter
,
T.
Todorov
,
D. B.
Mitzi
, and
S.
Guha
,
Appl. Phys. Lett.
98
,
051912
(
2011
).
19.
S. Y.
Chen
,
A.
Walsh
,
X. G.
Gong
, and
S. H.
Wei
,
Adv. Mater.
25
,
1522
(
2013
).
20.
A. J.
Jackson
and
A.
Walsh
,
J. Mater. Chem. A
2
,
7829
(
2014
).
21.
J. M.
Skelton
,
A. J.
Jackson
,
M.
Dimitrievska
,
S. K.
Wallace
, and
A.
Walsh
,
APL Mater.
3
,
041102
(
2015
).
22.
P.
Sarker
,
M. M.
Al-Jassim
, and
M. N.
Huda
,
J. Appl. Phys.
117
,
035702
(
2015
).
23.
E.
Chagarov
,
K.
Sardashti
,
A.
Kummel
,
Y.
Lee
,
R.
Haight
, and
T.
Gershon
,
J. Chem. Phys.
144
,
104704
(
2016
).
24.
J. J. S.
Scragg
,
L.
Choubrac
,
A.
Lafond
,
T.
Ericson
, and
C.
Platzer-Bjorkman
,
Appl. Phys. Lett.
104
,
041911
(
2014
).
25.
S.
Schorr
,
Sol. Energy Mat. Sol. Cells
95
,
1482
(
2011
).
26.
G.
Kresse
and
J.
Hafner
,
Phys. Rev. B
47
,
558
(
1993
).
27.
G.
Kresse
and
J.
Hafner
,
Phys. Rev. B
49
,
14251
(
1994
).
28.
G.
Kresse
and
J.
Furthmüller
,
Comput. Mater. Sci.
6
,
15
(
1996
).
29.
G.
Kresse
and
J.
Furthmüller
,
Phys. Rev. B
54
,
11169
(
1996
).
30.
P. E.
Blochl
,
Phys. Rev. B
50
,
17953
(
1994
).
31.
G.
Kresse
and
D.
Joubert
,
Phys. Rev. B
59
,
1758
(
1999
).
32.
See supplementary material at http://dx.doi.org/10.1063/1.4959591 for detailed description of phonon spectra, calculated thermodynamic curves, and simulation parameters.
33.
A.
Togo
,
F.
Oba
, and
I.
Tanaka
,
Phys. Rev. B
78
,
134106
(
2008
).
34.
A.
Togo
and
I.
Tanaka
,
Scr. Mater.
108
,
1
(
2015
).
35.
J.
Heyd
,
G. E.
Scuseria
, and
M.
Ernzerhof
,
J. Chem. Phys.
118
,
8207
(
2003
).
36.
J.
Heyd
and
G. E.
Scuseria
,
J. Chem. Phys.
121
,
1187
(
2004
).
37.
J.
Heyd
,
G. E.
Scuseria
, and
M.
Ernzerhof
,
J. Chem. Phys.
124
,
219906
(
2006
).
38.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
39.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
78
,
1396
(
1997
).
40.
R.
Haight
,
X.
Shao
,
W.
Wang
, and
D. B.
Mitzi
,
Appl. Phys. Lett.
104
,
033902
(
2014
).
41.
S.
Chen
,
X. G.
Gong
,
A.
Walsh
, and
S.-H.
Wei
,
Appl. Phys. Lett.
94
,
041903
(
2009
).
42.
S.
Bag
,
O.
Gunawan
,
T.
Gokmen
,
Y.
Zhu
,
T. K.
Todorov
, and
D. B.
Mitzi
,
Energy Environ. Sci.
5
,
7060
(
2012
).
43.
M.
Ichimura
and
Y.
Nakashima
,
Jpn. J. Appl. Phys.
48
,
090202
(
2009
).
44.
H.
Katagiri
,
K.
Saitoh
,
T.
Washio
,
H.
Shinohara
,
T.
Kurumadani
, and
S.
Miyajima
,
Sol. Energy Mater. Sol. Cells
65
,
141
(
2001
).
45.
B.
Shin
,
O.
Gunawan
,
Y.
Zhu
,
N. A.
Bojarczuk
,
S. J.
Chey
, and
S.
Guha
,
Prog. Photovoltaics: Res. Appl.
21
,
72
(
2013
).
46.
J. E.
Moussa
,
P. A.
Schultz
, and
J. R.
Chelikowsky
,
J. Chem. Phys.
136
,
204117
(
2012
).
47.
J. W.
Ochterski
, “
Thermochemistry in Gaussian
,” http://www.gaussian.com/g_whitepap/thermo.htm, 17 June 2015.
48.
J. J.
Scragg
,
P. J.
Dale
,
D.
Colombara
, and
L. M.
Peter
,
ChemPhysChem
13
,
3035
(
2012
).
49.
W.-C.
Hsu
,
B.
Bob
,
W.
Yang
,
C.-H.
Chung
, and
Y.
Yang
,
Energy Environ. Sci.
5
,
8564
(
2012
).
50.
A.
Redinger
,
K.
Hönes
,
X.
Fontané
,
V.
Izquierdo-Roca
,
E.
Saucedo
,
N.
Valle
,
A.
Pérez-Rodríguez
, and
S.
Siebentritt
,
Appl. Phys. Lett.
98
,
101907
(
2011
).

Supplementary Material