Thermoelectric properties and figure of merit of Cu₂ZnSnS₄

Abstract

Thermoelectric (TE) power is a leading technology for addressing energy waste, and the quaternary sulphide Cu₂ZnSnS₄ (CZTS) is of considerable interest for lightweight and flexible thin-film devices. However, the large variation in the figure of merit, ZT, obtained in experimental studies highlights the need for a more comprehensive characterisation of the material properties. In this work, we apply a fully ab initio modelling approach to obtain high-quality reference predictions of the thermoelectric properties and ZT of CZTS. We find that CZTS has large bulk κ_latt, due to long phonon lifetimes, which can be minimised by nanostructuring. We predict significantly higher electrical conductivity, Seebeck coefficients and power factors with p-type (hole) doping compared to n-type (electron) doping. The electrical properties are less sensitive to crystallite size than the κ_latt, and we predict an industry-standard high-temperature ZT > 1 should be achievable for nanostructured p-type CZTS with crystallite sizes below 10 nm. Comparing our predictions to measurements also suggests that modelling heavily-doped CZTS or CZTS-based solid solutions may yield valuable insight into how the ZT can be optimised further.

---

1 Introduction

The world has historically responded to increasing demand for energy by burning fossil fuels. However, fossil fuels are a finite resource, and fossil-fuel derived energy produces greenhouse gas emissions that contribute heavily to anthropogenic climate change. There is therefore a pressing need to switch to clean and renewable sources of energy. Over 60% of the energy produced globally is lost as heat,¹ and improving the efficiency of energy-intensive processes is therefore a key enabler of the clean-energy transition.

Recovering waste heat to electricity using thermoelectric power is a promising approach to addressing energy waste.² Thermoelectric generators (TEGs) are solid-state devices that extract electrical energy from a temperature gradient.³ TEGs have no mechanical components and are thus easy to maintain, and do not produce any noise, heat, chemical waste or greenhouse gases during operation.^4,5 TEGs have proven applications in the aerospace industry, and can potentially be scaled to a wide range of applications from powering small devices (e.g. pacemakers), to improving the efficiency of automobile engines, to heat recovery from large industrial facilities.^6–8 The market for thermoelectric energy harvesting was valued at US $320 m in 2017 and is forecast to increase to $1.5 bn by 2028.²

TEGs utilise the Seebeck effect in a thermoelectric (TE) material, where a temperature difference induces a potential difference quantified by the Seebeck coefficient S:³

Positive and negative S indicate dominant p-type (hole) and n-type (electron) carriers, respectively, and a working TEG requires materials of both type to form a couple. The efficiency of a TE material can be defined by the dimensionless figure of merit ZT:²

where σ is the electrical conductivity, κ = κ_latt + κ_el is the thermal conductivity from lattice vibrations (phonons) and charge carriers, and T is the absolute temperature. The numerator S²σ is the thermoelectric power factor (PF), which is related to the maximum power output of the device. The ZT is optimised by maximising the PF and minimising the κ, and a ZT ≥ 1 is considered the minimum requirement for an industrially-viable TEG.^2,9

There is a strong correlation between the electrical properties S, σ and κ_el through the carrier concentration n, which makes maximising the ZT challenging.^8,10,11 Increasing n increases the σ, but may also increase the κ_el and degrade the S. The n should therefore be optimised to maximise the PF but without significantly increasing the electronic thermal conductivity. On the other hand, the κ_latt is largely independent of the electronic structure but is linked to the crystal structure, chemical bond strength, and structural dynamics.^8,12 All four properties in the ZT equation are also temperature dependent, and, in particular, the κ_latt typically falls sharply with T, resulting in many TE materials performing best at high temperature.

Current high-performance TE materials include PbTe, with a ZT of up to ∼2.2 at 915 K,¹³ and Bi₂Te₃, with a ZT of up to 1.35 from 300–400 K.¹⁴ Bi₂Te₃ is commonly used for the dual application of solid-state refrigeration, while PbTe has been used in aerospace applications.^15,16 Due to the scarcity of tellurium and toxicity of lead, which limit scalability, increase costs, and potentially cause problems with end-of-life disposal, there has been significant research into earth-abundant, non-toxic alternatives. Metal chalcogenide systems often show the required balance of favourable electrical properties and intrinsically-low κ_latt.^2,10,17 For example, orthorhombic SnSe has a high-temperature bulk ZT of up to 2.6,^18,19 which can be increased to ZT ≈ 3 in polycrystalline samples.²⁰ More complex quaternary chalcogenides have also shown potential as high-performance TEs due to the possibility of low κ_latt,^21–24 particularly when used in nanoparticulate or thin-film forms.^21,23

Among these, the quaternary metal chalcogenide Cu₂ZnSnS₄ (CZTS) is well known as a thin-film photovoltaic material,^25,26 and has recently attracted attention as a prospective TE due to its abundant constituents, low toxicity, and low cost.²⁷ However, studies on pristine CZTS have thus far indicated performance well short of the ZT ≥ 1 required to be industrially viable. For example, Yang et al. and Liu et al. reported maximum ZT of 0.026 and 0.039 at 700 K, respectively, for CZTS nanocrystals and consolidated powders.^12,28 This modest performance can be improved substantially through chemical doping. For example, Nagaoka et al. achieved a very large ZT = 1.6 at 800 K in Na-doped CZTS single crystals,²⁹ while Liu et al. found that substituting Zn with Fe and Cd in CZTS thin films yielded a significant improvement in the PF and an estimated maximum ZT of 0.18–0.69 at 550 K.³⁰ The latter result is particularly interesting given that thin-film TEs require less material, which reduces costs, and have the advantages of being lightweight and flexible.^31,32 From a performance perspective, thin films can also show high ZT through a reduction in the κ_latt due to the boundary scattering introduced by small crystallite sizes.^8,33

Despite these successes in preparing and characterising CZTS-based TEs, the considerable variation in experimental studies highlights the need for an improved understanding of the intrinsic TE properties of CZTS and the level of doping and nanostructuring required to optimise its performance. We have recently developed a fully ab initio workflow for predicting the thermoelectric properties and ZT, which we have tested against multiple classes of TE materials including chalcogenides and oxychalcogenides.^34–36 In this work, we apply our methodology to obtain high-quality reference predictions for the properties and ZT of pristine CZTS as a function of temperature and doping level. Our modelling indicates that p-type CZTS should show improved electrical properties and higher ZT compared to the n-doped material. We predict a relatively large κ_latt of 9.4 W m⁻¹ K⁻¹ at 300 K, which can be reduced by as much as 80% by nanostructuring to the ∼10 nm crystallite size reported in some experiments on thin films. On the other hand, the electrical properties are much less sensitive to nanostructuring, and we predict a maximum ZT, ZT_max, of 1.01 at 1000 K for p-type CZTS with a crystallite size of 10 nm. This is more than double the bulk ZT_max = 0.45, and can be increased by ∼25% with a further reduction of the crystallite size to 5 nm. This indicates that targeting thin films with small crystallite sizes and/or exploring other techniques such as defect engineering to minimise the κ_latt should provide a facile route to maximising the ZT of CZTS-based TEs.

2 Methodology

Calculations on bulk CZTS were performed using pseudopotential plane-wave density-functional theory (DFT) with the Vienna Ab initio Simulation Package (VASP) code.³⁷

Electron exchange and correlation were modelled using the PBEsol exchange–correlation functional.³⁸ The ion cores were modelled using projector augmented-wave (PAW) pesudopotentials^39,40 with valence configurations of: Cu – 3d¹⁰4s¹; Zn – 3d¹⁰4s²; Sn – 5s²5p²; and S – 3s²3p⁴. Explicit convergence testing determined that a plane-wave cutoff of 600 eV for the valence electrons, in conjunction with Γ-centered 3 × 3 × 4 and 4 × 4 × 1 Monkhorst–Pack k-point grids⁴¹ for the primitive and conventional cells, respectively, was sufficient to converge the total energy to <1 meV atom⁻¹ and the external pressure to <1 kbar.

The initial structure of CZTS in the kesterite phase was obtained from the materials project,⁴² symmetrised using spglib,⁴³ and fully optimised to tight tolerances of 10⁻⁸ eV and 10⁻² eV Å⁻¹ on the electronic total energy and ionic forces, respectively. When required, the following transformation matrix was used to convert from the conventional to the primitive cell:

The second-order and third-order force constants (harmonic/anharmonic FCs) were determined using the supercell finite-displacement method implemented in the Phonopy and Phono3py codes.^44,45 The harmonic FCs were computed from a 4 × 4 × 2 expansion of the conventional unit cell with 512 atoms and a displacement step of 10⁻² Å. An atom-projected phonon density of states (PDoS) curve was obtained by interpolating the phonon frequencies onto a uniform Γ-centered q-point mesh with 12 × 12 × 12 subdivisions and using the linear tetrahedron method for Brillouin zone (BZ) integration. A phonon dispersion (band structure) was obtained by interpolating the frequencies along strings of q-points passing through all the high-symmetry points in the tetragonal I4̄ BZ. The anharmonic FCs were determined in a 2 × 2 × 1 expansion of the conventional unit cell with 64 atoms and a step size of 3 × 10⁻² Å. The third-order FCs were used in conjunction with the harmonic phonon frequencies and eigenvectors to determine the phonon linewidths and lifetimes using perturbation theory.⁴⁵ The phonon lifetimes were then combined with the harmonic modal heat capacities and group velocities to calculate the lattice thermal conductivity within the single-mode relaxation time approximation (SM-RTA), with modal properties evaluated on a 12 × 12 × 12 q-point mesh.⁴⁵

