Unraveling the role of chemical composition in the lattice thermal conductivity of oxychalcogenides as thermoelectric materials

Abstract

Oxychalcogenides represent a large chemical space with potential application as thermoelectric materials due to their low thermal conductivity. However, the nature of this behaviour is still under debate. Understanding the origin of the anharmonicity of these materials is key to developing models that improves the efficiency of thermoelectric materials. In this work, we combine machine learning with first principles calculations to explore oxychalcogenides materials. Machine learning not only accelerates the prediction of the lattice thermal conductivity for large chemical spaces with high accuracy, but also catalyzes the development of design principles to discover new thermoelectric materials. Using this approach, lattice thermal conductivity has been directly connected to the effect of each species in the material, using atomic projections of the scattering rates. The role of the monovalent atom and the lone pair electron for the trivalent cation are discussed in detail. Based on this knowledge, it is possible to connect complex properties such as lattice thermal conductivity with a more manageable synthetic variable such as chemical composition. Using this strategy, we propose promising new oxychalcogenides such as BiOAgSe, which subsequently has been confirmed as having ultra-low lattice thermal conductivity and BiOAg_0.5Cu_0.5Se which presents a lower κ_l and promising properties for thermoelectricity.

---

Jose J. Plata

Jose J. Plata is an Assistant Professor in the School of Chemistry at the University of Seville. Previously, Jose held a postdoctoral position at Duke University (2013–2017) and a Marie Skłodowska-Curie individual fellowship (2018–2020). He earned his PhD in Theoretical and Computational Chemistry from the University of Seville in 2013. His research focuses on the development of high-throughput frameworks for the prediction of materials properties, utilizing state-of-the-art computational techniques.

1 Introduction

Global energy demands require the development of new strategies to improve the efficiency of our energy grid and, simultaneously, accelerate the transition to renewable energies. Thermoelectric, TE, materials are one of the best available solutions when it comes to reducing the dependency on fossil fuels and making other green energy resources more cost efficient. However, the discovery and optimization of new thermoelectric materials face many obstacles. The efficiency of a thermoelectric materials, measured by its figure of merit ZT, rarely surpassed 1 during the last 50 years of the 20th century. This is partially due to the relationship between the transport properties that define ZT,where σ is the electrical conductivity, κ_e and κ_l are the electronic and lattice thermal conductivities, T is the temperature and S is the Seebeck coefficient. Maximizing ZT is a difficult and complex task; Wiedemann–Franz law requires σ to be proportional to κ_e and Pisarenko relation limits the enlargement of S and σ simultaneously.^1,2 For this reason, one of the most common strategies for maximizing ZT is looking for new thermoelectric materials with low κ_l.

The experimental measurement of κ_l or any other of the transport properties involved in ZT represent another obstacle. Most of these properties are extremely sensitive to synthesis conditions, hampering the systematic exploration of large chemical spaces. Despite of these difficulties, thermoelectric materials are experiencing a renaissance during last two decades.³ Atomistic simulations and computational materials science have propelled the discovery of many TE materials with ZT values above 2.⁴ Nevertheless, there is plenty of room for improvement in this area. For instance, the accurate prediction of κ_l requires large and expensive ab initio molecular dynamics, MD,⁵ or the calculation, at least, of 3rd-order interatomic forces constants, IFCs.⁶

Machine learning, ML, has recently emerged as the ideal strategy for reducing computational costs when calculating κ_l. Accurate ML-based potentials have reduced the computational cost of approaches based on MD⁷ and new methods have been developed to extract nth-order IFCs using ML at a reduced cost.^8,9 These new approaches have changed the paradigm when exploring new TE materials. Exploring large chemical spaces is no longer a problem, and the scientific community is focusing on extracting design principles for these large amount of data.¹⁰ However, this last step cannot be underestimated. Connecting real space variables such as chemical composition with reciprocal-space variables which determine transport properties is fundamental for developing new and more efficient TE materials.¹¹

Oxychalcogenides represent one of these large chemical spaces that need to be explored in more detail.¹² Since BiOCuSe was reported as a promising TE material in 2011,^13,14 a large amount of data has been published with respect to their TE performance.^15–17 For instance, different strategies have been proposed to increase their power factor, such as doping,¹⁸ alloying¹⁹ or engineering grain boundaries.²⁰ Torimayama et al. have recently demonstrated that electron transport properties of these materials can be improved through n-type doping instead of p-type doping.²¹ However, the key property that make oxychalcogenides excel as TE materials is their low lattice thermal conductivity. Interestingly, this fact has been attributed to different factors, and there is still some controversy. For instance, Ding et al. described how anharmonic vibrations and structural scattering of phonons are partially caused by in-layer and interlayer off-phase vibrational modes.²² Saha et al. compared the vibrational properties of BiOCuSe and LaOCuSe, finding that the main reason for their different κ_l is the mass difference between Bi and La.^23,24 Additionally, Vaqueiro et al. combined DFT calculations with analysis of in situ neutron diffraction data to demonstrate that weak bonding of copper atoms play a crucial role, reducing the lattice thermal conductivity of these compounds.²⁵ Although all these works present solid proofs, most of them are based on the calculation of the harmonic force constants to study harmonic (phonon dispersion curves) or quasi-harmonic (Grüneisen parameters) properties. Lattice thermal conductivity is a property that stems from the anharmonicity of the material, thus, features such as anharmonic force constants and scattering rates are essential to understanding its behavior. However, until recently, exploring highly computationally demanding properties such a κ_l in large chemical spaces was an unaffordable task.

