Introduction

The thermoelectric effect enables the direct conversion of a temperature gradient into an electric voltage and vice versa1. Hence, there have been many studies to discover new materials that show high-energy conversion efficiency2,3,4. The performance of thermoelectric materials can conveniently be evaluated by the dimensionless figure of merit ZT, which is defined as

$$ZT=\frac{{S}^{2}\sigma T}{{\kappa }_{{\rm{elec}}}+{\kappa }_{{\rm{latt}}}}$$
(1)

where S is the Seebeck coefficient, σ is the electrical conductivity, T is the temperature, κelec is the electronic thermal conductivity, and κlatt is the lattice thermal conductivity. Achieving high Seebeck coefficient and electrical conductivity with low electronic and lattice thermal conductivity is therefore the formula for success. In practice, there are limits in the accessible property space for semiconductors, including the positive correlation between Seebeck coefficient, electrical conductivity, and electronic thermal conductivity5. Thus, low lattice thermal conductivity is one of promising design strategies for higher performing thermoelectric materials.

We have been inspired by observations of the thermoelectric effect in halide perovskites6,7,8, where a record ZT of 0.14 was recently achieved9. However, there is a strong desire to move away from heavy and toxic elements such as Pb. Cs3Cu2I5 is an alternative metal halide that has drawn attention for optoelectronics. Since the report of its application in blue light emitting diodes (LEDs)10, the optical properties of this material has been intensively investigated11,12,13. Cs3Cu2I5 has an unusual crystal structure where [Cu2I5]3− clusters are surrounded and isolated by Cs+ cations. Unlike perovskites, for example, this Pnma structure has no through-bond connectivity between neighboring anionic fragments. The atomic structure of Cs3Cu2I5 was recently linked to its optical properties including self-trapped excitons and a large Stokes shift in the photoluminescence11,13. Complex crystal structures have been associated with low thermal conductivity in the past14. Indeed, a value of 0.15 W m−1 K−1 was recently measured in [Mn(C2H6OS)6]I4, which shares a low-dimensional crystal structure15.

In this study, we assess the thermoelectric potential of Cs3Cu2I5. We calculate the thermal transport using anharmonic lattice dynamics, which highlights an ultra-low conductivity for the perfect crystal that is below the limit of standard models. We further investigate electron transport through an ab initio description of carrier scattering. The same low-dimensional crystal structure that gives rise to ultra-low thermal transport supports efficient electron transport with a room temperature electron mobility of 18.2 cm2 V−1 s−1. We predict that Cs3Cu2I5 exhibits the characteristics of an ideal thermoelectric semiconductor—a phonon-glass electron-crystal—with an accessible ZT of greater than 2.

Results

Basic bulk properties

The ground-state Pnma crystal structure of Cs3Cu2I5 is illustrated in Fig. 1a. Cs atoms reside at 4c and 8d Wyckoff positions, Cu atoms reside at 4c, and I atoms reside at three distinct 4c and 8d. An associated electron localization function (ELF) plot could be found in Supplementary Fig. 1. The atomic structure is similar to low-dimensional perovskite-inspired materials such as Cs4PbBr6 where metal halide anion clusters are surrounded and isolated by cation counter parts (effective 0D connectivity)16,17. However, due to distinct size and orbital configuration of Cu(I) 3d104s0, it does not form octahedral clusters like Pb(II) 5d106s218. The calculated equilibrium lattice constants of a0 = 10.06 Å b0 = 11.53 Å and c0 = 14.15 Å match well with values determined from X-ray diffraction (\({a}_{0}^{\exp }=10.19\) Å \({b}_{0}^{\exp }=11.66\) Å and \({c}_{0}^{\exp }=14.40\)Å)10. By fitting the energy-volume relationship to the third-order Birch-Murnaghan equation-of-state19, we obtained a bulk modulus of B0 = 10.6 GPa for Cs3Cu2I5, which is even softer than Cs4PbBr6 (B0 = 12.1 GPa)20. For comparison, the bulk modulus of Si is ~100 GPa.

Fig. 1: Crystal and band structures.
figure 1

