Ab initio thermodynamic model of Cu$_2$ZnSnS$_4$

Thin-film solar cells based on the semiconductor Cu$_2$ZnSnS$_4$ (CZTS) are a promising candidate for Terawatt-scale renewable energy generation. While CZTS is composed of earth abundant and non-toxic elements, arranged in the kesterite crystal structure, there is a synthetic challenge to produce high-quality stoichiometric materials over large areas. We calculate the thermodynamic potentials of CZTS and its elemental and binary components based on energetic and vibrational data computed using density functional theory. These chemical potentials are combined to produce a thermodynamic model for the stability of CZTS under arbitrary temperatures and pressures, which provide insights into the materials chemistry. CZTS was shown to be thermodynamically stable with respect to its component elements and their major binary phases binaries under modest partial pressure of sulfur and temperatures below 1100K. Under near-vacuum conditions with sulfur partial pressures below 1 Pa decomposition into binaries including solid SnS becomes favourable, with a strongly temperature-dependent stability window.


Introduction
Inorganic thin-film solar cells consist of several materials (a combination of metallic, semiconducting and insulating compounds) arranged in a particular order to exploit the photovoltaic effect. The deposition and optimisation of each layer requires a specific set of conditions and sometimes chemical treatments to 'activate' their performance. The role of compositional and structural variations in limiting performance for a range of solar cell technologies is known. 1 The development of processing and annealing conditions has been largely empirical in the past; however, the importance of chemical thermodynamics in this area is beginning to be recognised. 2 Cu 2 ZnSnS 4 (CZTS), a quaternary chalcogenide semiconductor, was introduced as a photovoltaic material in 1988 by Ito and Nakazawa. 3 Several crystal structures are known, but while early work assumed the stannite structure (space group I42m), the lowest-energy structure is now known to be the kesterite structure (space group I4, Fig 1). Computational work has shown this to be a few meV per formula unit lower in ground-state energy than the stannite and CuAu-derived structures. 4 All low-energy crystal structures are related to the face-centred cubic zincblende lattice, with Cu, Zn and Sn distributed over one sublattice and S filling a second sublattice.
In recent years CZTS and Se-containing variations have come under particular attention as a candidate for large-area thin-film cells, with a current record light-to-electricity conversion efficiency of 12.6 % in a Se-dominated cell. 5 The record efficiency for Se-free CZTS is 8.4 %. 6 Its distinct advantages over competing technologies are the combination of † Electronic Supplementary   a direct optical bandgap around the "optimum" 1.5 eV and its abundant, inexpensive elemental components. [3][4][5]7,8 The reserves and production rates of these materials suggest that CZTS is a strong contender for global-scale generation compared to peers including Cu(In,Ga)(S,Se) 2 , CdSe and CdTe. 9 The long-term requirements for such generation are expansive, forming part of "country-sized renewable facilities". 10 Large-scale production requires a pragmatic process, preferably one which is adaptable to "roll-to-roll" processing. Rapid reactions and modest pressures are therefore of particular interest, as is the avoidance of exotic and dangerous substances.
Laboratory studies have already demonstrated the complexity of the phase diagram, with secondary phases or partial disproportionation commonly observed, and off-stoichiometry compositions employed in order to manipulate this. [11][12][13][14] These are commonly expressed in terms of chemical potential, giving insight into the transitions between phases but not the corresponding physical conditions. In particular, the materials chemistry of the system is known to be sensitive to the partial pressure of the chalcogen atmosphere. This pressure has been manipulated by supplying S/Se solids, H 2 S gas and/or SnS(e) solids. [15][16][17][18][19] A recent paper provided Sn and Se powders with the intent that they would form reactive gas-phase Sn-Se compounds. 20 It should be noted that SnS is itself a semiconductor that has been attracting interest for application as an absorber layer in thin-film solar cells. 21,22 Weber et al. studied the relationship between the composition of CZTS and temperature, finding a significant shift in composition and loss of Sn at temperatures above around 500 • C under 10 −2 Pa of S, with some SnS evaporation at temperatures as low as 350 • C. 16 As well as two reactions involving a ternary phase, they propose the quaternary decomposition to binaries and sulfur vapours: where the SnS is lost to the vapour phase. Scragg et al. studied the SnS -S interaction experimentally and concluded with the aid of kinetic modelling that this is a two-step reaction in which solid SnS is formed in an equilibrium reaction, liberating sulfur vapours before evaporating to SnS vapour: 23 and SnS(s) SnS(g).
They found the partial pressure of S 2 to be critical in the region 10 −4 mbar (10 −2 Pa), predicting a stability envelope given sufficient SnS vapour to prevent irreversible evaporation.
A typical chalcogenide photovoltaic device consists of a metallic back electrical contact, an active p-type absorber layer, an n-type buffer layer and a transparent oxide electrical contact. For kesterite-based devices this is usually Mo|CZTS|CdS|In 2 O 3 :Sn. Evidence of chemical reactivity and compositional gradients is frequently found at each of the interfaces. 2 In particular, significant amounts of MoS 2 is formed at the Mo-CZTS interface during the annealing step. 24 It is known that such reactions change the currentvoltage characteristics of the solar cell, but the specific processes occurring and how to control them are not understood. There are significant opportunities for materials chemists to characterise the structure and properties of photovoltaic systems such as these.
While reference thermochemical data is available for wellstudied semiconductors such as ZnS and SnS 2 , for CZTS and related compounds such information is currently unknown. 25,26 For CZTS it is difficult to define a reference experimental sample; a recent study employing scanning transmission electron microscopy suggested substantial cation inhomogeneity in CZTS with features of the order 1nm, while atom probe tomography of CZTSe has revealed a co-existing ZnSe network with features of the order 10nm. 27,28 Significant cation disorder is expected in kesterites following high-temperature annealing due to the gain in configurational entropy, and controlled by largely-undocumented cooling rates. 29 In this study, we combine quantum and statistical mechanics to compute a range of thermodynamic potentials for pure kesterite CZTS and its elemental and binary components. This consistent set of data, with energies and vibrations computed using density functional theory (DFT), is used to predict the stability window for the material with respect to different processing conditions. The database is freely available and the model will be extended to include other materials systems and processing scenarios. ‡ It is not practical to examine the effect of long-range disorder with ab initio calculations, and if these are suspected to provide a significant thermodynamic driving force in CZTS formation then alternative approaches will be needed. The stannite phase was briefly examined and the data is included in the ESI † , but as the chemical potential is consistently within a few meV of the kesterite phase it does not meaningfully affect any equilibrium results. A configurational entropy term from inter-mixing is possible, but would require a greater understanding of any phase-coexistence.