Here, we present a new approach in which the low thermal conductivity of oxychalcogenides is directly connected to their chemical composition. Machine learning is combined with first principles calculations to accelerate the calculation of κ_l for oxychalcogenides. These results are analyzed through atomic-projected scattering rates-frequency heat maps, unraveling the role of each element in the anharmonicity and lattice thermal conductivity of these materials.

2 Methodology

2.1 Geometry optimization

The atoms and the lattice of the ground state structures were fully relaxed at 0 K with VASP package,^26,27 using projector-augmented wave (PAW) potentials.²⁸ The exchange-correlation functional proposed by Perdew–Burke–Ernzerhof (PBE)²⁹ was combined with Grimme D3 van der Waals corrections³⁰ to obtain the energies. D3 corrections are critical for a good description of slab interactions and lattice parameters (Table 1). Core and valence electrons were selected following standards proposed by Calderon et al.³¹ Plane waves with a cut-off of 500 eV were used for valence electrons and the primitive cells were optimized with a 5 × 5 × 3 k-point mesh. Wavefunction was considered as converged when the energy difference between two electronic steps was smaller than 10⁻⁹ eV. Geometry and lattice vectors were fully relaxed, using an 8 atoms conventional cell, until forces over all atoms were smaller than 10⁻⁷ eV Å⁻¹. An additional support grid for the evaluation of the augmentation charges was included to reduce the noise in the forces.

Table 1 Calculated lattice parameters a and c for all the materials explored in this work. Average lattice thermal conductivity, κ_l, and its components in xx and zz directions at 300 K. Experimental values are included in parenthesis

Material	a (Å)	c (Å)	κ l (W K^−1 m^−1)	κ ^xxl (W K^−1 m^−1)	κ ^zzl (W K^−1 m^−1)
BiOCuSe	3.91	8.97	1.26	1.66	0.47
	(3.92)^32	(8.91)^32	(1.03)^14
BiOCuS	3.85	8.59	1.83	2.20	1.08
	(3.87)^32	(8.56)^32
LaOCuS	3.94	8.45	4.41	5.40	1.64
	(4.00)^33	(8.51)^33
AcOCuS	4.045	8.64	1.70	2.44	0.24
LaOAgS	4.03	9.04	2.45	3.28	0.80
	(4.05)^34	(9.04)^34
LaOLiS	3.98	8.70	7.83	9.76	3.97
BiOAgSe	3.99	9.54	0.87	1.09	0.43
	(3.99)^35	(9.46)^35

---

2.2 Supercell single-point calculations and force constants

Interatomic force constants, IFCs, were calculated using the HiPhive package, which combines the forces calculated for random atomic distortion in supercells with machine learning regression.⁸ The forces were calculated in 5 × 5 × 3 supercell (600 atoms) using the same setup as the one used for the geometry optimizations. The amplitude of the distortions applied to the atoms plays an important role in the calculation of the IFCS so a 2 steps approach was designed.³⁶ First, small random distortions were generated for all the atoms of 3 supercells and second and third-order IFCs were calculated using the HiPhive package. Then, eight new distorted supercells were created superimposing normal modes with random phase factors and amplitudes corresponding to 300 K. The recursive feature elimination RFE, algorithm was used to obtain the force constants from the DFT forces. Even though RFE is more expensive than ordinary least-square regression, this algorithm requires fewer number of structures to converge. Additionally, reducing the number of parameters via RFE also simplifies the model, keeping only the most relevant interaction terms. IFCs were calculated including cutoffs for second, third and fourth-order terms. In order to ensure transferability across compounds, cutoffs were determined based on coordination shells. We have used a wrapper code for the hiPhive program that automates the distorted supercell creation, force calculation using VASP, and the construction of the ML IFCs, which is available for public use.^36,37

2.3 BTE solver

The ShengBTE package was used to solve the Boltzmann Transport Equation, BTE, and obtain lattice thermal conductivity.³⁸ BTE is solved through the iterative approach, which predicts more accurate values than the relaxation time approximation. Lattice thermal conductivity is calculated as a 3 × 3 tensor, κ^αβ_l, so the anisotropy of the system can be analyzed. Oxychalcogenides studied in this work belong to tetragonal space group so only the diagonal values present non-null values with κ^xx = κ^yy Average lattice thermal conductivity, κ_l, will be defined as the average of the diagonal values of the κ^αβ_l tensor. Total scattering rates, Γ_tot, were computed including isotopic and three-phonon scattering. Memory demand and the convergence of κ_l with the number of q-points were balanced using a gaussian smearing of 0.1 eV and a dense mesh of 12 × 12 × 12. Scattering rates heat maps and their projections were obtained combining ShengBTE and Phonopy³⁹ through a high-throughput tool which is available for public use.⁴⁰

3 Results

3.1 BiOCuSe

BiOCuSe is a layered oxychalcogenide that belongs to the tetragonal structural prototype ZrCuSiAs. Its 3D structure alternates fluorite-like [Bi₂O₂]²⁺ slabs with antifluorite-like [Cu₂Se₂]²⁻ slabs along the c-axis.¹² In order to give a good description of the long range interactions between layers, Grimme D3 van der Waals corrections³⁰ have been included, obtaining very good agreement with both, a and c, experimental lattice parameters for BiOCuSe and the other oxychalcogenides included in this work (Table 1). Dispersion curves extracted from ML-regression are also in good agreement with the results obtained using the finite differences approach, full-DFT, where the IFCs are determined directly from the displacement of one or two specific atoms in the supercell (Fig. 1a).