The three electrical properties, viz. the Seebeck coefficient, electrical conductivity and electronic thermal conductivity, were calculated using the AMSET code.⁴⁶

An accurate electronic bandgap was obtained using the HSE06 hybrid functional⁴⁷ and used to correct an electronic band structure obtained using PBEsol on a uniform k-point mesh with 2× the number of points as that used for the geometry optimisation. This k-point sampling was subsequently increased by a factor of ∼20× using Fourier interpolation during the transport calculations.

The electron scattering rates were estimated for four mechanisms, viz. acoustic-deformation potential (ADP), polar-optic phonon (POP), ionised impurity (IMP) and piezoelectric (PIE) scattering, requiring the calculation of a number of additional materials properties. The acoustic deformation potentials were calculated from single-point calculations on deformed structures generated using AMSET. The high-frequency dielectric constant ε_∞ and atomic Born effective charges Z* were calculated using density-functional perturbation theory (DFPT).⁴⁸ The Γ-point phonon frequencies, the elastic constants C_ij, the ionic contribution to the dielectric constant ε_ionic, and the piezoelectric tensor elements e⁽⁰⁾_ij were determined using a combination of the DFPT and finite-difference routines in VASP. Finally, the static dielectric constant ε_S was then calculated as ε_S = ε_∞ + ε_ionic, and the phonon frequencies and Z* were used to calculate the polar-optic phonon frequency ω_po.

3 Results and discussion

3.1 Structure and phonon spectra

Kesterite CZTS crystallises in the tetragonal I4̄ spacegroup (no. 82) with sixteen atoms in the conventional unit cell (Fig. 1).

Fig. 1 Optimised structure of Cu₂ZnSnS₄ obtained in this work. Atom colors: Cu – blue; Zn – grey; Sn – orange; S – yellow. This image was prepared using the VESTA software.⁴⁹

Our optimized lattice parameters are compared to experimental measurements^50–53 and previous calculations^54–57 in Table S1 (ESI†). The calculated lattice constants a = b and c are within ∼1% of the measurements, and the calculated volume V is within 2–3%. Given that our optimized structures are “athermal”, obtained at 0 K without zero-point corrections, and that we would expect to see some degree of thermal expansion in experimental measurements performed at finite temperature, this agreement is more than acceptable.

Our lattice parameters are very close to those obtained in the calculations in ref. 55, which used the same PBEsol functional. The studies in ref. 54, 56 and 57 all used the local density approximation (LDA), which we would typically to expect to show a small “overbinding” compared to PBEsol. This is generally the case, except for the inexplicably large overestimate of the experimental lattice parameters reported in ref. 54. Ref. 54 also performed calculations with the HSE hybrid functional (it is not clear whether they used the HSE03 or HSE06 variant^47,58) and obtained similar lattice parameters to the present study.

The calculated phonon dispersion and atom-projected density of states (PDoS) of CZTS are shown in Fig. 2. The primitive unit cell of CZTS has n_a = 8 atoms, resulting in 3n_a = 24 phonon branches at each wavevector q in the phonon dispersion. All of the frequencies in the dispersion are real, confirming that the kesterite structure is dynamically stable, as expected.²⁵ The overall form of the DoS is comparable to that obtained in our previous calculations.⁵⁹ The dispersion curve is equally divided between bands of low/mid- and high-frequency phonon modes separated by a “phonon bandgap” of ∼2 THz. The PDoS shows that the two groups of modes are predominantly associated with motion of the metal and S atoms, respectively. The splitting in the dispersion can be explained by the difference in mass between Cu/Zn and Sn, and between the metal atoms and the chalcogen.⁶⁰

Fig. 2 Calculated phonon dispersion and density of states (DoS) of Cu₂ZnSnS₄. The DoS curve shows the total DoS in black and projections onto the Cu, Zn, Sn and S atoms in blue, grey, orange and yellow, respectively.

3.2 Lattice thermal conductivity

Using the single-mode relaxation-time approximation to the phonon Boltzmann equations, the lattice thermal conductivity is given by:⁴⁵

The macroscopic κ_latt is an average over contributions κ_qj from individual phonon modes qj with one of N_q wavevectors q and band indices j = 1 − 3n_a. The size of the κ_qj are determined by the heat capacities C_qj, group velocities ν_qj and lifetimes τ_qj.⁴⁵

We note that the κ_latt is a tensor, and the three diagonal components κ_αβ correspond to transport along the a = b (κ_xx = κ_yy) and c directions (κ_zz). To compare to experiments on polycrystalline films or pressed pellets with randomly-oriented crystallites, it is useful to compute the scalar average κ_latt, which we denote in Roman type to distinguish it from the tensor κ_latt:

The computed axial κ_αα and average κ_latt are compared in Fig. 3(a). We predict thermal conductivities of 9.43 and 2.24 W m⁻¹ K⁻¹ at T = 300 and 1200 K, respectively, the former of which is several times larger than orthorhombic SnS (κ_latt = 2.15 W m⁻¹ K⁻¹ at 300 K),⁶¹ Bi₂S₃ (0.9 W m⁻¹ K⁻¹)⁶² and Bi₂SO₂ (2.62 W m⁻¹ K⁻¹).³⁶ In contrast to the layered (oxy)chalcogenides SnS and Bi₂SO₂,^36,61 our calculations predict very minimal anisotropy in the κ_latt of CZTS, which might be attributed to the 3D tetrahedral bonding network.

Fig. 3 Calculated lattice thermal conductivity κ_αβ of Cu₂ZnSnS₄. Plot (a) shows the κ_αβ as a function of temperature, while plots (b) and (c) show the harmonic function (κ/τ^CRTA)_αβ and weighted-average lifetime τ^CRTA_αβ defined in eqn (6) in the text. On each plot we show the principal xx = yy and zz components of the tensors, corresponding to transport along the a = b and c directions, together with the scalar averages computed using the equivalents of eqn (5).

To analyse the κ_latt, we used the scalar-average modal contributions κ_qj to compute the cumulative thermal conductivity as a function of frequency at T = 300, 600 and 900 K (Fig. S1, ESI†). At all three temperatures, approximately 80% of the κ_latt is through modes with frequencies f_qj ≤ 2 THz, while most of the remaining 20% is through modes with 2 < f_qj ≤ 4 THz. With reference to the phonon spectrum in Fig. 2, this indicates that the majority of the transport is via the acoustic modes, with some contribution from the low-lying optic modes, and that the higher-frequency Cu/Zn- and S-based bands make a negligible contribution to the thermal transport.

To understand the origin of the high κ_latt of CZTS compared to other leading chalcogenide TEs, we used the constant relaxation-time approximation (CRTA) model developed in our previous work.^61,63,64 In this model, the κ_latt is written as the product of a harmonic function κ/τ^CRTA and a weighted-average lifetime τ^CRTA, which are indicative of differences in the group velocities and lifetimes, respectively:

We note that we have defined the τ^CRTA as tensors, which allows us to account for possible differences in the lifetimes of the dominant heat-carrying modes along the a = b and c axes.

The xx = yy and zz components of the κ/τ^CRTA and τ^CRTA tensors are shown together with the averages in Fig. 3(b) and (c). At T = 300 K we obtain an averaged κ/τ^CRTA of 0.439 W m⁻¹ K⁻¹ ps⁻¹, which can be compared to values of 0.718, 0.422 and 1.46 W m⁻¹ K⁻¹ ps⁻¹ for SnS, Bi₂S₃ and Bi₂SO₂.^36,61,62 We have previously shown that low κ/τ^CRTA can be linked to inhomogeneous chemical bonding, indicators of which are large primitive cells and low-symmetry spacegroups.⁶⁴ CZTS and SnS have the same-sized primitive cell (n_a = 8), and CZTS has a higher-symmetry spacegroup (tetragonal I4̄ vs. orthorhombic Pnma for SnS), which suggests the quaternary composition significantly reduces the group velocities. Indeed, the κ/τ^CRTA are comparable to Bi₂S₃, despite the fact that the latter has a larger, lower-symmetry primitive cell (Pnma, n_a = 20). On the other hand, the κ/τ^CRTA of Bi₂SO₂ are much larger, which we attribute to the smaller primitive cell (n_a = 5) and stronger ionic bonding.³⁴

Our calculated average lifetime τ^CRTA is 21.5 ps at 300 K, which is 7–12× longer than SnSe (3 ps⁶¹), Bi₂S₃ (2.14 ps⁶²) and Bi₂SO₂ (1.8 ps³⁶). This more than compensates for the low group velocities and produces a large overall κ_latt.

To establish the origin of the longer phonon lifetimes, we write the τ_qj as the product of an averaged three-phonon interaction strength P̃ and a weighted two-phonon joint density of states (w-JDoS) function N₂(q,ω):⁴⁵