Thermodynamic framework
Classical thermodynamics is used here to predict heats of formation and relative phase stabilities. A consistent approach is used to calculate the key thermodynamic potentials: internal energy U and enthalpy H. By calculating U as a function of temperature and pressure, the heat capacity C v and entropy S are also derived, and hence the Helmholtz free energy A and Gibbs free energy G. G may be seen as the 'key' to phase stability, as this potential is minimised at equilibrium.
Overall enthalpy changes ∆H are calculated from the molar enthalpies of componentsĤ i following Hess's law: where ∆n i is the stoichiometry change associated with component i. Likewise, Gibbs free energy changes are a sum of ‡ The project is hosted at http://github.com/WMD-Bath/CZTS-model and a current snapshot is included in the ESI.

Computational details
species chemical potentials µ i : These species-wise energies require a consistent reference point. Here the reference is a state in which all electrons are non-interacting (i.e., infinitely separated). We make the common assumption that chemistry is governed entirely by the electrostatic interactions and kinetic energy of electrons and nuclei, within the Born-Oppenheimer approximation. Gravity and nuclear forces are neglected. The chemical potentials can be separated into a ground-state contribution and a vibrational contribution; in this work: where a superscript '0' indicates conditions of absolute zero temperature and pressure. E DFT is the (athermal) ground state energy from density functional theory calculations, and E ZP is the energy from zero-point vibrations. Electronic excitations are neglected.
Expressions for these thermochemical potentials as functions of temperature and pressure, while aligning suitable reference energies, are at the heart of ab initio thermodynamics. [30][31][32][33] The forms used here were derived and applied in previous work by the authors 34 : for ideal gaseŝ while for incompressible solids θ is used to denote an intermediate reference state; thermodynamic properties are often provided relative to a standard temperature and pressure. This enables the use of standard heat capacities C p and C v , or tabulated enthalpies relative to some arbitrary state. (Note that C v and C p are used interchangeably for incompressible solids.) We emphasise that we cannot comment on the rates of reactions from our thermodynamic treatment. It is possible, for example, that kinetic barriers for a particular phase separation would be prohibitive for it to complete in the timescale of typical processing or annealing conditions.