Fig. 1 (a) Comparison of BiOCuSe phonon dispersion curves obtained using the machine-learned force constant potential, ML, via HiPhive with those from the full-DFT method via Phonopy. (b) Comparison of experimental κ_l (orange)¹⁴ with calculated results using a full-DFT approach (green),⁴¹ and ML-based model (blue) between 300 and 900 K.

Root mean square error, RMSE, of the forces rapidly converges with respect to the number of supercells included in the training (Fig. S1a†). Average lattice thermal conductivity, κ_l, follows the same trend and it is aligned with previous values calculated using more computationally-demanding approaches.⁶ Cutoff-interaction distances of 9.75, 6.89, and 3.71 Å for the second-, third-, and fourth-order force constants were found to be sufficient for convergence of the ML-based method. These cutoffs were consistently modified for other materials based on coordination shells. Lattice thermal conductivity was also converged with respect to the q-points grid size (Fig. S1b†). Using a ML-model to extract IFCs and compute κ_l drastically reduces the computational cost without reducing the accuracy of the predictions. Small differences with previous calculated κ_l for BiOCuSe are due to the use of different exchange correlation functionals.⁴¹ Kumar et al. obtained slightly larger lattice parameters using a PBE + U functional which usually is connected to less rigid bonds and lower κ_l. Both, ML-model and full-DFT results, slightly overestimate experimental values,¹⁴ which can be justified because of the polycrystalinity of the samples (Fig. 1b). While calculated κ_l values correspond to single-crystal and free-defects solids, experimental measurements were obtained from polycrystalline samples and most likely in the presence of point defects. Both, grain boundaries and point defects act as scattering centers, reducing the life time of the phonons and κ_l, with respect to single crystal samples.

3.2 Effects of chalcogen substitution

It has already been reported that κ_l follows a monotonous trend down the chalcogen group for BiOCQ (Q = S, Se, Te).^17,42 Both, experimental measurements and calculations have reported that κ_l at 300 K is reduced around 25% when S is substituted by Se and between 30–40% when Se is substituted by Te.^17,42 This is not always the case. For instance, I–III–VI₂ chalcopyrite semiconductors follows a nonmonotonous variation of κ_l, with selenides being the materials with the lowest values.³⁶ Substituting Se by S brings strong modifications in the phonon density of states, DOS (Fig. 2a and b). While vibrational modes of the [Cu₂S₂]²⁻ slab are strongly localized on Cu at lower frequencies (around 2 THz) or S (around 6–8 THz), [Cu₂Se₂]²⁻ slab presents more delocalized vibrations between Cu and Se, participating both atoms in vibrational modes at 2 and 4–5 THz. There are not significant differences between the group velocities (Fig. S2†) of BiOCuSe and BiOCuS. Although there is an important shift to lower energies in the vibrational modes in which the chalcogen atoms participate, when S is substituted by Se, the curvature of the dispersion curves is not strongly modified (Fig. S2a†). Thus, group velocities are not significantly changed either, and their differences in κ_l should be connected to the anharmonicity of the systems. Although they can look similar, there are important differences in the scattering rates of both compounds in the 2–5 THz region of the spectra, where acoustic-optical phonon scattering takes place (Fig. 2c). Increasing the atomic mass of the chalcogen delocalizes and reduces the frequencies of the vibrational modes at which the [Cu₂Q₂]²⁻ slabs are vibrating. This shift maximized the overlap of these optical vibrational modes with the acoustic modes that are mainly localized at the Bi atoms, increasing their scattering rates and reducing κ_l. The same behavior can be extracted from the average accumulative κ_l where there are not big differences in κ_l between BiOCuSe and BiOCuS up to 2 THz approximately (acoustic region) (Fig. 2d). The lattice thermal conductivity of both compounds changes in the region where acoustic-optic scattering takes places.

Fig. 2 (a and b) Total and projected phonon density of states, DOS, of BiOCuSe and BiOCuS. (c) Total scattering rates, Γ_tot for BiOCuSe (green) and BiOCuS (blue) at 300 K. (d) Cumulative κ_l for BiOCuSe (green solid line) and BiOCuS (blue dashed line) at 300 K.

3.3 Effects of the trivalent cation

