Suppression of lattice thermal conductivity by mass-conserving cation mutation in multi-component semiconductors

Direct heat to electricity conversion in a thermoelectric device represents an important process for energy generation and efficiency.¹ The physical principles and applications of thermoelectric power conversion are well established; however, current commercial devices are based on heavy metal tellurides, which have issues associated with cost, toxicity, and element availability. A new generation of materials is required to support growth of this technology, ideally with performance superior to existing compounds.² 

There have already been decades of research focused on discovering materials with enhanced performance. A set of standard guidelines have emerged in the field and are widely adopted. These include high atomic weight (which led to PbTe and Bi₂Te₃),³ loosely bound “rattling” atoms (e.g., clathrates),⁴ low dimensionality (e.g., nanorods and quantum dots),⁵ and phonon localization (e.g., superlattices and alloys).⁶ Since hundreds of materials have been systematically screened, tried, and tested, new ideas are required to move this field forward.^7–9 

A critical factor underpinning thermoelectric performance is thermal conductivity, which represents a thermal loss mechanism that should be minimised for high-performance. All crystals have phonons (extended lattice vibrations), but the thermal conductivity is linked to their lifetime, which can vary by several orders of magnitude between materials.¹⁰ Both electrons and phonons can play a role in thermal conductivity, but for dielectrics and intrinsic semiconductors the latter are usually dominant. An ideal thermoelectric material will offer “phonon-glass electron-crystal” behaviour:¹¹ a crystalline system where phonons have short lifetimes (low thermal conductivity), while electron carriers have long lifetimes (high electrical conductivity).

In this study we attempt to answer a simple question: can the lattice thermal conductivity of tetrahedral semiconductors be tuned without changing the local structure or average mass of the compound? The topic was inspired by the early work of Pamplin¹² and Goodman¹³ on cation cross-substitution to produce multi-component semiconductors, which was later applied to screening materials for photovoltaics,¹⁴ topological conductivity,¹⁵ and spintronics.¹⁶ We consider the mutation from ZnS → CuGaS₂ → Cu₂ZnGeS₄ where the average atomic mass and crystal structures are conserved. A systematic analysis, from direct computation of the phonon energies and lifetimes, demonstrates a large decrease in phonon lifetimes on transition from the binary to ternary compounds, but no further enhancement is found in the quaternary system. Analysis of the results suggests avenues for lowering the thermal conductivity of semiconducting crystals.

Lattice thermal conductivity from first-principles. Within the harmonic approximation, phonons propagate indefinitely and the lattice thermal conductivity is formally infinite. Even in the absence of structural defects and chemical impurities, the interaction between phonon modes leads to scattering processes, of which anharmonic three-phonon interactions are considered to be dominant. As stated by Ziman in his seminal monograph: “Knowledge of the magnitude of the anharmonic terms which generate the [phonon-phonon scattering] processes is scanty, and can only be obtained by roundabout arguments from other general macroscopic phenomena.”¹⁷ Fifty years later, the study of many-phonon processes remains a daunting task.

By solving the phonon Boltzmann equation within the relaxation-time approximation (RTA), the lattice thermal conductivity tensor (κ) can be expressed succinctly as a sum over phonons of band index λ and wavevector q,

where C_V is the isochoric modal heat capacity, ν is the group velocity, and ν_qλτ_qλ = Λ_qλ, the phonon mean free path. The single-mode relaxation time τ can be deduced from the imaginary part of the phonon self-energy, computed within many-body perturbation theory.

At high temperatures, C_V approaches a constant, so the two critical parameters for lattice thermal conductivity are ν and Λ. The velocity ν is determined simply by the phonon dispersion with respect to reciprocal-space wavevector (⁠$∂ ω ∂ q$⁠). Λ depends on ν, but also on the relaxation time, which is closely related to the phonon-phonon interaction strength and the distribution of frequencies in the phonon density of states, the latter determining the number of possible energy-conserving scattering events.

It has only recently become possible to directly compute many-phonon interactions from first-principles.^18–20 Even though a robust infrastructure now exists, there is a large computational cost to performing the simulations for complex crystals. By calculating the three-phonon interaction strength (ϕ_λλ′λ″), τ can be computed up to second-order within many-body perturbation theory. We employ the Phonopy and Phono3py packages^20,21 using VASP,²² as the force calculator and the PBEsol²³ exchange-correlation treatment within Kohn-Sham density functional theory (DFT). PBEsol was chosen as it has been found to be an excellent compromise between accuracy and computational cost for lattice dynamics.²⁴ For each material we consider the harmonic phonon dispersion, thermal expansion within the quasi-harmonic approximation (QHA), and finally lattice thermal conductivity within the RTA. The technical setup and procedure is provided in the computational details section. It should be noted that an alternative approach to capture anharmonic interactions is through molecular dynamics simulations where, for example, the Green-Kubo method can be used to probe phonon lifetimes and thermal conductivity.²⁵ 