Computational details
Ground-state total energies and forces were computed in DFT calculations with the FHI-aims quantum chemistry code. 35,36 These were used to optimise the initial crystal and molecular structures, without symmetry constraints, before modelling vibrational properties. The PBEsol exchange-correlation (XC) functional was employed; this functional employs the generalised gradient approximation (GGA) and is optimised for solid-state calculations. 37 Evenly-spaced k-point grids were used, with scaled sampling following the procedure of Moreno and Soler. 38 In general a 10Å reciprocal-space cutoff is sufficient for semiconductors and insulators, while a higher cutoff is helpful for achieving convergence in total energies for metals. The recommended 'tight' set of numerically-tabulated atom-centred basis functions was used throughout, except for the cases of Zn metal in which an extended set of 13 basis functions was employed, and Sn metal in which a full 'tier 2' basis set was used (17 functions per atom).
Vibrational calculations were performed with Phonopy ‡ , a code which implements the Parlinski-Li-Kawazoe "direct method" for computing solid-state phonons within the harmonic approximation. 39,40 From a primitive crystal structure, a mixture of finite displacements and analytical gradients is used to construct a dynamical matrix of force constants. Forces are obtained from DFT calculations on large periodic cells; generally it is necessary to form supercells in order to avoid self-interaction of displaced atoms. Supercell sizes for this study are listed in Table 1.
In this work the electronic structure was iterated until the analytical forces were converged to within 10 −5 eV/Å. Finite displacements of 0.01Å were used, and symmetry employed to reduce the number of calculations where possible. From the dynamical matrix, a set of frequencies is calculated, forming a vibrational model of the system within the harmonic approximation and defining E ZP . By applying Bose-Einstein statistics, a relationship is formed between temperature and vibrational energy, yielding U(T ), and hence C v (T ) and S(T ). Note that the effect of pressure is not taken into account with this method; a logical extension for exploring the anharmonic effects of lattice expansion would be the quasi-harmonic method or thermodynamic integration from molecular dynamics simulations with an appropriate ensemble. ‡ Phonopy is an open source code developed by Atsushi Togo from the earlier package FROPHO. It is available from

Crystal structures
Ternary phases and metal alloys were disregarded at this initial stage. The structure for CZTS was drawn from previous work, and optimised for the basis set and XC functional used in this study. 41 The initial structure of Cu was obtained from the Inorganic Crystal Structure Database (ICSD) in the form of a simple face-centred cubic cell (collection code 64699). 42 For Zn an initial structure was drawn from the ICSD (collection code 64990) consisting of a 2-atom hexagonal-close-packed unit cell. Standard local geometry-optimisation algorithms (also tested in the plane-wave code VASP) struggled to find the energy-minimising geometry as Zn is soft in the c-axis. A series of fixed c-values were tested over increments of 0.025 A, while relaxing the a and b parameters to find optimal a and c values; b, α, β , γ and the atomic positions were then fixed to the P6 3 /mmc space group. For β -Sn the lattice constant was drawn from the ICSD (collection code 40039) and a 2atom face-centred-cubic cell constructed. The α-S structure was based on a previous DFT study and is a large triclinic cell containing 32 S atoms; this in turn is a symmetry reduction from an orthorhombic conventional cell. 43 The structures of Cu 2 S phases have been previously studied using X-ray crystallography and density functional theory; while the high-temperature and cubic structures contain partially-occupied sites, low-Cu 2 S has a well-defined, albeit large, 144-atom monoclinic unit cell. 44 SnS 2 is the binary phase corresponding to the formal oxidation state in CZTS (i.e. Sn(IV)). A 3-atom hexagonal unit cell was obtained from the ICSD (collection code 100612). An 8-atom orthorhombic structure for the stable Pnma phase of SnS was drawn from recent work on the phase stability of this material. 43 ZnS is encountered in (and gives its name to) both the zincblende and wurtzite crystal structures. The starting point was 2-atom zincblende primitive cell with a lattice parameter a = 5.4053 A from reference data. 45 A corresponding set of calculations were performed for the wurtzite phase of ZnS; the results are not considered here as the zinc blende ground state is more stable and preferred for modelling, but the results are available as part of the ESI. † The initial and optimised structural parameters for all other phases are included in Table 2.