While there is consensus about the role of the chalcogen atom and how its chemical substitution can slightly modify κ_l for BiOCuQ (Q = S, Se, Te) compounds, the role of the Bi is certainly under debate. Some authors point to the Bi atom's lone pair electron as the main reason for the low thermal conductivity of BiOCuSe,^15,43 but there are also studies that suggest that the Cu-Q vibrations are the ones which contribute to reduce κ_l.²⁵ Partially, these discrepancies come from using only harmonic (phonon dispersion curves) or quasiharmonic (Gruneisen parameters) features to explore a property that is strongly coupled to the anharmonicity of the system. Analyzing the scattering rates is fundamental to understanding the origin of the low lattice thermal conductivity of oxychalcogenides. In order to understand the role of Bi, previous works have compared the vibrational, elastic and electronic properties of the BiOCuQ materials with their La counterparts.^23,24 However, that substitution introduces more than one change in the system: (i) the lone pair is removed and (ii) the mass of the trivalent atom is drastically reduced. Here, the lattice thermal conductivity and anharmonicity of BiOCuS, LaOCoS and AcOCuS are compared in order to decouple the role of the mass and lone pair of the trivalent atom. Substituting Bi by Ac, lone pair electron is removed, but the atomic mass is increased, which complements the information that can be extracted from substituting Bi by La. In addition to the phonon DOS and the calculation of the cumulative κ_l (Fig. 3), scattering rates–frequencies heat maps have been built, including its projection over each element (Fig. 4). These representations not only connect phonon frequencies with a direct measurement of the anharmoncity of the system, but they also unravel the role of each species. Phonon density of states, DOS, of LaOCuS and AcOCuS present similarities with BiOCuSe DOS (Fig. 3a and b). All of them present a region around 2 THz where optics bands are overlapping with acoustic bands, participating Cu and the trivalent cation. However, DOS cannot explain the differences in the in-plane (xx) and perpendicular (zz) cumulative κ_l. Due to the layered nature of these materials, κ^zz_l which is perpendicular to the slab stacking, is smaller than κ^xx_l (Fig. 3c and d). Cumulative κ^zz_l follows the same trend for all three materials, 100% of κ^zz_l is stored up to 2 THz, approximately. The acoustic-optic scattering region around 2 THz includes atoms from both different layers which can be related to the interlayer off-phase vibrational modes described by Ding.²² Larger κ^zz_l values are obtained for LaOCuS, then BiOCuS and finally AcOCuS which points to the atomic mass as the main responsible of this trend. However, the behavior is more complex for κ^xx_l. Cumulative κ^xx_l presents two main contributions, one in the 0–3 THz range and another one around 8–13 THz (Fig. 3c and d). First large contribution to κ^xx_l follows the same trend than κ^zz_l, ranking their values based on the atomic mass of the trivalent atom. In addition to larger group velocities, heat maps present lower scattering rates in the 0–3 THz range for LaOCuS, which is linked to a lower anharmonicity and larger κ_l (Fig. 4). Second contribution to κ^xx_l is entirely due to optical phonons and it is slightly smaller than the first one, but still constitutes a 30–45% of the total value. While acoustic phonons with long mean free paths usually control κ_l, this is one of the few reported examples in which high-frequency optical modes strongly modify κ_l.⁴⁴ Even more interestingly, optical modes contributions to κ^xx_l are different for each material. BiOCuS presents the lowest κ^xx_l because the optical mode contribution to κ_l is proportionally smaller than in LaOCuS and AcOCuS. Due to these high frequency contributions, κ^xx_l is not ranked by the atomic mass of the trivalent atom. This behavior is directly connected to the scattering rates of the vibrational modes at high frequencies (Fig. 4). LaOCuS presents the lower scattering rates for frequencies larger than 7–8 THz, which is consistent with larger contributions to κ_l. On the other hand, BiOCuS exhibits the highest scattering rates at high frequencies, so their contributions to κ_l is reduced. Projected scattering rates show how this high-frequency contribution to κ_l is mainly localized on O atoms vibrations. Thus, the differences in the scattering rates of the high frequency phonons stem from the different nature of the M-O bonds (M = La, Ac, Bi). If BiOCuS is the material with the highest scattering rates at high frequencies, atomic mass cannot be the only variable to consider and its lone pair electron must play a role. It has been proven that the presence of stereochemically active lone pair electrons, LPEs, enhances the anharmonicity of a system and reduce κ_l.⁴⁵ This reduction is due to nonlinear repulsive electrostatic force between LPEs and other neighboring bonds, which lowers the lattice symmetry and promote the appearance of more scattering processes.^46,47 Following the revised lone pair model,⁴⁸ it has been experimentally confirmed the presence of a LPE in BiOCuSe.⁴⁹ Sallis et al. studied the electronic structure of BiOCuSe and LaOCuSe using K-edge X-ray emission spectroscopy, XES, X-ray absorption spectroscopy, XAS and density functional theory. They found that LPE in the BiOCuSe is the result of the formation of a Bi 6s–O 2p bonding state just below the valence band and the stabilization of the Bi 6s–O 2p anti-bonding state due to its hybridization with Bi 6p empty states. This hybridization has also been described through first principles calculations.^50,51 Here, it is important to stress the relevance of the concept of LPE as a bonding state, in order to understand their role reducing κ_l. Although, optical phonons contributions to κ^xx_l are mainly projected on O atoms, the scattering rates of these vibrational modes are larger than in LaOCuS and AcOCuS because of the Bi–O bond nature, which is singularly characterized by the presence of a LPE.

Fig. 3 (a and b) Total and projected phonon density of states, DOS, of LaOCuS and AcOCuS. (c and d) Cumulative κ^xx_l and κ^zz_l for BiOCuSe, LaOCuS, and AcOCuS at 300 K.

Fig. 4 Scattering rates–frequency heat maps for BiOCuSe (top left), LaOCuS (top mid), and AcOCuS (top right) at 300 K. Heat map color represents the density of vibrational modes with similar scattering rates and frequencies. Projections on each element have been depicted below the total contributions.

3.4 The role of Cu

Some previous works have described Cu as the main contributor to the low lattice thermal conductivity of BiOCuQ family.²⁵ Vaqueiro et al. used neutron diffraction experiments and DFT calculations to demonstrate that low-energy vibrational modes on Cu and its weak bonding with the chalcogen are essential to understand their low κ_l. However, is this a singular interaction between Cu and the chalcogens or is it possible to obtain similar results with other elements? Here, the vibrational properties and κ_l of LaOMS (M = Cu, Ag, Li) have been explored for a better understanding of the role of the monovalent cation. While the substitution of Bi by Ac or La do not produce large modification in the phonon DOS of the compounds, substituting Cu by Li or Ag lead to severe modifications of the vibrational structure (Fig. 5a and b).