Zincblende ZnS. ZnS can crystallise in two polymorphs, but here we consider only zincblende (sphalerite), the most common face-centered-cubic form with space group type $F 4 ̄ 3m$ (T_d symmetry). The primitive cell contains two lattice sites, and hence there are 6 phonon modes. At the zone-centre (q = 0), the T₂ (IR and Raman active) optic branch is split into lower-energy transverse (TO) and higher-energy longitudinal (LO) modes.

In the first reported Raman measurement of ZnS crystals, the room temperature TO and LO modes were measured at 276 and 351 cm⁻¹, respectively.²⁷ The calculated values from the harmonic phonon dispersion (Figure 1) are 275 and 335 cm⁻¹, in good agreement. In addition, the phonon dispersion across the Brillouin zone agrees very well with neutron scattering measurements.²⁶ The thermal expansion, calculated within the QHA, is plotted in Figure 2. The known negative thermal expansion at low temperatures is reproduced and the room temperature value of 2 × 10⁻⁶ K is in good agreement with both experiment and reported lattice-dynamics calculations.²⁸ 

From binary to ternary: CuGaS₂. The crystal structure of CuGaS₂ is a simple 1a × 1a × 2a supercell expansion of sphalerite with the Cu and Ga atoms aligned along the (201) planes. The chalcopyrite mineral structure with space group type $I 4 ̄ 2d$ (D_2d symmetry) has eight atoms in the primitive cell and 24 phonon modes. Even at the Brillouin zone centre, the phonon structure is complicated with A₁ + 2A₂ + 3B₁ + 4B₂ + 7E vibrational branches, of which only the A₂ modes are IR and Raman inactive.

The harmonic phonon dispersion (Figure 1) is in good agreement with previous calculations.^31,32 The measured upper optic branch runs from 261 to 387 cm⁻¹ (collated in Ref. 31), which matches well with our calculated values of 255 to 377 cm⁻¹. The thermal expansion behaviour is found to be similar to ZnS, including the low-temperature lattice contraction. The linear expansion coefficient α = 25.4 × 10⁻⁶ K⁻¹ at 300 K is in agreement with the experimental measurements of 23 – 27 × 10⁻⁶ K⁻¹.^33,34

From ternary to quaternary: Cu₂ZnGeS₄. The kesterite crystal structure adopted by Cu₂ZnGeS₄ is closely related to chalcopyrite, with the addition of alternating Cu–Zn and Cu–Ge (001) planes. The space group type is $I 4 ̄$ (S₄ symmetry) and there are again eight atoms in the primitive cell, giving rise to 24 phonon modes. The zone-centre irreducible representations are 3A + 7B + 7E. All modes Raman active, while only the B and E modes are IR active.

The phonon dispersion is shown in Figure 1 (Multimedia view), and the optic branches run from 78 to 373 cm⁻¹. The measured Raman modes run from 72 to 409 cm⁻¹, which were also found to be in good agreement with prior DFT calculations using the PBE functional.³⁵ 

Anharmonicity and thermal conductivity. To assess the modal contributions to the thermal transport, the room-temperature (300 K) modal thermal-conductivity tensors κ_λ were integrated over the phonon Brillouin zone using the linear tetrahedron method. The frequency derivatives of the resulting accumulation functions were then computed, yielding a function analogous to a density of states (DoS) that quantifies the contribution of modes in different frequency bands to the overall transport. A similar procedure was followed for integrating the tensor products ν_qλ ⊗ ν_qλ, the phonon lifetimes τ_λ, and the averaged phonon-phonon interaction strengths P_qj,λ (see Ref. 20 for further details). With reference to the phonon DoS, these functions allow the modal contributions to the thermal conductivities to be interpreted in terms of mode lifetimes and group velocities (cf. Eq. (1)).

The four sets of calculated functions are overlaid on the phonon DoS curves in Figure 3. In all three materials, the majority of the heat transport is through modes up to ∼3 THz, which is due to a combination of long lifetimes and high group velocities. All three compounds also show a smaller secondary contribution to the thermal conductivity from modes between 3 and 4 THz with a comparable group velocity, but relatively shorter lifetimes. We note that the modal heat capacities are saturated at 300 K, and the contribution of this term to κ essentially mirrors the phonon DoS.

The data in Figure 3 show that the reduction in thermal conductivity on going from ZnS to CGS and CZGS is attributable both to a reduction in the mode lifetimes and a reduction in the group velocity, the latter of which can be attributed to the weaker bonding. For the most part, long-lived modes are associated with regions of the DoS in which the interaction strength is small. It is worth noting that the interaction strength is considerably larger in ZnS than in CGS and CZGS, despite its longer phonon lifetimes and high thermal conductivity. This can be explained in terms of the conservation of energy: the more complex structures of CGS and CZGS lead to a broadening of the phonon DoS, particularly the high-frequency optic branches, which in turn provides more energy-conserving scattering channels that outweigh the weaker phonon-phonon interactions.