where ω_qj is the angular phonon frequency and N₂(q,ω,T) counts the number of energy- and momentum-conserving scattering channels for a phonon with wavevector q and frequency ω at temperature T. (This is sometimes termed the “scattering phase space” in the literature.) The N₂ can be split into separate w-JDoS functions for collision (type 1) and decay (type 2) events and averaged over q to obtain a function of frequency only:

Fig. 4(a) shows the linear fitting^63,64 to obtain the P̃ at 300, 600 and 1200 K, and Fig. 4(b) shows the N̄₂ calculated at 300 K. We obtain a scaled P̃ × (3n_a)² = 9.9 × 10⁻⁹ eV², which is ∼1.5× weaker than the 2.5 × 10⁻⁸ eV² obtained in similar calculations on SnS.⁶¹ While it is difficult to compare the w-JDoS functions quantitatively, the N̄₂ of SnS peaks at 4.5 THz with a scaled value of N̄₂/(3n_a)² = 1.04 THz⁻¹,⁶¹ where the equivalent scaled N̄₂ for CZTS is around 2.5× smaller at 0.43 THz⁻¹. Based on these differences in P̃ and N̄₂, we would estimate CZTS to have 6× longer lifetimes, which is close to the 7× difference in the τ^CRTA. We can therefore attribute the longer lifetimes in CZTS compared to SnS to roughly equal contributions from weaker three-phonon interactions and a smaller scattering phase space.

Fig. 4 Analysis of the phonon lifetimes in Cu₂ZnSnS₄ using eqn (7) in the text. (a) Scalar average lattice thermal conductivity κ as a function of the three-phonon interaction strength P̃⁻¹ at T = 300, 600 and 1200 K, indicating the P̃ that reproduce the calculated κ. (b) Averaged weighted two-phonon joint density of states (w-JDoS) functions N̄⁽¹⁾₂ and N̄⁽²⁾₂ defined in eqn (8), calculated at T = 300 K.

A likely consequence of the long phonon lifetimes are long phonon mean-free paths (MFPs) Λ_qj:

Λ_qj = ν_qjτ_qj (9)

Long phonon MFPs would make the κ_latt highly sensitive to the crystallite size L. The effect of finite L can be modelled by adding a boundary-scattering contribution to the phonon lifetimes given by:

Fig. 5 shows the predicted κ_latt as a function of crystallite size at T = 300, 600 and 1200 K. At all three temperatures, the thermal conductivity drops sharply as L is reduced from ∼1 μm to 1 nm. At 300 K, we predict a κ_latt of 8.31 W m⁻¹ K⁻¹ for L = 1 μm, falling to 4.95, 1.59 and 0.3 W m⁻¹ K⁻¹ for L = 100, 10 and 1 nm. Based on these values, nanostructuring on the order of 10 nm could reduce the κ_latt by ∼80%.

Fig. 5 Dependence of the scalar average lattice thermal conductivity κ of Cu₂ZnSnS₄ on the crystallite size L at T = 300, 600 and 1200 K estimated using the boundary-scattering model in eqn (10) in the text.

The reported feature size of CZTS thin films is typically L < 1 μm,^65–70 and the sensitivity of the κ_latt to nanostructuring is clearly seen when comparing our predictions to measurements (Table 1). For example, ref. 12 reports a κ_latt of ∼1 W m⁻¹ K⁻¹ at 700 K for CZTS nanocrystals with L ≈ 10 nm. This is well reproduced by our calculations with L = 10 nm (1.1 W m⁻¹ K⁻¹), and corresponds to a 72% reduction of the predicted bulk κ_latt = 3.9 W m⁻¹ K⁻¹.

Table 1 Comparison of the calculated scalar average lattice thermal conductivity κ_latt of Cu₂ZnSnS₄ (CZTS), calculated using eqn (5) in the text, to experimental measurements on pristine CZTS samples.^12,71,72 For each experimental measurement we show the predicted κ_latt obtained with the reported grain or crystallite size, the crystallite size L for which the prediction best matches the measured κ_latt (where applicable), and the predicted κ_latt for bulk CZTS

	T [K]	Crystallite or grain size L	κlatt [W m^−1 K^−1]
a Determined from X-ray diffraction data.b Determined by fitting the measured κlatt to the Callaway model.
Exp.^71	294	24^a nm, 22.1/7.47^b nm	0.9
Calc.	300	5–20 nm	1.02–2.38
Calc.	300	4 nm	0.87
Calc. (bulk)	300	—	9.43
Exp.^71	294	24–28^a nm, 4.25/4.54^b nm	4.0
Calc.	300	60 nm	4.05
Calc.	300	5–30 nm	1.02–2.94
Calc. (bulk)	300	—	9.43
Exp.^72	623	800 nm	0.50
Calc.	620	800 nm	4.06
Calc.	620	2–3 nm	0.44–0.58
Calc. (bulk)	620	—	4.38
Exp.^12	700	10.6 ± 1.9 nm	0.97
Calc.	700	10 nm	1.11
Calc. (bulk)	700	—	3.87

---

However, there is some ambiguity in the experimentally-reported feature size, in that grains visible at longer length scales may be formed of smaller crystallites. Microscopy typically measures the grain size, whereas techniques based on X-ray diffraction or analysis of the κ_latt probe the crystallite size. The interfaces between the crystallites in a larger grain would impede phonon transport, and the L in eqn (10) can thus be equated to the crystallite size. A further complication is that stacking and twinning faults are common in CZTS, and, depending on the preparation conditions, can occur with densities corresponding to nm length scales.^73–75 Moreover, intrinsic defects, including e.g. vacancies and antisites, and strain promoted by small crystallites, as well as extrinsic defects introduced to optimise the charge-carrier concentration, could lead to enhanced phonon scattering and further reductions in the κ_latt. All of these features could in principle have a similar effect to the interfaces between crystallites, corresponding to nm-scale L in the boundary-scattering model that are well below the measured crystallite size.

With this in mind, ref. 71 reports considerable variation in the crystallite sizes determined using X-ray diffraction (XRD) and by fitting the κ_latt to the Callaway model, such that in one of the samples the Callaway model suggested an L around 5× smaller than measured by XRD. As suggested in the discussion, this likely indicates that the “effective” crystallite size for phonon transport is limited by the presence of defects. Our boundary-scattering model predicts a crystallite size of 4 nm for the as-deposited sample in this paper, which is comparable to, albeit lower than, the value obtained by fitting the κ_latt. On the other hand, we predict the ∼4.5× larger κ_latt obtained after sulphurization would be indicative of a 60 nm crystallite size, which is much larger than determined from both XRD measurements and from fitting the κ_latt. Given the off-stoichiometry of the experimental samples, we suspect this could indicate the formation of a secondary phase with higher thermal conductivity.

Small crystallites and/or the presence of defects could explain the κ_latt of 0.5 W m⁻¹ K⁻¹ reported for the hot-pressed powders in ref. 72. The very low thermal conductivity suggests much smaller crystallites than the reported 800 nm grains, and the L ≈ 2–3 nm estimated from the present modelling is more in keeping with results from other experimental studies.^12,72

In all cases, however, the calculated κ_latt for bulk CZTS is considerably larger than the measured thermal conductivities. Therefore, while the large impact of crystallite size and, potentially, of defects complicates the comparison between our predictions and experimental measurements, it is clear that the κ_latt in CZTS is highly sensitive to the (effective) crystallite size, and that optimising this, through whatever means, is likely to be key to mitigating the relatively high bulk thermal conductivity compared to other (oxy-)sulphide TEs.

3.3 Electronic structure and transport properties

The calculated electronic band structure and density of states (DoS) of CZTS, obtained using the HSE06 hybrid functional, is shown in Fig. 6. The conduction band minimum (CBM) and valence band maximum (VBM) both lie at the Γ point, giving a direct bandgap of E_g = 1.28 eV.

Fig. 6 Electronic band structure and density of states (DoS) of Cu₂ZnSnS₄ calculated with the HSE06 hybrid functional.

Experimental measurements of the bandgap range from 1.42–1.6 eV,^76–78 and our predicted value thus underestimates these by 0.14–0.32 eV (10–20%; Table S2, ESI†). Our prediction is, however, very close to the measured photoluminescence peak in the experiments in ref. 79 (1.3 eV). We also note that the E_g has been reported to be highly sensitive to stoichiometry.⁸⁰ Previous theoretical studies using hybrid DFT have reported values between 1.26–1.50 eV.^81–85 Our values are an excellent match to the corrected GGA calculation in ref. 85, with a difference of just 2 meV (0.16%), but show a larger discrepancy of 0.22 eV (14.7%) to the HSE06 result in ref. 83. However, with E_g > 1 we would expect limited thermal generation of (intrinsic) carriers, and thus that the transport properties would be dominated by the extrinsic carriers at typical doping levels above 10¹⁸ cm⁻³.^86–88 We therefore do not expect the underestimation of the bandgap relative to experiments to significantly impact our predictions of the transport properties.

The electrical conductivity, Seebeck coefficient and electronic thermal conductivity can be determined from the spectral conductivity Σ(ε,T) and generalised moments of the transport distribution function Lⁿ(ε_F,T) given by:⁴⁶

e is the elementary charge, k and j are the electron wavevectors and band indices, respectively, and ν_kj, τ_kj and ε_kj are the electron group velocities, lifetimes and band energies. f ⁰(ε,T) is the Fermi–Dirac distribution function given by:where ε_F is the Fermi energy and k_B is the Boltzmann constant. The σ, S and κ_el are obtained from the Lⁿ as:

σ(ε_F,T) = L⁰(ε_F,T) (14)

The τ_kj can be computed from first-principles using electron–phonon coupling calculations. However, these are relatively expensive, and prohibitively so for complex materials with large unit cells. For semiconductors, similar accuracy can often be obtained by estimating the electron scattering rates (inverse lifetimes) within the momentum relaxation-time approximation (MRTA) using phenomenological models for the scattering matrix elements.⁴⁶ In the present work, we use the expressions for the scattering rates Γ_kj due to acoustic deformation potential (ADP), polar-optic phonon (POP), ionised impurity (IMP) and piezoelectric (PIE) scattering defined in ref. 46. The overall τ_kj are then calculated by combining the four rates using Matthiessen's rule:

τ_kj⁻¹(T) = Γ^ADP_kj(T) + Γ^POP_kj(T) + Γ^IMP_kj(T) + Γ^PIE_kj(T) (17)

This approach has been shown to give good predictions for a variety of semiconductors,⁴⁶ and has also been carefully validated against several well-studied families of TE materials including the Sn/Ge chalcogenides^34,35 and the bismuth oxychalcogenides.³⁶ We note that the parameters used in the scattering models are determined from first-principles calculations and are not, for example, inferred from similar materials or fitted from experimental measurements.

Under the assumption that doping to adjust the carrier concentration does not affect the host band structure, the ε_F can be adjusted to include a nominal extrinsic carrier concentration, in addition to the thermally-generated carrier populations, allowing the properties to be predicted as a function of n and T. This is termed the “rigid band approximation” (RBA). While RBA calculations can provide a useful guide for the carrier concentrations required to optimise the electrical properties,^34,36,89 the approximation can break down at large n due to modulation of the electronic structure induced by heavy doping. However, to investigate this explicitly requires calculations on atomistic models with dopants, and possibly a more sophisticated treatment of the electron–phonon coupling to compute the τ_kj, which we consider beyond the scope of the present study.

Finally, we also note that the three electrical properties are tensors, and, as for the κ_latt, the average of the three diagonal elements are most representative of pressed pellets and thin films (cf. eqn (5)). We denote these quantities in Roman type as σ, S and κ_el respectively.

Fig. 7 shows the predicted average σ, S, S²σ (PF) and κ_el for both p- and n-type CZTS with n = 10¹⁶–10²² cm⁻³ at a fixed T = 600 K. At this temperature, the σ is negligible up to n ≈ 10¹⁸ cm⁻³ for the n-type material and n ≈ 10¹⁹ cm⁻³ for the p-type material, after which the electrical conductivity rises substantially with n according to:

σ(n,T) = neμ(n,T) (18)

where μ is the carrier mobility and in general depends on both n and T. Whereas the p-type σ increases to the largest n = 10²² cm⁻³ in our calculations, the n-type σ reaches a maximum around 10²¹ cm⁻³ before decreasing sharply to ∼40% of its peak value at 10²² cm⁻³.

Fig. 7 Calculated electrical properties of p- and n-type Cu₂ZnSnS₄ as a function of extrinsic carrier concentration n (“doping level”) at a fixed temperature T = 600 K: (a) electrical conductivity σ, (b) absolute Seebeck coefficient |S|, (c) power factor S²σ (PF), and (d) electronic thermal conductivity κ_el. For each property we show the scalar average values obtained using the equivalents of eqn (5) in the text.

To understand this, we analysed the carrier mobilities and scattering rates from the four scattering mechanisms, and the combined (average) μ in eqn (18) (Fig. S2–S4, ESI†). For p-type doping, the μ is limited by POP scattering except at larger n → 10²² cm⁻³ where IMP scattering becomes competitive (Fig. S2, ESI†). The μ^POP is roughly constant for n ≤ 10¹⁸ cm⁻³, above which it drops sharply with n. In contrast, the μ^POP for n-type CZTS is an order of magnitude lower than for the p-type material at small n and remains lower up to around n ≥ 10²¹ cm⁻³ (Fig. S3, ESI†). However, IMP and PIE scattering become significant at lower n, and the higher scattering rates of all three mechanisms explain the unusual behaviour of the σ of n-type CZTS at large n (Fig. S4, ESI†).

For p-type doping, the S show an almost monotonic decrease with n, tending to zero at the largest n = 10²² cm⁻³ we examined. The absolute |S| with n-type doping show similar behaviour at low n, but the decrease flattens above ∼10¹⁹ cm⁻³, which we attribute to the unusual behaviour of the conductivity. The decrease in |S| with n is typical of the inverse relationship expected from:

where h is the Planck constant and is the Seebeck effective mass. The is related to the band effective mass by:⁹⁰where N_V is the valley degeneracy. A multi-valley electronic structure with large N_V is beneficial to a high |S|, and in particular can offset the impact of the increased n in cases where filling the valence or conduction bands involves additional bands in the transport (so-called “band convergence”).¹⁰ With reference to the band structure in Fig. 6, the valence-band maximum (VBM) is triply degenerate whereas the conduction-band minimum (CBM) is a single isolated band. This means that the larger |S| predicted for p-type doping is likely due at least in part to the higher N_V of the CBM. In terms of band convergence, we predict that n = 10²⁰ and n = 10²¹ cm⁻³ correspond to filling the valence band to −0.24 and −0.84 eV below the VBM, respectively. At the latter energy (larger n) there are additional bands beyond the feature around k = Γ, indicating the potential for band convergence to improve the electrical properties. On the other hand, even the largest n = 10²² cm⁻³ we tested, which corresponds to the addition of 1.56 electrons per unit cell, is insufficient to fill beyond the isolated first conduction band, so it would likely be more challenging to improve the n-type electrical properties with band convergence.

The combination of the rise in σ and decrease in |S| produces a peak in the S²σ. For p-type CZTS we obtain a relatively sharp maximum around 10²¹ cm⁻³, whereas for n-type doping we obtain a broader maximum with similar PFs from ∼10¹⁹ to 2 × 10²⁰ cm⁻³. The much smaller |S| predicted for n-type CZTS results in a smaller predicted maximum PF of 0.58 mW m⁻¹ K⁻² (σ = 1218 S cm⁻¹, S = −69 μV K⁻¹) compared to 1.53 mW m⁻¹ K⁻² for p-type doping (σ = 1275 S cm⁻¹, S = 110 μV K⁻¹). We therefore predict p-type CZTS to show significantly better performance at most n, with an approximately three times higher maximum PF.

At large n the electronic thermal conductivity κ_el closely follows the σ according to the Wiedemann–Franz relation:

κ_el(n,T) = Lσ(n,T)T (21)

where L is the Lorentz number and may not be constant.⁹¹ Increasing the n to optimise the PFs thus has the undesirable side effect of increasing the κ_el, with non-negligible values of 1.5 and 1.68 W m⁻¹ K⁻¹ predicted at the n that maximise the PFs.

The temperature dependence of the electrical transport properties in the degenerate regime at n = 5 × 10²⁰ cm⁻³, chosen to be between the n where the maximum p- and n-type PFs are predicted to occur, is shown in Fig. S5 (ESI†). At this n, the electrical conductivity shows the metallic-like decrease with temperature expected for a degenerate semiconductor. Analysis of the temperature-dependent mobilities for the four scattering mechanisms indicates that this is primarily driven by an increase in the POP scattering rates (Fig. S6 and S7, ESI†). On the other hand, the absolute Seebeck coefficients show a significant rise with T, with a predicted 4–5× increase in the averaged |S| for both p- and n-type CZTS from 300–1000 K. For p-type CZTS, this results in a peak in the PF between 700–750 K, followed by a slow decay at higher temperature, while for the n-type material we predict a monotonic rise in the PF that levels off at higher T. The electronic thermal conductivity peaks around 400 K, which is largely due to the fall in electrical conductivity with temperature outweighing the T term in eqn (21).

In contrast to the κ_latt, the σ for both p- and n-type CZTS is predicted to be strongly anisotropic (Fig. S8, ESI†), with 1.8 and 1.2× larger electrical conductivity along the a = b axes, respectively, at the n where the maximum PFs are obtained. For both doping types, there is also an almost constant offset in the |S_αβ| along the two independent directions, which becomes significant at the larger n that maximise the power factors. As a result, the maximum PFs are predicted to be 5.5× larger along the a = b direction for p-type CZTS, and 2.1× larger for n-type CZTS. These translate to modest 1.8 and 1.2× enhancement of the PFs along the a = b direction compared to the scalar averages. While this is interesting, when set against the large reductions in the κ_latt we predict can be obtained by nanostructuring this finding is unlikely to be of practical use in optimising CZTS-based thermoelectric devices.