Fig. 5 (a and b) Total and projected phonon density of states, DOS, of LaOLiS and LaOAgS. (c) Group velocities, v_λ for LaOLiS (blue), LaOCuS (green), and BiOCuS (orange). (d) Cumulative κ_l for LaOLiS (blue solid line), BiOAgS (green dashed line) and BiOCuS (orange dashed line).

Li atoms contributions to the DOS are shifted more than 7 THz, if compared with the projection of Cu atoms in LaOCuS DOS. On the other hand, Ag substitution shift the monovalent contribution to the DOS to even lower frequencies, making it as the main contributor to the acoustic modes. The absence of a heavy monovalent atom such as Cu or Ag also increases the group velocities of the phonons at low frequencies for LaOLiS (Fig. 5c). This behaviour has a clear impact in the cumulative κ_l of LaOLiS which grows rapidly up to 6 W K⁻¹ m⁻¹ below 4 THz (Fig. 5d). Silver substitution produces the opposite effect on κ_l because acoustic-optical scattering region is pushed to even lower frequencies when compared to LaOCuS. That trend can be visualized in the scattering rates heat maps where LaOAgS scattering rates at low frequencies are higher than in LaOCuS and LaOLiS. As expected, La and Ag contributions dominates in this region so the anharmoncitiy comes from interlayer scattering (Fig. 6).

Fig. 6 Scattering rates–frequency heat maps for LaOLiS (top left), LaOAgS (top mid), and BiOAgSe (top right). Heat map color represents the density of vibrational modes with similar scattering rates and frequencies. Projections on each element have been depicted below the total contributions.

3.5 Designing new oxychalcogenides

Understanding the mechanisms that governs κ_l in oxychalcogenides become fundamental for the design of new materials withe potential TE applications. Larger contribution to κ_l corresponds to scattering processes between acoustic and low frequency optical phonons. This scattering can be maximized with heavy monovalent and trivalent atoms and using Se or Te. In addition, second contribution to κ_l at higher frequencies can be minimized using trivalent atoms with LPE that increases the scattering rates of high frequency phonons, mainly localized on O atoms. Taking both considerations into account, BiOAgSe is proposed as an oxychalcogenide with potential ultra-low κ_l. To the best of our knowledge, there are very few recent reports in which BiOAgSe has been synthesized.^35,52,53 Some of these works highlight their optoelectronic properties and potential photoconversion applications.^35,52 Although BiOAgTe should present also a very low κ,^54,55, selenides are usually preferred instead of tellurides because tellurium is a very volatile element which hampers the stability and durability of the material. BiOAgSe has been very recently reported as a promising candidate as TE materials due to its tunable electrical conductivity,⁵⁶ and, during the preparation of this manuscript, its κ_l has been reported using the relaxation time approximation, RTA.^53,57 However, very little has been described about the origin of its ultra-low κ_l. Phonon DOS already shows how there are some important differences with respect previous materials (Fig. 7a). BiOAgSe exhibits κ_l lower than 1.2 W K⁻¹ m⁻¹ at 300 K in both xx and zz directions and an average κ_l of 0.87 W K⁻¹ m⁻¹ (Fig. 7b and c). These values are slightly larger than the ones predicted by Li et al.⁵⁷ because the RTA tends to underestimate κ_l and in very good agreement with recently calculated values obtained by Zhang et al.⁵³ solving the BTE iteratively (Fig. 7b). These low values are due to the effective overlap between acoustic and optical modes mainly projected on the Bi and Ag atoms (Fig. 6 and 7a). Contribution to κ_l at high frequencies is also very small (around 0.3 W K⁻¹ m⁻¹) (Fig. 7c), and it is minimized with respect to other oxychalcogenides explored in this work due to the high scattering rates of the high frequency phonons (Fig. 6). In addition to the single-crystal lattice thermal conductivity calculated by Li et al. and the analysis of its origen presented in this work, lattice thermal conductivity for polycrystalline BiOAgSe is included (Fig. 7d). The effect of nanostructuring on the thermal conductivity is based on the decomposition of the contributions to κ by the phonon mean free path.⁵⁸ In this approach, the value of κ_l corresponding to a particular particle size L is calculated as the cumulative contributions for all mean free paths up to L. Large scattering rates or very short lifetimes produces short mean-free-path so small grains are needed to major reductions of κ_l. However, some reduction can be achieved if the average grain size is reduced to 100 nm, where the value of κ_l is below 0.8 W K⁻¹ m⁻¹. These reduction of κ_l is in agreement with the experimental measurements of polycrystalline samples with an average size of few hundreds of nanometers (0.62 W K⁻¹ m⁻¹ at 300 K).⁵³

Fig. 7 (a) Total and projected phonon density of states, DOS, of BiOAgSe and cumulative κ^xx_l and κ^zz_l at 300 K. (b) Average κ_l of BiOAgSe between 300 K and 1000 K. Open points are used for experimental data⁵³ and filled points for calculated values. The predicted values of this work (blue) are compared with the results calculated by Li⁵⁷ (red) and Zhang⁵³ (purple). (d) Cumulative average lattice thermal conductivity at 300 K from mean-free path contributions up to distance L for BiOAgSe, indicating the effect that nanostructuring would have on their κ_l.