Figure 4 shows the calculated thermal-conductivity curves for the three materials. The measured thermal conductivities for ZnS range from 360 Wm⁻¹ K⁻¹ at 30 K to 27 Wm⁻¹ K⁻¹ at 300 K,³⁶ the latter of which is within a factor of two of our calculated value. For CuGaS₂ and Cu₂ZnGeS₄, the isotropic average of the thermal conductivity $κ= 1 3 ( κ x x + κ y y + κ z z )$ is presented in Figure 4. The anisotropy in the thermal conductivity of these tetragonal crystals is small, but non-negligible. A higher conductivity is found in the ab plane (⁠$κ xx CGS$ = $κ xx CZGS =10.2$ Wm⁻¹ K⁻¹ at T = 300 K) than along the c axis (⁠$κ zz CGS =8.8$ Wm⁻¹ K⁻¹ ; $κ xx CZGS =8.0$ Wm⁻¹ K⁻¹ at T = 300 K).

Although average mass is conserved for the three compounds considered here, the mass variance caused by the 2Zn → Cu + Ga mutation may result in additional scattering analogous to the anharmonic scattering found in disordered alloys. We have tested this hypothesis by modifying ZnS with heavy and light isotopes distributed in the same arrangement as the ternary and quaternary structures (labelled “mv” in Figure 4). The 300 K thermal conductivity is reduced by 18 % for the natural isotope variation, 38% for CuGaS₂-type variation, and 46% for Cu₂ZnGeS₄-type variation. This mass variance contribution is significant, but a much weaker effect than the changes in the force constants caused by the chemical substitutions. In the actual calculations of CuGaS₂ and Cu₂ZnGeS₄, the thermal conductivity is reduced by 78% in comparison to the value for ZnS at 300 K.

As noted above, acoustic phonon modes are responsible for conducting most of the heat, owing to their large group velocity and longer lifetime. While optical phonons (less disperse bands with lower ν) do not directly contribute to thermal transport, they are important indirectly in determining mode lifetimes by their involvement in scattering processes. A clear change between the phonon DoSs on transition from the binary to multernary materials is that the optic branch of the DoS is widened. These modes mostly involve motion of the anion sub-lattice — animations of the modes are provided as a multimedia view in Figure 1 — and the distribution of environments found when multiple cations are introduced causes a spread in frequency. The result is a higher probability of three-phonon interactions, limiting the lifetimes, and this, together with a reduction in the group velocity, serves to reduce the thermal conductivity.

In the 1970 overview by Spitzer¹⁰ on the thermal conductivity of semiconductors, an important observation was made: “It is found that lattice thermal conductivity may be correlated rather reliably with crystal structure. In general, increasing coordination of the ions is associated with decreasing thermal conductivity.” We have further shown that even for a fixed coordination number, the composition can have a large impact. Now that the underlying contributions to thermal transport can be separated and quantified, deeper insights can be obtained into the structure-property relationships that can be used towards the rational engineering of thermal transport.

Computational details: For the density functional theory calculations, a kinetic energy cutoff of 450 eV for the plane wave-basis set was combined with a reciprocal-space sampling equivalent to 8 × 8 × 8 k-points for the zincblende primitive cell (i.e., 4 × 4 × 4 for CuGaS₂ and Cu₂ZnGeS₄, respectively). Projector augmented-wave (PAW) datasets³⁷ were employed. A tolerance of 10⁻⁸ eV was applied during the electronic minimization, and geometry optimisations were carried out to a force threshold of 10⁻² eV/Å. The precision of the charge-density grids was automatically chosen to avoid aliasing errors, and an additional support grid with 8 × the number of points was used to evaluate the forces during the single-point force calculations for the lattice-dynamics calculations. The PAW projection was applied in reciprocal space (LREAL = .FALSE. in VASP).

For the lattice-dynamics calculations, second- and third-order force constants were computed in 2 × 2 × 2 expansions of the primitive cells. Additional calculations were carried out on cells with expansions and contractions of ±3% about the athermal equilibrium volume, in steps of 1%. The third-order calculations were performed using the 300 K volumes predicted from the quasi-harmonic approximation.

Phonon DoS curves and thermodynamic partition functions were computed using Γ-centred q-point grids with 48 × 48 × 48 subdivisions to integrate the phonon Brillouin zones. The phonon lifetimes for ZnS and CuGaS₂/Cu₂ZnGeS₄ were calculated using Γ-centred q-point grids with 48 × 48 × 48 and 16 × 16 × 16 subdivisions, respectively, and the tetrahedron method was used for interpolation.

The calculations in this work used the facilities of the Supercomputer Center at the Institute for Solid State Physics, University of Tokyo. Some calculations were also performed on the UK Archer HPC facility, accessed through membership of the UK HPC Materials Chemistry Consortium, which is funded by EPSRC Grant No. EP/L000202. We also made use of the Balena HPC facility at the University of Bath, which is maintained by Bath University Computing Services. J.M.S. is funded by an EPSRC Programme Grant (Grant No. EP/K004956/1). A.W. acknowledges support from the Royal Society and the ERC (Grant No. 277757).

Data Access Statement. The crystal structures and phonon data reported in this work are available in an online repository at https://github.com/WMD-group/Phonons, which can be processed using the Phonopy and Phono3py packages available from http://atztogo.github.io/phonopy and http://atztogo.github.io/phono3py. The animations were made using ascii-phonons, available from https://github.com/ajjackson/ascii-phonons.