a Crystal structure, b electronic band structure (the valence band maximum is set to 0 eV), and c phonon band structure of Cs3Cu2I5. Atom projected density of states are plotted on the right side of electronic and phonon band structures, where Cs, Cu, and I contributions are shown in green, blue, and purple lines, respectively.

The electronic band structure is shown in Fig. 1b. This has been calculated using density functional theory (DFT) with the hybrid HSE06 functional using Wannier interpolation21,22,23. The upper valence band, formed of Cu 3d and I 5p orbitals, is found to be flat in reciprocal space (i.e., short-range interactions in real space). At the same time, the lower conduction band is dispersive and formed of Cu 4s and I 5s (see Supplementary Fig. 2). The associated electron effective mass is 0.3 me, which is comparable to other wide band gap semiconductors such as In2O3 and ZnO. The behavior contrasts with other 0D metal halides that have a flat conduction band dispersion (e.g., Cs4PbBr6) where the lower conduction band is formed from Pb 6p)20. The dispersive nature and isotropic overlap of the Cu 4s orbitals contributes to the low electron effective mass in Cs3Cu2I5. For Cs3Cu2I5, the direct band gap of 3.82 eV at Γ point is obtained from HSE06 calculation, which matches well with measured values for the optical absorption onset12.

The phonon dispersion of Cs3Cu2I5 is shown in Fig. 1c. The lattice vibrations were calculated using the finite-displacement method and the PBEsol functional due to its good description of the crystal structure and phonon properties of semiconductors24,25,26. Since the primitive cell contains 40 atoms, there are 120 phonon bands where 3 acoustic branches are followed by a series of localized optic modes up to 7 THz. We find no imaginary frequencies, which confirms dynamic stability for this phase of Cs3Cu2I5. There is a dense block of optic modes that run from 0.5–2.5 THz involving the collection motion of group atoms, such as relative tilting of [Cu2I5]3− clusters. As a result, acoustic branches are strongly suppressed < 0.5 THz. The presence of this dense block of low-energy phonons just above the acoustic branch creates many energy and momentum allowed scattering pathways that significantly reduce the lifetime of the heat carrying acoustic modes. Such a high density of low-frequency modes also explains the small bulk modulus. Inspection of the eigenvectors at Γ point shows that acoustic modes are simple translation of the entire lattice sites, while the highest frequency optic mode involves Cu–Cu stretching within each cluster. Zoomed electronic and phonon band structures are shown in Supplementary Fig. 3.

Ultra-low lattice thermal conductivity

Next we consider the lattice thermal conductivity (κlatt) by taking into account three-phonon scattering. This approach should provide an upper limit to the thermal transport as additional scattering events can be active in real crystals. Three-phonon scattering is subject to energy and momentum conservation rules for phonon creation and annihilation; however, the flat nature of the phonon dispersion for Cs3Cu2I5 opens up many allowed channels. Owing to structural anisotropy, the conductivity along three principal directions (κxx, κyy, and κzz) differs, as shown in Fig. 2a. The largest values are found along the b axis, which is parallel to the base of the trigonal-planar Cu atoms, and has larger group velocities. The isotropically-averaged conductivity (κiso) is 0.018 W m−1 K−1 at room temperature. We tested both the single mode relaxation time approximation (RTA) and direct solution of linearized phonon Boltzmann equation (LBTE)27, which are in good agreement. The LBTE value is ~10% lower across the relevant temperature range (see Supplementary Fig. 4).

Fig. 2: Lattice thermal conductivity analysis.
figure 2

a Lattice thermal conductivity, calculated from direct solution of linearized phonon Boltzmann equation, of Cs3Cu2I5 as a function of temperature. The red (κxx), orange (κyy), and sky blue (κzz) dotted lines indicate direction dependent lattice thermal conductivity that are parallel to a-, b-, and c-axes of the primitive cell. The blue triangle markers (κiso) show the isotropic average of the thermal conductivity tensor (κiso = [κxx + κyy + κzz]/3). \({\kappa }_{\min }\) refers to the minimum lattice thermal conductivity in a diffuson-mediated transport model. b Phonon lifetime (τλ) at 300 K as a function of mode frequency (νλ). The color gradient is given by phonon band index.