Sulfur vapours
The thermochemistry of sulfur has been studied relatively lightly given its abundance and the scale of application. While sulfur vapours are known to consist of a series of rings from S 8 down to the dimer S 2 , data for the intermediate rings and equilibrium mixtures is relatively scarce. Standard thermochemical tables prefer to treat the gas phase as an ideal diatomic gas, as do other treatments of the CZTS equilibrium. 2,23,26 How- ever, the vapour phase of sulfur is thought to contain a mixture of the cyclic allomorphs, and even above 1200 K may only be around 70% S 2 . [46][47][48] In the condensed phase, the stable structure α-S is a molecular solid formed of packed S 8 rings. In the absence of a comprehensive model, calculations here consider α-S solid and the S 2 and S 8 vapours; the true behaviour of the mixture is expected to be somewhere between the effect of these pure species. DFT-optimised structures were used for the ground state energies of the S 2 dimer and S 8 ring. Temperaturedependent data was drawn from the NIST-JANAF thermochemical data tables; these are based on spectroscopy and assume ideal behaviour. 26

Results
We begin by reporting the vibrational properties for each solid compound of interest. The resulting thermodynamic potentials are then combined to assess the formation of CZTS with respect to its constituent elements and isovalent binary sulfides. Finally the decomposition liberating SnS and sulfur vapour is considered.

Lattice dynamics
Optimised ground-state lattice parameters are given in Table 2. The only significant shifts in structure are for α-S (a large, soft, molecular crystal) and to a lesser extent the c parameters of Zn metal and SnS 2 , which are weakly bound. Following the procedure outlined above, phonon densities of states and band structures were computed for CZTS, Cu, Sn, Zn, α-S, Cu 2 S, ZnS, and SnS 2 . While there have been isolated reports Table 2 Lattice parameters for CZTS, elemental and binary precursors, before and after unit cell optimisation with the PBEsol functional. No symmetry constraints were enforced except for the case of Zn metal (discussed in Section 2.3). Except where other references are given, initial structures were drawn from the Inorganic Crystal Structure Database and collection codes are given in their respective discussions above. 42    The phonon dispersion for CZTS is shown in Fig. 2, while the remaining curves are available in the ESI. † The general behaviour is similar to previous work by Khare et al., but lacks LO-TO splitting at the Γ point, which was not considered in this case. 49 The 16-atom unit cell of kesterite results in 48 (3N) modes. There are two blocks of bands, which spread from 0 to 171 cm −1 and from 251 to 350 cm −1 . The calculated phonon density of states (DOS) (Fig. 2) shows activity between around 270 and 250 cm −1 ; from experimental studies the Raman spectrum is known to contain two distinct peaks at about 286 and 337 cm −1 (A 1 modes). 11 The associated total energies and thermal properties have been packaged into a Python code for ease of use and manipulation. These are also available as supplementary information, and as part of a continuing project online. †

Standard thermodynamic properties of
The standard formation enthalpy of kesterite CZTS at 298.15 K and under 1 bar of pressure is calculated to be −3.83 eV per formula unit (−369.1 kJ mol −1 ). At standard conditions the effect of temperature and pressure is small (< 1 meV) as there is no gas component to the reaction. However, stability is determined by the Gibbs free energy of formation, ∆G f . The free energy of formation is plotted against temperature and pressure in Fig. 3. The effect of pressure is negligible due to assumed solid incompressibility and absence of a gas phase. It is clear that CZTS is thermodynamically stable with respect to its solid elemental precursors over all reasonable processing conditions. Considering the equilibrium with sulfur vapours: 2 Cu + Zn + Sn + 2 S 2 (g) Cu 2 ZnSnS 4 (13) 2 Cu + Zn + Sn + 1 2 S 8 (g) Cu 2 ZnSnS 4 (14) we see a stronger interaction at high temperatures and low pressures (Fig. 4). This is driven by entropy, as it is more entropically favourable for sulfur to enter a low-pressure environment. The effect is greatest for S 2 , suggesting an instability at high temperatures over 1300 K at low pressures, whereas if only S 8 is to be considered then the formation appears irreversible even under relatively extreme conditions. Given that the actual composition of sulfur vapours is known to shift towards S 2 at high temperatures, Fig. 4a is more appropriate in this regime. 48