Nanostructuring is not the only approach for reducing the lattice thermal conductivity of solids. Alloying also constitutes a common strategy to reduce κ_l.⁵⁹ This strategy has been implemented in BiOAgSe replacing one of the two Ag atoms by a Cu atom in its unit cell. Here, the lattice thermal conductivity of BiOAg_0.5Cu_0.5Se has been calculated following a similar methodology. Phonon dispersion curves present strong imaginary frequencies which point out that this material is not dynamically stable at 0 K (Fig. S3†). Dynamical stability at finite temperatures has been analyzed using the self-consistent phonon, SCPH, approach implemented in the hiPhive code.^8,9 Using this technique, it is possible to include the anharmonicity effects in the dynamical matrix, obtaining temperature dependent frequencies for each vibrational mode. It can be noticed that imaginary frequencies are completely removed between 200 and 300 K (Fig. 8a), demonstrating the dynamic stability of this material at least above 300 K. BiOAg_0.5Cu_0.5Se presents a κ_l of 0.78 W K⁻¹ m⁻¹ at 300 K which is a 10% smaller than κ_l for BiOAgSe (Fig. 8b). This reduction is mainly due to the alloying effect that promotes scattering process in the [AgCuSe₂]²⁻ layer. This behavior can be demonstrated analyzing cumulative κ_l and the scattering rates (Fig. 8b). First, there are very little changes between κ^zz_l for BiOAgSe and BiOAg_0.5Cu_0.5Se, so the reduction of κ_l is not due a modification of the interaction between the layers. Moreover, the contribution of high frequency vibrational modes to κ^xx_l, which are localized in the Bi–O bonds, is similar or both materials (Fig. 8b). Main differences between both materials comes from the contribution of optical acoustic modes and low energy optical modes to κ^xx_l. This contribution is drastically reduced in BiOAg_0.5Cu_0.5Se because of the presence of an additional atom with a different mass, which increases the scattering rates of the vibrational modes between 0 and 5 THz (Fig. 8c). Similarly to BiOAgSe, polycrystalline samples present lower κ_l, with values around 0.7 W K⁻¹ m⁻¹ at 300 K for an average grains size between 100 and 500 nm. However, a reduction of the thermal conductivity is not the only benefit that is obtained when alloying is applied to thermoelectric materials. The other properties included in ZT, S and σ, can be also tuned by different synthetic variables and strategies such as doping, nanostructuring or alloying. For instance, Xie et al. have found that Cu_1−xAg_xGaTe₂ alloys present an unusual non-parabolic band structure, which is crucial for obtaining a high Seebeck coefficient.⁶⁰

Fig. 8 (a) Phonon dispersion curves for BiOAg_0.5Cu_0.5Se at 300 K. (b) Total and projected phonon density of states, DOS, of BiOAg_0.5Cu_0.5Se and cumulative κ^xx_l and κ^zz_l for BiOAg_0.5Cu_0.5Se at 300 K. (c) Total scattering rates, Γ_tot for BiOAg_0.5Cu_0.5Se at 300 K. (d) Cumulative average lattice thermal conductivity at 300 K from mean-free path contributions up to distance L for BiOAg_0.5Cu_0.5Se, indicating the effect that nanostructuring would have on their κ_l.

Moreover, the substitution of Ag by Cu favors the presence of carriers, which is beneficial to obtaining higher power factors.^61,62

Although the calculation of ZT is out of the scope of this work, it is possible to anticipate that BiOAg_0.5Cu_0.5Se is a potential candidate as thermoelectric material.

4 Conclusions

In this work, we have demonstrated how machine learning not only accelerates the prediction of the lattice thermal conductivity for large chemical spaces with high accuracy but also catalyzes the development of design principles to discover new thermoelectric materials. Lattice thermal conductivity of seven oxychalcogenides with very different compositions have been predicted and analyzed. While most previous theoretical work discussed the origin of the low κ_l of oxychalcogenides based on harmonic properties, we introduce a new approach based on the analysis of phonon DOS, cumulative κ_l, and most importantly scattering rates heat maps and their atomic projections. Chemical substitution for all sites have been considered in order to unravel the role of each chemical specie. Two main contributions to κ_l have been described. Low frequency phonons contribution to κ_l can be minimized with heavy trivalent and monovalent cations. Heavy atoms with similar vibrational frequencies promotes the scattering between acoustic and low-frequency optical modes. We also found that high-frequency optical modes have an important contribution to κ^xx_l. Atomic projections of the scattering rates heat maps show that these modes are mainly localized on the O atoms. Although O atoms are the main contributors to these high-frequency phonons, their scattering rates are sensitive to the trivalent atom that is part of the [M₂O₂]²⁺ slab. It has been demonstrated that the atomic mass of the trivalent atom is not the driven force that modifies the scattering rates of these modes. However, it is the element with the LPE, the one that increases the scattering rates and reduces the contribution of these modes to κ_l. Finally, we have followed these design principles based on chemical composition to propose BiOAgSe and BiOAg_0.5Cu_0.5Se as two of the best oxychalcogenides for TE applications.

Author contributions

