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

In semiconductors almost all heat is conducted by phonons (lattice vibrations), which is limited by their quasi-particle lifetimes. Phonon-phonon interactions represent scattering mechanisms that produce thermal resistance. In thermoelectric materials, this resistance due to anharmonicity should be maximised for optimal performance. We use a first-principles lattice-dynamics approach to explore the changes in lattice dynamics across an isostructural series where the average atomic mass is conserved: ZnS to CuGaS$_2$ to Cu$_2$ZnGeS$_4$. Our results demonstrate an enhancement of phonon interactions in the multernary materials, and confirm that lattice thermal conductivity can be controlled independently of the average mass and local coordination environments.

In semiconductors almost all heat is conducted by phonons (lattice vibrations), which is limited by their quasi-particle lifetimes. Phonon-phonon interactions represent scattering mechanisms that produce thermal resistance. In thermoelectric materials, this resistance due to anharmonicity should be maximised for optimal performance. We use a first-principles lattice-dynamics approach to explore the changes in lattice dynamics across an isostructural series where the average atomic mass is conserved: ZnS to CuGaS 2 to Cu 2 ZnGeS 4 . Our results demonstrate an enhancement of phonon interactions in the multernary materials and confirm that lattice thermal conductivity can be controlled independently of the average mass and local coordination environments. C  Direct heat to electricity conversion in a thermoelectric device represents an important process for energy generation and efficiency. 1 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. 2 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 2 Te 3 ), 3 loosely bound "rattling" atoms (e.g., clathrates), 4 low dimensionality (e.g., nanorods and quantum dots), 5 and phonon localization (e.g., superlattices and alloys). 6 Since hundreds of materials have been systematically screened, tried, and tested, new ideas are required to move this field forward. [7][8][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. 10 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: 11 a a.walsh@bath.ac.uk 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 12 and Goodman 13 on cation cross-substitution to produce multi-component semiconductors, which was later applied to screening materials for photovoltaics, 14 topological conductivity, 15 and spintronics. 16 We consider the mutation from ZnS → CuGaS 2 → Cu 2 ZnGeS 4 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." 17 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][19][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 P and P3 packages 20,21 using VASP, 22 as the force calculator and the PBEsol 23 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. 24 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. 25 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 F43m (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 2 (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 −1 , respectively. 27 The calculated values from the harmonic phonon dispersion (Figure 1) are 275 and 335 cm −1 , in good agreement. In addition, the phonon dispersion across the Brillouin zone agrees very well with neutron scattering measurements. 26 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 −6 K is in good agreement with both experiment and reported lattice-dynamics calculations. 28 From binary to ternary: CuGaS 2 . The crystal structure of CuGaS 2 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 I42d (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 1 + 2A 2 + 3B 1 + 4B 2 + 7E vibrational branches, of which only the A 2 modes are IR and Raman inactive.
The harmonic phonon dispersion (Figure 1) is in good agreement with previous calculations. 31 From ternary to quaternary: Cu 2 ZnGeS 4 . The kesterite crystal structure adopted by Cu 2 ZnGeS 4 is closely related to chalcopyrite, with the addition of alternating Cu-Zn and Cu-Ge (001) planes. The space group type is I4 (S 4 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 −1 . The measured Raman modes run from 72 to 409 cm −1 , which were also found to be in good agreement with prior DFT calculations using the PBE functional. 35 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 q j, λ (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 −1 K −1 at 30 K to 27 Wm −1 K −1 at 300 K, 36  ). These are compared to similar functions computed from the trace of the group-velocity tensor products ((d)-(f)), the phonon lifetimes ((g)-(i)), and the averaged phonon-phonon interaction strengths ((j)-(l)). On each plot, the functions are overlaid on the background phonon density of states (DoS) curves. A large value of a frequency band with a small DoS implies that those modes have a high value of the accumulated quantity (e.g., higher group velocity, longer lifetime).
latter of which is within a factor of two of our calculated value. For CuGaS 2 and Cu 2 ZnGeS 4 , 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 (κ CGS xx = κ CZGS xx = 10.2 Wm −1 K −1 at T = 300 K) than along the c axis (κ CGS zz = 8.8 Wm −1 K −1 ; κ CZGS xx = 8.0 Wm −1 K −1 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 2 -type variation, and 46% for Cu 2 ZnGeS 4 -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 2 and Cu 2 ZnGeS 4 , 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 10 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 2 and Cu 2 ZnGeS 4 , respectively). Projector augmented-wave (PAW) datasets 37 were employed. A tolerance of 10 −8 eV was applied during the electronic minimization, and geometry optimisations were carried out to a force threshold of 10 −2 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 2 /Cu 2 ZnGeS 4 were calculated using Γ-centred q-point grids with