Given the large predicted effect of crystallite size on the κ_latt, we performed a calculation of the electrical properties at a reference n = 5 × 10²⁰ cm⁻³ and T = 600 K including a boundary-scattering term analogous to eqn (10). Fig. 8 shows the predicted PF and κ_el as a function of L for p- and n-type doping. Decreasing the grain size begins to impact the electrical properties from around L = 100 nm, which is 1–2 orders of magnitude lower than the L = 1–10 μm below which the κ_latt starts to drop (cf. Fig. 5). This indicates that the electron MFPs are significantly shorter than the phonon MFPs. Decreasing the L from a bulk 1 mm to 1 nm reduces the p-type σ and κ_el by ∼50% and reduces the S by 4%. On the other hand, the n-type σ/κ_el are predicted to decrease to around 25% of the bulk values at L = 1 nm, but the |S| is predicted to increase by 50% from −44 μV K⁻¹ at L = 1 mm to −65 μV K⁻¹ at L = 1 nm. For p-type doping, reducing the L has a similar relative impact on both the PF and κ_el, with both falling by around 10, 20 and 50% at L = 10, 5 and 1 nm. For n-type doping, at the larger L = 5–10 nm the fall in the σ is partially offset by the increase in |S| so that the PFs are reduced by 3.5 and 11% of the bulk values, while at L = 1 nm we predict a similar 50% reduction to p-type CZTS. However, the κ_el is predicted to show significant decreases of ∼30, 50 and 75% at L = 10, 5 and 1 nm. Taken together with the predicted behaviour of the κ_latt, this analysis suggests a crystallite size in the range of 5–10 nm is likely to yield an optimal ZT.

Fig. 8 Dependence of the scalar average power factor S²σ (a) and electronic thermal conductivity κ_el (b) of p- and n-type Cu₂ZnSnS₄ on the crystallite size L at n = 5 × 10²⁰ cm⁻³ and T = 600 K estimated using an equivalent boundary-scattering model to eqn (10) in the text.

Table 2 compares our calculated bulk electrical properties to selected experimental measurements on pristine p-type CZTS.^12,86,92 Our calculations provide a relatively close match to the measurements in ref. 12, although reproducing the σ requires a 20% underestimate of the S, reproducing the S requires an order of magnitude smaller σ, and reproducing the PF requires a 1.8× larger σ and a 1.3× smaller S. There is also an order of magnitude variation in the n for the different predictions, although as the n was not measured we cannot determine which is closer to the experiments. For the other two measurements, we find the calculated n that reproduce the σ overestimate the S by a factor of 2–3×, whereas the n that reproduce the S overestimate the σ by 35–50×. Moreover, while our calculations can reproduce the measured PFs, they underestimate the σ and overestimate the S compared to the measurements. This may be because, while the calculations include a boundary-scattering model to limit the carrier mean-free paths, they do not take into account other effects that could arise in small crystallites such as carrier trapping at defects and energy filtering.⁹³ Similarly, as alluded to during the discussion of the lattice thermal conductivity in the previous section, it is possible that strain, defects and/or dopants could impact the electrical properties in ways that cannot be reproduced by the boundary-scattering approach. These phenomena could easily account for the order-of-magnitude difference between the σ measured in the studies in ref. 86 and 92 and the values we predict based on the measured Seebeck coefficients.

Table 2 Comparison of the predicted scalar average electrical conductivity σ, Seebeck coefficient S and power factor S²σ for Cu₂ZnSnS₄ to selected experimental measurements on pristine CZTS samples.^12,86,92 The crystallite sizes L reported in the experimental studies and used as the boundary-scattering parameter in eqn (10) in the calculations are shown in parentheses. For each measurement we provide predictions at the three n that best reproduce the σ, S and PF

	n [cm^−3]	T [K]	σ [S cm^−1]	S [μV K^−1]	S^2σ (PF) [mW m^−1 K^−2]
a Grain size estimated based on electron microscopy.b Measurements compared to predictions for bulk CZTS as crystallite size has negligible impact at large L (cf. Fig. 8).
Expt^12 (10.6 ± 1.9 nm)	—	700	0.37	990	3.6 × 10^−2
Calc. (10 nm, best sigma)	2.9 × 10^17	700	0.35	783	2.1 × 10^−2
Calc. (10 nm, best S)	2.3 × 10^16	700	0.03	990	2.8 × 10^−3
Calc. (10 nm, best PF)	5 × 10^17	700	0.69	724	3.6 × 10^−2
Expt^86 (32.7 nm)	1.5 × 10^18	573	1.58	296	1.4 × 10^−2
Calc. (33 nm, best sigma)	8.7 × 10^17	573	1.54	644	6.4 × 10^−2
Calc. (33 nm, best S)	5.7 × 10^19	573	75.7	298	0.67
Calc. (33 nm, best PF)	1.1 × 10^17	573	0.2	819	1.4 × 10^−2
Expt^92^a (2 μm)	—	663	12.4	170	3.6 × 10^−2
Calc. (bulk,^b best sigma)	9.3 × 10^18	663	11.7	482	0.27
Calc. (bulk,^b best S)	5.3 × 10^20	663	458	173	1.37
Calc. (bulk,^b best PF)	4.3 × 10^17	663	0.64	738	3.5 × 10^−2

---

Finally, it is of interest to compare the predicted electrical properties to similar calculations on orthorhombic SnS and Bi₂SO₂.^34–36 We previously predicted a maximum PF of 1.76 mW m⁻¹ K⁻² for p-type Pnma SnS at T = 880 K and n = 3.16 × 10¹⁹ cm⁻³ (σ = 226 S cm⁻¹, S = 279 μV K⁻¹),^34,35 and a maximum PF of 0.41 mW m⁻¹ K⁻² for n-type Bi₂SO₂ at a higher T = 900 K and n = 5 × 10¹⁹ cm⁻³ (σ = 120 S cm⁻¹, S = −186 μV K⁻¹).³⁶ At these temperatures, we predict that p-type CZTS can attain a maximum PF 1.58 mW m⁻¹ K⁻², which is lower than, but comparable to, p-type SnS and higher than n-type Bi₂SO₂.

3.4 Thermoelectric figure of merit

We now combine the calculated electrical transport properties and lattice thermal conductivity to predict the thermoelectric figure of merit ZT of p- and n-type CZTS as a function of both temperature and carrier concentration. Given that our analysis of the κ_latt indicated that nanostructuring should significantly enhance the ZT, and that the majority of experiments on CZTS have focussed on thin films, we discuss here the averaged ZT calculated from the scalar average S, σ, κ_el and κ_latt (cf. eqn (5)). Nanostructuring to crystallite sizes on the order of 10 nm appears to be readily achievable in experiments,^65–70 while effective crystallite sizes below 10 nm may be accessible through e.g. defect engineering. We therefore provide predictions for bulk CZTS and for crystallite sizes of 10 and 5 nm using the boundary scattering model in eqn (10) for the κ_latt and the equivalent for the S, σ and κ_el.

Fig. 9 shows the predicted ZT as a function of n and T for bulk and nanostructured p- and n-type CZTS. In general, the largest ZT are predicted at high T and large n, due to minimising the κ_latt and maximising the S²σ (cf. Fig. 6(a) and 7(c)). Both p- and n-type CZTS are predicted to show the maximum ZT, ZT_max, at the highest T = 1260 K in our calculations and at doping levels close to the n at which the peak PFs are obtained, viz. 10²¹ and 10²⁰ cm⁻³, respectively. As expected from our analysis of the electrical transport properties, p-type CZTS is predicted to perform significantly better than n-type CZTS, with double the bulk ZT_max (0.66 vs. 0.33), but neither bulk ZT_max meets the industry requirement of ZT > 1. However, we predict that significantly larger ZT_max can be obtained by nanostructuring to reduce the large bulk κ_latt. We predict that nanostructuring to 10 nm effectively doubles the p- and n-type ZT_max to 1.28 and 0.56, respectively, while further reducing the crystallite size to 5 nm will raise the ZT_max to 1.56 and 0.62.

Fig. 9 Predicted scalar average thermoelectric figure of merit ZT of (a), (c) and (e) p-type and (b), (d) and (f) n-type Cu₂ZnSnS₄ (CZTS) as a function of extrinsic carrier concentration (“doping level”) n and temperature T. Predictions are shown for (a) and (b) bulk CZTS and nanostructured CZTS with crystallite sizes of (c) and (d) L = 10 nm and (e) and (f) L = 5 nm.

To better define potential application areas, Table 3 reports the predicted ZT_max attainable for p-type CZTS at 400, 600 and 1000 K, corresponding roughly to the three heat-recovery scenarios outlined in ref. 94. The predicted ZT_max = 0.22 and 0.50 at 400 and 600 K are below the industry-standard ZT > 1, whereas at 1000 K we predict that a ZT_max ≈ 1 could be achieved by nanostructuring to 10 nm, and increased by 21% to 1.15 by reducing the crystallite size to 5 nm. Data for n-type CZTS are provided in Table S3 (ESI†), and we predict ZT_max of 0.35 and 0.38 for crystallite sizes of 10 nm and 5 nm, respectively, at T = 1000 K.

Table 3 Predicted maximum scalar average figure of merit ZT_max of p-type Cu₂ZnSnS₄ (CZTS) at T = 400, 600 and 1000 K, corresponding roughly to the three heat-recovery scenarios outlined in ref. 94. Predictions are given for bulk CZTS and with nanostructuring to crystallite sizes of L = 10 and 5 nm. For each prediction we show the corresponding extrinsic carrier concentration (“doping level”) n, electrical conductivity σ, Seebeck coefficient S, power factor S²σ (PF), and electronic, lattice and total thermal conductivity κ_el, κ_latt and κ_tot