The authors confirm contribution to the paper as follows: conceptualization: J. J. P. and A. M. M.; data curation: E. J. B and J. S.; investigation: E. J. B, J. S. and F. L. P.; visualisation: E. J. B.; writing – original draft: J. J. P. R., A. M. M. and J. F. S. All authors reviewed the results and approved the final version of the manuscript.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was funded by the Ministerio de Ciencia e Innovación (PID2019-106871GB-I00). E. J. B., J. J. P., J. S., A. M. M., and J. F. S. thankfully acknowledge the computer resources at Lusitania and the technical support provided by Cénits COMPUTAEX and Red Española de Supercomputación, RES (QS-2021-1-0027). The authors thank Fredrik Eriksson, Petter Rosander and Erik Fransson for discussions on hiPhive self-consistent phonons routines.

Notes and references

1. H.-S. Kim, Z. Gibbs, Y. Tang, H. Wang and G. Snyder, APL Mater., 2015, 3, 041506.
2. L. Xi, S. Pan, X. Li, Y. Xu, J. Ni, X. Sun, J. Yang, J. Luo, J. Xi, W. Zhu, X. Li, D. Jiang, R. Dronskowski, X. Shi, G. Snyder and W. Zhang, J. Am. Chem. Soc., 2018, 140, 10785–10793.
3. T. Zhu, Y. Liu, C. Fu, J. Heremans, J. Snyder and X. Zhao, Adv. Mater., 2017, 29, 1605884.
4. P. Gorai, V. Stevanović and E. S. Toberer, Nat. Rev. Mater., 2017, 2, 17053.
5. C. Carbogno, R. Ramprasad and M. Scheffler, Phys. Rev. Lett., 2017, 118, 175901.
6. J. J. Plata, P. Nath, D. Usanmaz, J. Carrete, C. Toher, M. de Jong, M. D. Asta, M. Fornari, M. Buongiorno Nardelli and S. Curtarolo, npj Comput. Mater., 2017, 3, 45.
7. C. Verdi, F. Karsai, P. Liu, R. Jinnouchi and G. Kresse, npj Comput. Mater., 2021, 7, 156.
8. F. Eriksson, E. Fransson and P. Erhart, Adv. Theory Simul., 2019, 2, 1800184.
9. E. Fransson, F. Eriksson and P. Erhart, npj Comput. Mater., 2020, 6, 135.
10. Y. Iwasaki, I. Takeuchi, V. Stanev, A. Kusne, M. Ishida, A. Kirihara, K. Ihara, R. Sawada, K. Terashima, H. Someya, K.-I. Uchida, E. Saitoh and S. Yorozu, Sci. Rep., 2019, 9, 2751.
11. H. Zhu, C. Xiao and Y. Xie, Adv. Mater., 2018, 30, 1802000.
12. S. Luu and P. Vaqueiro, J. Materiomics, 2016, 2, 131–140.
13. Y. Liu, L.-D. Zhao, Y. Liu, J. Lan, W. Xu, F. Li, B.-P. Zhang, D. Berardan, N. Dragoe, Y.-H. Lin, C.-W. Nan, J.-F. Li and H. Zhu, J. Am. Chem. Soc., 2011, 133, 20112–20115.
14. F. Li, J.-F. Li, L.-D. Zhao, K. Xiang, Y. Liu, B.-P. Zhang, Y.-H. Lin, C.-W. Nan and H.-M. Zhu, Energy Environ. Sci., 2012, 5, 7188–7195.
15. Y.-L. Pei, J. He, J.-F. Li, F. Li, Q. Liu, W. Pan, C. Barreteau, D. Berardan, N. Dragoe and L.-D. Zhao, NPG Asia Mater., 2013, 5, e47.
16. M.-K. Han, Y.-S. Jin, B. Yu, W. Choi, T.-S. You and S.-J. Kim, J. Mater. Chem. A, 2016, 4, 13859–13865.
17. H. Zhu, T. Su, H. Li, C. Pu, D. Zhou, P. Zhu and X. Wang, J. Eur. Ceram. Soc., 2017, 37, 1541–1546.
18. T.-L. Chou, G. Tewari, D. Srivastava, A. Yamamoto, T.-S. Chan, Y.-Y. Hsu, J.-M. Chen, H. Yamauchi and M. Karppinen, Mater. Chem. Phys., 2016, 177, 73–78.
19. C. Barreteau, D. Bérardan, L. Zhao and N. Dragoe, J. Mater. Chem. A, 2013, 1, 2921–2926.
20. G. Li, S. Hao, S. Morozov, P. Zhai, Q. Zhang, W. Goddard and G. Snyder, ACS Appl. Mater. Interfaces, 2018, 10, 6772–6777.
21. M. Toriyama, J. Qu, G. Snyder and P. Gorai, J. Mater. Chem. A, 2021, 9, 20685–20694.
22. J. Ding, B. Xu, Y. Lin, C. Nan and W. Liu, New J. Phys., 2015, 17, 083012.
23. S. Saha, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 92, 041202.
24. S. Saha and G. Dutta, Phys. Rev. B, 2016, 94, 125209.
25. P. Vaqueiro, R. Al Orabi, S. Luu, G. Guélou, A. Powell, R. Smith, J.-P. Song, D. Wee and M. Fornari, Phys. Chem. Chem. Phys., 2015, 17, 31735–31740.
26. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558–561.
27. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
28. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
29. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
30. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
31. C. E. Calderon, J. J. Plata, C. Toher, C. Oses, O. Levy, M. Fornari, A. Natan, M. J. Mehl, G. L. W. Hart, M. Buongiorno Nardelli and S. Curtarolo, Comput. Mater. Sci., 2015, 108, 233–238.
32. A. Kusainova, P. Berdonosov, L. Akselrud, L. Kholodkovskaya, V. Dolgikh and B. Popovkin, J. Solid State Chem., 1994, 112, 189–191.
33. Y. Takano, C. Ogawa, Y. Miyahara, H. Ozaki and K. Sekizawa, J. Alloys Compd., 1997, 249, 221–223.
34. M. Palazzi, C. C. P. Laruelle and J. Flahaut, The Rare Earths in Modern Science and Technology, Boston, MA, 1982, pp. 347–350.
35. J. Gamon, D. Giaume, G. Wallez, J.-B. Labégorre, O. Lebedev, R. Al Rahal Al Orabi, S. Haller, T. Le Mercier, E. Guilmeau, A. Maignan and P. Barboux, Chem. Mater., 2018, 30, 549–558.
36. J. J. Plata, V. Posligua, A. M. Márquez, J. F. Sanz and R. Grau-Crespo, Chem. Mater., 2022, 34, 2833–2841.
37. J. J. Plata, V. Posligua, A. M. Márquez, J. Fdez Sanz and R. Grau-Crespo, hiPhive Wrapper, 2021, https://github.com/NewMaterialsLab/hiPhive_wrapper_NML.
38. W. Li, J. Carrete, N. A. Katcho and N. Mingo, Comput. Phys. Commun., 2014, 185, 1747–1758.
39. A. Togo and I. Tanaka, Scr. Mater., 2015, 108, 1–5.
40. E. J. Blancas, J. J. Plata, A. M. Márquez and J. Fdez Sanz, Scattering Rates Heat Map Plotter, 2022, https://github.com/ErnstBlancas/PRmaps.
41. S. Kumar and U. Schwingenschlögl, Phys. Chem. Chem. Phys., 2016, 18, 19158–19164.
42. H. S. Ji, A. Togo, M. Kaviany, I. Tanaka and J. H. Shim, Phys. Rev. B, 2016, 94, 115203.
43. M. Nielsen, V. Ozolins and J. Heremans, Energy Environ. Sci., 2013, 6, 570–578.
44. Z. Li, H. Xie, S. Hao, Y. Xia, X. Su, M. G. Kanatzidis, C. Wolverton and X. Tang, Phys. Rev. B, 2021, 104, 245209.
45. C. Chang and L.-D. Zhao, Mater. Today Phys., 2018, 4, 50–57.
46. D. Morelli, V. Jovovic and J. Heremans, Phys. Rev. Lett., 2008, 101, 035901.
47. E. Skoug and D. Morelli, Phys. Rev. Lett., 2011, 107, 235901.
48. A. Walsh, D. Payne, R. Egdell and G. Watson, Chem. Soc. Rev., 2011, 40, 4455–4463.
49. S. Sallis, L. Piper, J. Francis, J. Tate, H. Hiramatsu, T. Kamiya and H. Hosono, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 085207.
50. E. Stampler, W. Sheets, M. Bertoni, W. Prellier, T. Mason and K. Poeppelmeier, Inorg. Chem., 2008, 47, 10009–10016.
51. H. Hiramatsu, H. Yanagi, T. Kamiya, K. Ueda, M. Hirano and H. Hosono, Chem. Mater., 2008, 20, 326–334.
52. A. Baqais, A. Curutchet, A. Ziani, H. Ait Ahsaine, P. Sautet, K. Takanabe and T. Le Bahers, Chem. Mater., 2017, 29, 8679–8689.
53. C. Zhang, J. He, R. McClain, H. Xie, S. Cai, L. Walters, J. Shen, F. Ding, X. Zhou, C. Malliakas, J. Rondinelli, M. Kanatzidis, C. Wolverton, V. Dravid and K. Poeppelmeier, J. Am. Chem. Soc., 2022, 144, 2569–2579.
54. M. Mukherjee and A. Singh, ACS Appl. Mater. Interfaces, 2020, 12, 8280–8287.
55. J. He, Z. Yao, V. Hegde, S. Naghavi, J. Shen, K. Bushick and C. Wolverton, Chem. Mater., 2020, 32, 8229–8242.
56. J. Li, C. Zhang, Y. Yan, J. Yang, B. Shi, Y. Wang and Z. Cheng, Comp. Mater. Sci., 2020, 171, 109273.
57. J. Li, W. Zhai, C. Zhang, Y. Yan, P.-F. Liu and G. Yang, Mater. Adv., 2021, 2, 4876–4882.
58. C. Dames and G. Chen, Thermoelectrics Handbook: Macro to Nano, Boca Raton, FL, 2006, p. 421.
59. H. Xie, H. Wang, Y. Pei, C. Fu, X. Liu, G. Snyder, X. Zhao and T. Zhu, Adv. Funct. Mater., 2013, 23, 5123–5130.
60. H. Xie, Y. Liu, Y. Zhang, S. Hao, Z. Li, M. Cheng, S. Cai, G. Snyder, C. Wolverton, C. Uher, V. Dravid and M. Kanatzidis, J. Am. Chem. Soc., 2022, 144, 9113–9125.
61. C. Wang, Q. Ma, H. Xue, Q. Wang, P. Luo, J. Yang, W. Zhang and J. Luo, ACS Appl. Energy Mater., 2020, 3, 11015–11023.
62. Y. Cao, X. Su, F. Meng, T. Bailey, J. Zhao, H. Xie, J. He, C. Uher and X. Tang, Adv. Funct. Mater., 2020, 30, 2005861.

---

Footnote

† Electronic supplementary information (ESI) available. See https://doi.org/10.1039/d2ta02180c

---