Stability of CZTS relative to binary sulfides
The binary sulfides are of interest both in terms of routes to forming CZTS and possible disproportionation reactions. Note that for a stoichiometric mixture, there is no dependence on the chemical potential of elemental sulfur: and hence the temperature-pressure dependence is again very mild (Fig. 5). However, the overall free energy change is con-  siderably reduced to the order −44 kJ mol −1 ; this is logical as the binary phases are themselves stable with respect to their elemental precursors. If we instead consider decomposition to Sn metal, due to an instability of SnS 2 , then a dependence appears as sulfur is released (Fig. 6) following the reaction Cu 2 ZnSnS 4 Cu 2 S + ZnS + Sn + S 2 (g).
This relationship is even stronger in the event of partial sulfur loss to form the divalent tin monosulfide, SnS: Cu 2 ZnSnS 4 Cu 2 S + ZnS + SnS(s) + 1 2 S 2 (g).
Note that Eqn. 17 appears identical to Eqn. 2; the only difference is that SnS is here defined as being the bulk solid, while it is understood that in the actual mechanism the SnS likely forms a reactive surface and may even be adsorbed to a CZTS or ZnS bulk phase. Application of the model at high (∼1 bar) partial pressures of sulfur suggest that such decomposition is not expected below temperatures of around 1300K. Figure 7 shows that CZTS becomes unstable following this process at around 600K and above for very low partial pressures of S 2 . The stability window is found to be almost identical to that predicted by kinetic modelling of experimental data (see figure 5 of Ref. 23). The agreement is especially interesting given that the kinetic model of Scragg et al. is based on SnS vapours (which were introduced as a gas stream in the accompanying experiment), while the result is reproduced here in an ab initio thermodynamic model with no such phase. 23 This model does not consider the evaporation of solid SnS, and as such is equivalent to the two-step model (Eqns. 2-3) given a saturated SnS vapour phase. It is reasonable to expect that the direct formation of SnS vapour would offer a higher entropy gain, lowering the free energy further and hence promoting decomposition; this would be equivalent to same two-step model with a very low partial pressures of SnS.
Direct comparison to experimental syntheses of CZTS is difficult as stability curves have not been measured directly, but it is possible to interpret experimental results where the conditions are reported clearly. At typical formation temperatures of 700-800 K we find that decomposition is only expected where there is a significant absence of sulfur in the atmosphere. Given that sulfur solids have a vapour pressure of the order 100 kPa in this temperature range, this would only be expected to occur where sulfur is limited or vapours are removed by a vacuum pump. 25 Ericson et al. successfully produced CZTS films by reactive sputtering followed by annealing in a static argon atmosphere of 35 kPa at 560 • ; given that they observed a correlation between sulfur loss and temperature during sputtering, we would assume that the sulfur was sufficiently mobile to form an equilibrium pressure during annealing. 19 Redinger et al. observed decomposition at 560 • C under "vacuum"; this is in agreement with our model provided that their vacuum pump maintained a sulfur partial pressure of around 0.1 Pa (10 −3 mbar) or less. 50

Conclusions
Based on first-principles total energy calculations, a thermodynamic model has been developed to describe the formation and stability of CZTS with respect to its elemental constituents and stoichiometric binary sulfides, as well as the tin monosulfide which is known to play an active role. Reactions involving solid and gaseous sulfur have been considered, the latter of which introduces a substantial temperature and pressure dependence. Temperature and pressure conditions have been related to the phase equilibrium; they indicate higher decomposition temperatures than those observed experimentally, but otherwise broadly similar behaviour. These results, which are ab initio except for the reference data for sulfur vapours, closely reproduce a previous model which was derived from experimental results and reference data. 23 It is clear that the pressure-temperature interaction is strong in nearvacuum conditions. This initial model is based on a number of approximations that could be removed in future work. These include the effects of thermal expansion and compressibility on the solid phases and the non-ideality of the S vapour. There is always a compromise between accuracy and computational cost, especially as future models will be extended to consider competitive ternary and quaternary phases, in addition to metallic alloys. One issue with rigorous validation of the model is the scarcity of experimental thermodynamic data so far. We suggest that high-temperature near-vacuum experiments need to control and report the annealing pressure carefully if they are to be reproducible and aid understanding of the phase equilibria of multi-component semiconductors such as Cu 2 ZnSnS 4 .