T [K]	Crystallite size L	ZTmax	n [cm^−1]	σ [S cm^−1]	S [μV K^−1]	S^2σ (PF) [mW m^−1 K^−2]	κ [W m^−1 K^−1]
							κel	κlatt	κtot
400	Bulk	0.06	2 × 10^21	21.0	78.1	1.28	1.84	6.91	8.75
400	10 nm	0.17	2.5 × 10^20	3.98	137	0.74	0.35	1.43	1.78
400	5 nm	0.22	10^20	1.70	185	0.59	0.15	0.93	1.08
600	Bulk	0.16	10^21	8.56	130	1.46	1.08	4.53	5.61
600	10 nm	0.40	2.5 × 10^20	2.49	201	1.00	0.32	1.19	1.52
600	5 nm	0.50	2.5 × 10^20	2.29	200	0.91	0.29	0.81	1.10
1000	Bulk	0.45	7.5 × 10^20	3.73	199	1.48	0.57	2.69	3.27
1000	10 nm	0.95	5 × 10^20	2.40	228	1.24	0.38	0.93	1.31
1000	5 nm	1.15	5 × 10^20	2.26	226	1.16	0.35	0.65	1.00

---

We emphasise again that our high-temperature predictions are for pristine CZTS. The quaternary composition of CZTS can make it prone to decomposition at elevated temperature through loss of SnS and S,^95–97 resulting in the formation of off-stoichiometric and secondary phases that are not accounted for in our calculations.

Our calculations yield favourable predictions of the maximum ZT between 0.025–0.039 reported in experimental studies on pristine CZTS (Table 4).^12,28,98 Our predictions match the PFs, total thermal conductivity κ_tot and ZT reported in ref. 12 and 28 near quantitatively. On the other hand, taking the relatively large 300 nm grain size reported in ref. 98 as a crystallite size, our calculations overestimate the PF and κ_tot by a factor of ∼6.7×. As discussed in Section 3.2, this can likely be attributed to polycrystalline grains with a much smaller crystallite size. The relatively good predictions for the figure of merit of CZTS suggests that the discrepancies in the predicted electrical conductivity and Seebeck coefficients highlighted in the previous section to some extent compensate one another to yield reasonable PFs and ZT. While none of these three studies report carrier concentrations, we estimate n between 5 × 10¹⁷ and 3 × 10¹⁸ cm⁻³, and we can tentatively ascribe the low ZT to the n being much lower than we predict are required to optimise the power factor.

Table 4 Comparison of our predicted scalar average figure of merit ZT for Cu₂ZnSnS₄ to selected experimental measurements on pristine^12,28,98 and doped CZTS.^29,30 For each entry we also provide the corresponding carrier concentration n, temperature T, crystallite size L, power factor S²σ and total thermal conductivity κ_tot

	n [cm^−3]	T [K]	Crystallite/grain size L	S^2σ (PF) [mW m^−1 K^−2]	κtot [W m^−1 K^−1]	ZT
a Electrical properties calculated based on bulk CZTS as crystallite size has negligible impact at large L (cf. Fig. 8).
Expt (Cu2ZnSnS4)^28	—	700	—	6.7 × 10^−2	1.21	3.9 × 10^−2
Calc.	1.15 × 10^18	700	10 nm	5.9 × 10^−2	1.11	3.7 × 10^−2
Expt (Cu2ZnSnS4)^12	—	700	10.6 ± 1.9 nm	3.6 × 10^−2	0.97	2.5 × 10^−2
Calc.	6.6 × 10^17	700	10 nm	3.9 × 10^−2	1.11	2.5 × 10^−2
Expt (Cu2ZnSnS4)^98	—	623	250–300 nm	2.1 × 10^−2	0.53	2.5 × 10^−2
Calc.	3.1 × 10^18	623	300 nm + bulk^a	0.14	3.66	2.4 × 10^−2
Expt (Na-doped Cu2ZnSnS4)^29	—	800	—	1.86	0.95	1.6
Calc.	2.3 × 10^20	800	1 nm	0.60	0.39	1.22
Calc.	3 × 10^20	800	3 nm	0.93	0.79	0.95
Expt (Cu2(Zn,Fe,Cd)SnS4)^30	—	550	1.3 μm	0.85	0.68	0.69
Calc.	1.2 × 10^21	550	1 μm + bulk^a	1.44	5.91	0.13
Calc.	2 × 10^20	550	3 nm	0.75	0.82	0.49

---

On the other hand, comparison to the Na-doped CZTS and Cu₂(Zn,Fe,Cd)SnS₄ solid solutions in ref. 29 and 30 is less favourable. Reproducing the ZT = 1.6 at 800 K measured in ref. 29 requires a large n = 2.3 × 10²⁰ cm⁻³ and a very small crystallite size of L = 1 nm, which results in a ∼59–68% underestimate of the measured PF and κ_tot. Though not typical, L in the range of 2.3–7.7 nm have been reported in a number of experiments.^99–103 With L = 3 nm the calculations predict a ZT of 0.95, albeit again with an underestimated PF and κ_tot. Similarly, to obtain a comparable figure of merit at 550 K to the ZT = 0.69 reported in ref. 30 requires much smaller crystallites than the reported grain size of 1.3 μm, and the calculations still underestimate the PF and overestimate the κ_tot by 11.8% and 20% respectively.

Taking these two sets of comparisons together, we are reasonably confident in our predictions for the figure of merit of pristine CZTS with modest n. However, our results suggest more in-depth modelling is needed to ascertain the effect of heavy doping to obtain larger ZT on the electronic structure and transport properties. For example, the underpredicted PF of the Cu₂(Zn,Fe,Cd)SnS₄ solid solutions and Na-doped CZTS could be due to the formation of resonant levels that are not captured with the rigid-band approximations.

Finally, it is also of interest to compare the predicted ZT to our previous calculations on SnS and Bi₂SO₂.^34–36 We predicted p-type Pnma SnS to reach a maximum ZT of 1.39 at 860 K (n = 3.2 × 10¹⁹ cm⁻³, PF = 1.76 mW m⁻¹ K⁻², κ_tot = 1.09 W m⁻¹ K⁻¹).³⁵ We predict that bulk p-type CZTS would have a much lower maximum ZT of 0.34 at this temperature, but that limiting the crystallite sizes to 10 and 5 nm would increase this to 0.76 and 0.92. In the latter case, we predict a PF of 1.11 mW m⁻¹ K⁻² and a κ_tot of 1.04 W m⁻¹ K⁻¹ (n = 4.04 × 10²⁰ cm⁻³). Compared to SnS, CZTS has a lower power factor and a much larger intrinsic κ_latt, which requires both a higher doping level and nanostructuring to offset in order to obtain a comparable ZT. On the other hand, we predicted n-type Bi₂SO₂ to reach a low ZT of 0.33 at 900 K (2.5 × 10¹⁹ cm⁻³, 0.41 mW m⁻¹ K⁻², 1.13 W m⁻¹ K⁻¹),³⁶ whereas we predict n-type CZTS to show a maximum ZT of 0.27 with a crystallite size to 10 nm and n = 1.87 × 10¹⁹ cm⁻³ (0.38 mW m⁻¹ K⁻², 1.25 W m⁻¹ K⁻¹). Nanostructured CZTS could therefore potentially be competitive with Bi₂SO₂ as an n-type sulphide TE. We note, however, that we predicted SnS and Bi₂SO₂ to have potential ZT of 1.78 and 2.53, respectively, with the less widely studied n- and p-type doping, which are both far higher than we predict is achievable with CZTS.

4 Conclusions

In summary, we have used ab initio modelling to explore the thermoelectric performance of the quaternary sulphide Cu₂ZnSnS₄. Analysis of the lattice thermal conductivity suggests that CZTS has unusually long phonon lifetimes due to weak anharmonicity and a small scattering phase space. This results in the κ_latt being highly sensitive to the crystallite size such that the bulk thermal conductivity can be reduced by >80% at the crystallite sizes of L < 10 nm reported in experiments. The calculated electrical properties indicate that p-type CZTS can support a considerably larger maximum power factor than can be obtained with electron doping, due to a higher Seebeck coefficient and a larger electrical conductivity at the n ≈ 10²¹ cm⁻³ required to optimise the p-type PF. The electrical properties are much less sensitive to the crystallite size than the κ_latt, with the PFs unaffected above L ≃ 100 nm. We predict that nanostructured p-type CZTS with a crystallite size of 10 nm should yield a high-temperature ZT > 1, with the potential for improvement through further reduction in the (effective) crystallite size, although similar calculations predict that larger maximum ZT can potentially be obtained with other (oxy-)sulphide TEs. The lower maximum ZT compared to p-type SnS is partly due to the lower power factor, and so exploring the selenide analogue Cu₂ZnSnSe₄, or other ternary or quaternary chalcogenides with the diamondoid structure, may prove fruitful. In addition, comparison of the predicted ZT to experiments suggests that heavy doping or alloying may have the dual function of reducing the κ_latt and improving the electrical transport properties beyond optimising the carrier concentration for a rigid band structure, and further computational studies to better understand these effects would also, in our view, be highly informative.

Author contributions