The lattice thermal conductivity is among the lowest reported for bulk crystals. For instance, κiso of CH3NH3PbI328, Sn2S329, and Cs2BiAgBr630 was reported to be 0.05, 0.06, and 0.07 W m−1 K−1, respectively, at similar levels of theory. This values is also lower than the minimum lattice thermal conductivity (\({\kappa }_{\min }\)) expected from diffuson-mediated thermal transport31, which can be estimated from:

$${\kappa }_{\min }=0.76{k}_{{\rm{B}}}\frac{{v}_{{\rm{s}}}}{{V}_{{\rm{a}}}^{2/3}}$$
(2)

where kB is Boltzmann constant, vs is the speed of sound, and Va is the volume of unit cell per atom. The model predicts \({\kappa }_{\min }\) of 0.1 W m−1 K−1 for Cs3Cu2I5. To understand the origin of such a low-lattice thermal conductivity, below the prediction of the analytical model, we analyzed the modal contributions to the net transport. Extremely short acoustic phonon lifetimes of sub 12 ps (Fig. 2b) result in short phonon mean-free paths of below 24 nm. Thus, we conclude that high-acoustic phonon-scattering rates are the main origin of the ultra-low lattice thermal conductivity in Cs3Cu2I5 despite their moderate group velocity and heat capacity. A full breakdown of the modal contributions is provided in Supplementary Fig. 5. The dense block of low-energy phonons have similarity to the states created by the CH3NH3 cation in CH3NH3PbI332. Although a similar computational approach has been shown to give good agreement with experimental measurements on the thermoelectric propoerties of metal halides such as CsSnI333, it should be noted that perturbation theory based on three-phonon scattering is pushed to its limits of applicability in such cases34,35. We propose that Cs3Cu2I5 is ideal for testing more sophisticated treatments of anharmonicity including self-consistent phonon theory36,37.

Electronic transport properties

Given the wide electronic band gap of 3.82 eV, pristine Cs3Cu2I5 is an insulator as there is no significant thermal population of electrons and holes. However, high carrier concentrations (~1021 cm−3) can be achieved in such wide band gap semiconductors by doping (e.g. Al or Ga doped ZnO38,39 and Zr or W doped In2O340,41), which is an active research topic in transparent conductors42. Motivated by the small electron effective mass, we performed electron transport simulations under n-type doping conditions, in the range 1016–1020 cm−3. Our preliminary calculations suggest that Ba\({\,}_{{\rm{Cs}}}^{+}\) has a shallow donor level, which could support an n-type doping density. The associated electrical transport properties are shown in Fig. 3. The electronic thermal conductivity remains lower than the lattice thermal conductivity for ne≤ 1018 cm−3, which means that thermal losses through phonons will be more critical for thermoelectric applications. Once the carrier concentration exceeds 1019 cm−3, the electronic component becomes dominant.

Fig. 3: Transport property analysis.
figure 3

Electronic transport properties of Cs3Cu2I5 as a function of temperature and carrier concentration. κ stands for the thermal conductivity, where κlatt is contribution of phonons and κelec is contribution of electrons. σ stands for the electrical conductivity, S stands for the Seebeck coefficient, and S2σ denotes the power factor.

The electron mobility of Cs3Cu2I5 is found to be limited by optical phonon scattering, followed by ionized impurity and acoustic deformation potential scattering (see Supplementary Fig. 6a), which is consistent with the behavior of other polar semiconductors43,44. The electron relaxation times and associated materials properties are given in Table 1. The electron mobility at room temperature is 18.2 cm2 V−1 s−1 for an n-type doping concentration of 4 × 1018 cm−3. The simpler constant-relaxation time approximation (CRT), assuming the standard τ = 10−14 s, predicts an overestimated value of 60.1 cm2 V−1 s−1 (see Supplementary Fig. 6b). The electrical conductivity increases from around 3.5 to 2.3 × 104 Sm−1, while the absolute value of Seebeck coefficient decreases from around 750 to 60 μV K−1 in the given doping range at room temperature. The power factor (S2σ) is maximized around 1020 cm−3 at a value of 150 μW m−1 K−2.

Table 1 Three types of electron scattering mechanism are considered as implemented in AMSET44: (i) acoustic deformation potential (ADP); (ii) ionized impurity (IMP); (iii) polar-optical phonon (POP). τT are the lifetimes for a given scattering process at T in unit of fs where electron concentration is assumed to be 4 × 1018 cm−3, \(\alpha_v^{\text{VBM}}\) & \({\alpha }_{v}^{{\rm{CBM}}}\) are the volume deformation potential of valence band maximum (VBM) and conduction band minimum (CBM) in unit of eV, B0 is the bulk modulus in unit of GPa, Z is the charge of electron donor dopant, ϵ0 is the static dielectric constant, ϵ is the electronic contribution to the static dielectric constant, and ωpo is the polar-optical phonon frequency in THz.

Thermoelectric performance prediction

By combining the results on thermal and electron transport, the predicted ZT of Cs3Cu2I5 as a function n-type doping concentration and temperature is shown in Fig. 4a. Here, the temperature range is set to 600 K as prediction for temperature region near the melting temperature of 663 K11 will not be reliable due to possibility of evaporation and mass loss, which is beyond limit of our simulation techniques. At 600 K maximum ZT is predicted to be 2.57 where a n-type doping concentration of 4 × 1018 cm−3 is assumed. Above this doping level, electronic thermal conductivity causes ZT to decrease. When the anisotropic nature of transport properties along the principal axes is considered, ZT can be maximized along c with a ZT of 3.35 due to a combination of lower lattice thermal conductivity and higher power factor. A full analysis of the anisotropy is found in Supplementary Fig. 7. The maximum ZT of Cs3Cu2I5 predicted is the highest for any metal halide compound. While the power factor is moderate, ZT is enhanced by the ultra-low lattice thermal conductivity. Although there have been attempts to make thermoelectric devices using metal halide perovksites8, the performance has been limited thus far. A value of ZT = 0.14 was achieved for CsSnI37, with a similar ZT = 0.13 for CH3NH3SnI36. On the other hand Pb-based halide perovskites have shown negligible thermoelectric response6,8. The predicted ZT of Cs3Cu2I5 is comparable to other emerging thermoelectric compounds. For instance, the maximum ZT of CsBi4Te6, Zn4Sb3, and, SnSe are 0.8 at 225 K45, 1.4 at 750 K46, 2.6 at 923 K47, respectively.

Fig. 4: ZT and efficiency prediction.
figure 4

a Predicted ZT values of Cs3Cu2I5 as a function of temperature and carrier concentration. b Estimated thermodynamic efficiency (η, brown triangles) and average ZT (\(\overline{ZT}\), dark green circles) at electron concentration (ne) of 4 × 1018 cm−3 over temperature range of TCTH. TH is fixed as 600 K. ZT at the given temperature is indicated by navy squares.

The efficiency of a thermoelectric device (η) is not determined by maximum ZT at a single temperature but depends on the average ZT (\(\overline{ZT}\)) over a wide temperature range following

$$\eta =\frac{{T}_{{\rm{H}}}-{T}_{{\rm{C}}}}{{T}_{{\rm{H}}}}\frac{\sqrt{1+\overline{ZT}}-1}{\sqrt{1+\overline{ZT}}+\frac{{T}_{{\rm{C}}}}{{T}_{{\rm{H}}}}}$$
(3)

where TH and TC are the temperature of hot side and cold side of the device, and \(\overline{ZT}=\frac{1}{({T}_{{\rm{H}}}-{T}_{{\rm{C}}})}\mathop{\int}\nolimits_{{T}_{{\rm{C}}}}^{{T}_{{\rm{H}}}}ZTdT\). To further evaluate the performance of n-type Cs3Cu2I5-based thermoelectric devices, we plot efficiency over temperature gradient where TH is fixed as 600 K in Fig. 4b. The average \(\overline{ZT}\) of 1.67 leads to a predicted efficiency of 14.9% over a 300–600 K temperature gradient, where ne is 4 × 1018 cm−3.