Conceptualisation – JMS; methodology – all authors; data curation – TS and JMF; formal analysis and investigation – all authors; writing (original draft) – TS; writing (review and editing) – all authors; funding acquisition – JMS. All authors have read and agreed to the published version of the paper.

Conflicts of interest

There are no conflicts to declare.

Data availability

Raw data from this study will be made available to download free of charge after publication from an online repository at https://doi.org/10.17632/4d699zwynb. Our analysis code is available on GitHub at https://github.com/skelton-group/ZT-Calc-Workflow.

Acknowledgements

JMS and JMF are supported by a UK Research and Innovation Future Leaders Fellowship (MR/T043121/1, MR/Y033973/1), and JMS previously held a University of Manchester Presidential Fellowship. This work used the ARCHER2 UK National Supercomputing Service via the UK Materials Chemistry Consortium (MCC; UK Engineering and Physical Sciences Research Council EP/R029431 and EP/X035859). We also made use of the UoM Computational Shared Facility (CSF), which is maintained by UoM Research IT.

Notes and references

1. A. Firth, B. Zhang and A. Yang, Appl. Energy, 2019, 235, 1314–1334.
2. R. Freer and A. V. Powell, J. Mater. Chem. C, 2020, 8, 441–463.
3. H. J. Goldsmid, The Seebeck and Peltier effects, Morgan & Claypool Publishers, New York, 2017, ch. 1, pp. 1-1–1-3.
4. C. Han, Q. Sun, Z. Li and S. X. Dou, Adv. Energy Mater., 2016, 6, 1600498.
5. L. E. Bell, Science, 2008, 321, 1457–1461.
6. J. He and T. M. Tritt, Science, 2017, 357, eaak9997.
7. H. Jouhara, A. Żabnieńska Góra, N. Khordehgah, Q. Doraghi, L. Ahmad, L. Norman, B. Axcell, L. Wrobel and S. Dai, Int. J. Thermofluids, 2021, 9, 100063.
8. C. Gayner and K. K. Kar, Prog. Mater. Sci., 2016, 83, 330–382.
9. Y. Liu, J. Zhi, W. Li, Q. Yang, L. Zhang and Y. Zhang, Molecules, 2023, 28, 5894.
10. G. Tan, L.-D. Zhao and M. G. Kanatzidis, Chem. Rev., 2016, 116, 12123–12149.
11. A. Kumar, S. Bano, B. Govind, A. Bhardwaj, K. Bhatt and D. K. Misra, J. Electron. Mater., 2021, 50, 6037–6059.
12. H. Yang, L. A. Jauregui, G. Zhang, Y. P. Chen and Y. Wu, Nano Lett., 2012, 12, 540–545.
13. K. Biswas, J. He, I. D. Blum, C.-I. Wu, T. P. Hogan, D. N. Seidman, V. P. Dravid and M. G. Kanatzidis, Nature, 2012, 489, 414–418.
14. Y. Lan, A. J. Minnich, G. Chen and Z. Ren, Adv. Funct. Mater., 2010, 20, 357–376.
15. T. Hendricks and W. T. Choate, Engineering Scoping Study of Thermoelectric Generator Systems for Industrial Waste Heat Recovery, 2006.
16. H. J. Goldsmid, Materials, 2014, 7, 2577–2592.
17. R. Freer, D. Ekren, T. Ghosh, K. Biswas, P. Qiu, S. Wan, L. Chen, S. Han, C. Fu, T. Zhu, A. K. M. A. Shawon, A. Zevalkink, K. Imasato, G. J. Snyder, M. Ozen, K. Saglik, U. Aydemir, R. Cardoso-Gil, E. Svanidze, R. Funahashi, A. V. Powell, S. Mukherjee, S. Tippireddy, P. Vaqueiro, F. Gascoin, T. Kyratsi, P. Sauerschnig and T. Mori, J. Phys.: Energy, 2022, 4, 022002.
18. L.-D. Zhao, S.-H. Lo, Y. Zhang, H. Sun, G. Tan, C. Uher, C. Wolverton, V. P. Dravid and M. G. Kanatzidis, Nature, 2014, 508, 373–377.
19. L.-D. Zhao, G. Tan, S. Hao, J. He, Y. Pei, H. Chi, H. Wang, S. Gong, H. Xu, V. P. Dravid, C. Uher, G. J. Snyder, C. Wolverton and M. G. Kanatzidis, Science, 2016, 351, 141–144.
20. C. Zhou, Y. K. Lee, Y. Yu, S. Byun, Z.-Z. Luo, H. Lee, B. Ge, Y.-L. Lee, X. Chen, J. Y. Lee, O. Cojocaru-Mirédin, H. Chang, J. Im, S.-P. Cho, M. Wuttig, V. P. Dravid, M. G. Kanatzidis and I. Chung, Nat. Mater., 2021, 20, 1378–1384.
21. M. Ibáñez, R. Zamani, A. LaLonde, D. Cadavid, W. Li, A. Shavel, J. Arbiol, J. R. Morante, S. Gorsse, G. J. Snyder and A. Cabot, J. Am. Chem. Soc., 2012, 134, 4060–4063.
22. W. G. Zeier, Y. Pei, G. Pomrehn, T. Day, N. Heinz, C. P. Heinrich, G. J. Snyder and T. Wolfgang, J. Am. Chem. Soc., 2013, 135, 726–732.
23. X. Zheng, Y. Liu, Y. Du, Y. Sun, J. Li, R. Zhang, Q. Li, P. Chen, G. Zhao, Y. Fang and N. Dai, J. Alloys Compd., 2018, 738, 484–490.
24. G. Jamwal, M. Warish, S. Muthiah, S. Chakravarty, N. Jakhar, A. Kandasami and A. Niazi, Inorg. Chem., 2022, 61, 16390–16404.
25. K. Pal, P. Singh, A. Bhaduri and K. B. Thapa, Sol. Energy Mater. Sol. Cells, 2019, 196, 138–156.
26. U. A. Shah, A. Wang, M. Irfan Ullah, M. Ishaq, I. A. Shah, Y. Zeng, M. S. Abbasi, M. A. Umair, U. Farooq, G.-X. Liang and K. Sun, Small, 2024, 20, 2310584.
27. K. Gupta, S. Gupta and Y. Batra, Mater. Sci. Eng., B, 2024, 303, 117291.
28. M.-L. Liu, I.-W. Chen, F.-Q. Huang and L.-D. Chen, Adv. Mater., 2009, 21, 3808–3812.
29. A. Nagaoka, K. Yoshino, T. Masuda, T. D. Sparks, M. A. Scarpulla and N. Kensuke, J. Mater. Chem. A, 2021, 9, 15595–15604.
30. Y. Liu, P. D. McNaughter, X. Liu, A. V. Kretinin, J. M. Skelton, F. Azough, D. J. Lewis and R. Freer, ACS Appl. Mater. Interfaces, 2024, 16, 11516–11527.
31. Z. Lin, C. Hollar, J. S. Kang, A. Yin, Y. Wang, H.-Y. Shiu, Y. Huang, Y. Hu, Y. Zhang and X. Duan, Adv. Mater., 2017, 29, 1606662.
32. C. Yang, D. Souchay, M. Kneiß, M. Bogner, H. M. Wei, M. Lorenz, O. Oeckler, G. Benstetter, Y. Q. Fu and M. Grundmann, Nat. Commun., 2017, 8, 16076.
33. R. Venkatasubramanian, E. Siivola, T. Colpitts and B. O’Quinn, Nature, 2001, 413, 597–602.
34. J. M. Flitcroft, I. Pallikara and J. M. Skelton, Solids, 2022, 3, 155–176.
35. M. Zhang, J. M. Flitcroft, S. K. Guillemot and J. M. Skelton, J. Mater. Chem. C, 2023, 11, 14833–14847.
36. J. M. Flitcroft, A. Althubiani and J. M. Skelton, J. Phys.: Energy, 2024, 6, 025011.
37. G. Kresse and J. Hafner, Phys. Rev. B:Condens. Matter Mater. Phys., 1993, 47, 558.
38. J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou and K. Burke, Phys. Rev. Lett., 2008, 100, 136406.
39. P. E. Blöchl, Phys. Rev. B:Condens. Matter Mater. Phys., 1994, 50, 17953.
40. G. Kresse and D. Joubert, Phys. Rev. B:Condens. Matter Mater. Phys., 1999, 59, 1758.
41. H. J. Monkhorst and J. D. Pack, Phys. Rev. B, 1976, 13, 5188.
42. A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder and K. A. Persson, APL Mater., 2013, 1, 011002.
43. A. Togo, K. Shinohara and I. Tanaka, Spglib: a software library for crystal symmetry search, arXiv, 2024, preprint, arXiv:1808.01590 DOI:10.48550/arXiv.1808.01590.
44. A. Togo and I. Tanaka, Scr. Mater., 2015, 108, 1–5.
45. A. Togo, L. Chaput and I. Tanaka, Phys. Rev. B:Condens. Matter Mater. Phys., 2015, 91, 094306.
46. A. M. Ganose, J. Park, A. Faghaninia, R. Woods-Robinson, K. A. Persson and A. Jain, Nat. Commun., 2021, 12, 2222.
47. A. V. Krukau, O. A. Vydrov, A. F. Izmaylov and G. E. Scuseria, J. Chem. Phys., 2006, 125, 224106.
48. M. Gajdoš, K. Hummer, G. Kresse, J. Furthmüller and F. Bechstedt, Phys. Rev. B:Condens. Matter Mater. Phys., 2006, 73, 045112.
49. K. Momma and F. Izumi, J. Appl. Crystallogr., 2011, 44, 1272–1276.
50. V. Kheraj, K. K. Patel, S. J. Patel and D. V. Shah, J. Cryst. Grow., 2013, 362, 174–177.
51. T. Chandel, V. Thakur, S. Halaszova, M. Prochazka, D. Haško, D. Velic and R. Poolla, J. Electron. Mater., 2018, 47, 5477–5487.
52. J. Guo, S. Sun, B. Liu, R. Hao and L. Sun, Optik, 2021, 242, 166998.
53. M. Sahu, V. R. M. Reddy, B. Kim, B. Patro, C. Park, W. K. Kim and P. Sharma, Materials, 2022, 15, 1708.
54. N. B. M. Amiri and A. Postnikov, Phys. Rev. B:Condens. Matter Mater. Phys., 2010, 82, 205204.
55. A. J. Jackson and A. Walsh, J. Mater. Chem. A, 2014, 2, 7829–7836.
56. S. P. Ramkumar, G. Petretto, W. Chen, H. P. C. Miranda, X. Gonze and G.-M. Rignanese, Phys. Rev. Mater., 2022, 6, 035403.
57. L. Boutahar, A. Benamrani, Z. Rouabah and S. Daoud, Ann. West Univ. Timisoara, Phys., 2023, 65, 160–170.
58. J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys., 2003, 118, 8207–8215.
59. J. M. Skelton, A. J. Jackson, M. Dimitrievska, S. K. Wallace and A. Walsh, APL Mater., 2015, 3, 041102.
60. F. A. Cotton, G. Wilkinson, C. A. Murillio and M. Bochmann, Advanced Inorganic Chemistry, Wiley, 6th edn, 1999.
61. J. M. Skelton, J. Mater. Chem. C, 2021, 9, 11772–11787.
62. J. Cen, I. Pallikara and J. M. Skelton, Chem. Mater., 2021, 33, 8404–8417.
63. J. Tang and J. M. Skelton, J. Phys.: Condens. Matter, 2021, 33, 164002.
64. S. K. Guillemot, A. Suwardi, N. Kaltsoyannis and J. M. Skelton, J. Mater. Chem. A, 2024, 12, 2932–2948.
65. S. A. Vanalakar, G. L. Agawane, S. W. Shin, M. P. Suryawanshi, K. V. Gurav, K. S. Jeon, P. S. Patil, C. W. Jeong, J. Y. Kim and J. H. Kim, J. Alloys Compd., 2015, 619, 109–121.
66. N. Kamoun, H. Bouzouita and B. Rezig, Thin Solid Films, 2007, 515, 5949–5952.
67. M. Abusnina, M. Matin, H. R. Moutinho, J. L. Blackburn, J. Alleman, C. DeHart, B. To and M. Al-Jassim, IEEE J. Photovolt., 2015, 5, 1470–1475.
68. D. M. M. Alam, Chalcogenide Lett., 2014, 11, 233–239.
69. M. Byeon, J.-S. Bae, T.-E. Hong, E.-D. Jeong, S. Kim and Y. Kim, Korean J. Mater. Res., 2013, 23, 7–12.
70. J. Zhang, B. Long, S. Cheng and W. Zhang, Int. J. Photoenergy, 2013, 2013, 1–6.
71. W. D. Thompson, A. Nandur and B. E. White, J. Appl. Phys., 2016, 119, 095108.
72. S. D. Sharma, B. Khasimsaheb, Y. Y. Chen and S. Neeleshwar, Ceram. Int., 2019, 45, 2060–2068.
73. N. Kattan, B. Hou, D. J. Fermín and D. Cherns, Appl. Mater. Today, 2015, 1, 52–59.
74. R. Ahmad, M. Brandl, M. Distaso, P. Herre, E. Spiecker, R. Hock and W. Peukert, CrystEngComm, 2015, 17, 6972–6984.
75. N. Song, M. Young, F. Liu, P. Erslev, S. Wilson, S. P. Harvey, G. Teeter, Y. Huang, X. Hao and M. A. Green, Appl. Phys. Lett., 2015, 106, 252102.
76. S. M. Pawar, B. S. Pawar, A. V. Moholkar, D. S. Choi, J. H. Yun, J. H. Moon, S. S. Kolekar and J. H. Kim, Electrochim. Acta, 2010, 55, 4057–4061.
77. H. Katagiri, K. Saitoh, T. Washio, H. Shinohara, T. Kurumadani and S. Miyajima, Sol. Energy Mater. Sol. Cells, 2001, 65, 141–148.
78. T. Chandel, V. Thakur, S. Halaszova, M. Prochazka, D. Haško, D. Velic and R. Poolla, J. Electron. Mater., 2018, 47, 5477–5487.
79. M. Altosaar, J. Raudoja, K. Timmo, M. Danilson, M. Grossberg, J. Krustok and E. Mellikov, Phys. Status Solidi A, 2008, 205, 167–170.
80. S. I. Swati, R. Matin, S. Bashar and Z. H. Mahmood, J. Phys.: Conf. Ser., 2018, 1086, 012010.
81. J. Paier, R. Asahi, A. Nagoya and G. Kresse, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 115126.
82. C. Persson, R. Chen, H. Zhao, M. Kumar and D. Huang, Electronic Structure and Optical Properties from First-Principles Modeling, Wiley, 2014, pp. 75–105.
83. S. Chen, X. G. Gong, A. Walsh and S.-H. Wei, Appl. Phys. Lett., 2009, 94, 041903.
84. C.-J. Tong, H. J. Edwards, T. D. C. Hobson, K. Durose, V. R. Dhanak, J. D. Major and K. P. McKenna, J. Phys. Chem. Lett., 2020, 11, 10463–10468.
85. M. Nishiwaki, K. Nagaya, M. Kato, S. Fujimoto, H. Tampo, T. Miyadera, M. Chikamatsu, H. Shibata and H. Fujiwara, Phys. Rev. Mater., 2018, 2, 85404.
86. K. V. Gunavathy, K. Tamilarasan, C. Rangasami and A. M. S. Arulanantham, Ceram. Int., 2020, 46, 28342–28354.
87. P. R. Ghediya, T. K. Chaudhuri and D. Vankhade, J. Alloys Compd., 2016, 685, 498–506.
88. S. Bouzida, M. Battas, E. B. Benamar, G. Schmerber, A. Dinia, M. Abd-Lefdil and M. Regragui, Mater. Res. Innovations, 2022, 26, 127–133.
89. A. Z. Khan, J. M. Flitcroft and J. M. Skelton, Mater. Adv., 2024, 5, 652–664.
90. Z. M. Gibbs, F. Ricci, G. Li, H. Zhu, K. Persson, G. Ceder, G. Hautier, A. Jain and G. J. Snyder, npj Comput. Mater., 2017, 3, 8.
91. G. S. Kumar, G. Prasad and R. O. Pohl, J. Mater. Sci., 1993, 28, 4261–4272.
92. B. D. Long, N. V. Khanh, D. N. Binh and N. H. Hai, Powder Metall., 2020, 63, 220–226.
93. D. Narducci, E. Selezneva, G. Cerofolini, S. Frabboni and G. Ottaviani, J. Solid State Chem., 2012, 193, 19–25.
94. S. LeBlanc, S. K. Yee, M. L. Scullin, C. Dames and K. E. Goodson, Renewable Sustainable Energy Rev., 2014, 32, 313–327.
95. A. Weber, R. Mainz and H. W. Schock, J. Appl. Phys., 2010, 107, 013516.
96. A. Redinger, D. M. Berg, P. J. Dale and S. Siebentritt, J. Am. Chem. Soc., 2011, 133, 3320–3323.
97. A. J. Jackson and A. Walsh, J. Mater. Chem. A, 2014, 2, 7829–7836.
98. S. D. Sharma and S. Neeleshwar, MRS Adv., 2018, 3, 1373–1378.
99. Q. Wen, Y. Li, J. Yan and C. Wang, Mater. Lett., 2015, 140, 16–19.
100. M. Boshta, S. Binetti, A. Le Donne, M. Gomaa and M. Acciarri, Mater. Technol., 2017, 32, 251–255.
101. L. Chen and C. Park, Korean J. Chem. Eng., 2017, 34, 1187–1191.
102. S. Rondiya, A. Rokade, A. Jadhavar, S. Nair, M. Chaudhari, R. Kulkarni, A. Mayabadi, A. Funde, H. Pathan and S. Jadkar, Mater. Renewable Sustainable Energy, 2017, 6, 8.
103. S. Y. Gezgin, A. Houimi and H. S. Kilic, Optik, 2019, 199, 163370.

---

Footnote

† Electronic supplementary information (ESI) available: Includes comparison of the calculated lattice parameters and electronic bandgap to previous experiments and computational studies, additional data for the predicted electrical properties, and predicted maximum figures of merit and associated material parameters for n-type Cu₂ZnSnS₄. See DOI: https://doi.org/10.1039/d5tc01537e

---