Beyond the performance of the bulk metal halide, existing crystal engineering approaches could also be applied to Cs3Cu2I5. Cation (e.g., Cs1−xRbx) and/or anion (e.g., I1−xBrx) substitution is one way to tune ZT, which has been widely adopted to modulate optoelectronic properties in metal halide perovskites11,48,49,50. Such ion substitution could further reduce the lattice thermal conductivity due to the additional alloy-scattering channels51,52. In practice, recent reports have shown improved values of ZT by partially replacing I in CsSnI3, and achieved ZT of 0.14 for CsSnCl3−xIx9 and 0.15 for CsSnBr3−xIx53.

Discussion

We have demonstrated the potential of a metal halide based on a low-dimensional bonding network for heat-to-electricity conversion. Cs3Cu2I5 is predicted to support efficient electron transport with slow heat transport. An ultra-low lattice thermal conductivity of < 0.1 W m−1 K−1 is predicted, among the lowest reported for a crystalline solid. Owing to its effective thermal insulation, the maximum ZT of n-type Cs3Cu2I5 is predicted to be 2.57 at 600 K, which is over 10 times larger than the current record of 0.15 in CsSnBr3−xIx and close to the champion value of 2.6 achieved in SnSe. By alloying or nanostructuring, the ZT may be even further enhanced. Our work therefore provides an alternative avenue for developing high-performing thermoelectric materials. Those materials with low-dimensional crystal structure and dispersive band structure (e.g., band edges with s-s orbital overlap) could be attractive candidates. The main bottleneck to realizing such a high ZT will be achieving appropriate carrier concentrations. Further investigations are required to enhance the electrical conductivity of multinary metal halides beyond their natural limits.

Methods

Thermal transport

The thermal physics was simulated from lattice dynamics calculations considering the phonon eigenvectors, eigenvalues, and lifetimes due to three-phonon scattering. First, harmonic phonon analysis was performed using the Phonopy54 code with a finite-displacement step of 0.01 Å. A 2 × 2 × 2 supercell expansion of the relaxed unit cell, containing 320 atoms, was employed. In total 36 symmetry-reduced displacements were considered. Then, anharmonic calculations were performed using the Phono3py55 code with a finite-displacement step of 0.03 Å. For the third-order force constants, 7818 symmetry-reduced displacements were considered in the 40 atom unit cell. Γ-centered q-meshes of 24 × 24 × 18 and 8 × 8 × 6 were used to compute the phonon density of states and lattice thermal conductivity. Symmetrization of the force constants was applied where the maximum drift value of the third order force constants is 0.73 eV Å−1 before symmetrization. All norms of the atomic forces sets and a convergence test on our thermal transport calculations over q-meshes can be found in Supplementary Fig. 8.

The lifetimes are solved through a direct solution of the linearized Boltzmann transport equation (LBTE)27. Then the macroscopic lattice thermal conductivity κlatt is obtained by summing over phonon band indices (νλ) and wavevectors (q), while normalizing by the volume of the unit cell:

$${\kappa }_{{\rm{latt}}}=\frac{1}{N{V}_{0}}\sum _{\lambda }{\kappa }_{\lambda }=\frac{1}{N{V}_{0}}\sum _{\lambda }{C}_{\lambda }{{\bf{v}}}_{\lambda }\otimes {{\bf{v}}}_{\lambda }{\tau }_{\lambda }$$
(4)

where N is the number of unit cells in the crystal (equivalent to the number of wavevectors included in the Brillouin zone summation) and V0 is the volume of the crystallographic unit cell. κλ contains the product of the modal heat capacity (Cλ), group velocity (vλ), and phonon mean-free path (vλ × τλ, where τλ is the phonon lifetime). Graphical analyses of the modal contributions to lattice thermal conductivity and the atomic force sets used for the thermal conductivity calculations were performed using the Phono3py-Power-Tool code56.

Electron transport

The transport of electrons was simulated by considering three scattering processes, with a characteristic scattering time following Matthiessen’s rule:

$$\frac{1}{{\tau }_{e}}=\frac{1}{{\tau }^{{\rm{ADP}}}}+\frac{1}{{\tau }^{{\rm{IMP}}}}+\frac{1}{{\tau }^{{\rm{POP}}}}$$
(5)

Each component of τe was calculated from Fermi’s golden rule:

$${\tilde{\tau }}_{n{\bf{k}}\to m{\bf{k}}+{\bf{q}}}^{-1}=\frac{2\pi }{\hslash }| {g}_{nm}({\bf{k}},{\bf{q}}){| }^{2}\delta ({\varepsilon }_{n{\bf{k}}}-{\varepsilon }_{m{\bf{k}}+{\bf{q}}})$$
(6)

where is the reduced Planck constant, ε is the electron energy, δ is the Dirac delta function and g is the coupling matrix element. Matrix elements due to the acoustic deformation potential (ADP), ionized impurity (IMP), and polar-optical phonon (POP) mechanisms are considered:

$${g}_{nm}^{{\rm{ADP}}}({\bf{k}},{\bf{q}})={\left[\frac{{k}_{{\rm{B}}}T{\alpha }_{v}^{2}}{{B}_{0}}\right]}^{\frac{1}{2}}\langle {\psi }_{m{\bf{k}}+{\bf{q}}}| {\psi }_{n{\bf{k}}}\rangle$$
(7)
$${g}_{nm}^{{\rm{IMP}}}({\bf{k}},{\bf{q}})={\left[\frac{{n}_{{\rm{ii}}}{Z}^{2}{e}^{2}}{{\epsilon }_{0}}\right]}^{\frac{1}{2}}\frac{\langle {\psi }_{m{\bf{k}}+{\bf{q}}}| {\psi }_{n{\bf{k}}}\rangle }{{\left|{\bf{q}}\right|}^{2}+{\beta }^{2}},$$
(8)
$${g}_{nm}^{{\rm{POP}}}({\bf{k}},{\bf{q}})={\left[\frac{\hslash {\omega }_{{\rm{po}}}}{2}\left(\frac{1}{{\epsilon }_{\infty }}-\frac{1}{{\epsilon }_{0}}\right)\right]}^{\frac{1}{2}}\frac{\langle {\psi }_{m{\bf{k}}+{\bf{q}}}| {\psi }_{n{\bf{k}}}\rangle }{\left|{\bf{q}}\right|}$$
(9)

The formalism is implemented in the ab initio scattering and transport (AMSET) package44. The required material parameters, including dielectric constants ϵ, phonon frequency ω, bulk modulus B0, deformation potential α, were all calculated from first-principles. The associated Seebeck coefficient, electrical conductivity, and electronic component of the thermal conductivity are calculated from the Onsager transport coefficients57. An interpolation factor of 10 was used during AMSET routine.

Density functional theory

All of the underlying total energy and electronic structure calculations were performed based on Kohn–Sham density functional theory58, where periodic boundary conditions are considered to represent the extended crystal. The Vienna Ab Initio Simulation Package (VASP)59,60 was used with Projector augmented-wave (PAW)61,62 method where the valence states of Cs, Cu, and I are treated explicitly by 9(5s25p66s1), 17(2p63d104s1), and 7(5s25p5) electrons, respectively.

The Perdew–Burke–Ernzerhof exchange-correlation functional revised for solids (PBEsol)63 with the plane-wave kinetic energy cutoff of 700 eV was used to optimize crystal structure including the internal positions. A 3 × 3 × 2 Γ-centered k-mesh was adopted for structure optimization, where the convergence criteria were set to 10−6 eV and 10−3 eV Å−1 for total energy and atomic forces, respectively. This crystal structure was then used consistently for all transport calculations.

For thermal transport, the same DFT setup was employed and the Brillouin zone was sampled with 2 × 2 × 1 k-mesh for the 2 × 2 × 2 supercell. For electron transport, the required Kohn–Sham wavefunctions (ψ) were generated using the hybrid functional HSE0664,65. Here we used a dense 6 × 6 × 4 Γ-centered k-mesh for better sampling of the electronic Brillouin zone, but lowered the kinetic energy cutoff to 400 eV in order to reduce the computational load. A convergence test on our electronic transport calculations over k-meshes and interpolation factors can be found in Supplementary Fig. 9.