Thermal transport in Ag₈TS₆ (T = Si, Ge, Sn) argyrodites: an integrated experimental, quantum-chemical, and computational modelling study

Abstract

Argyrodite-type Ag-based sulfides combine exceptionally low lattice thermal and high ionic conductivity, making them promising candidates for thermoelectric and solid-state energy applications. In this work, we studied Ag₈TS₆ (T = Si, Ge, Sn) argyrodite family by combining chemical-bonding analysis, lattice vibrational properties simulation, and experimental measurements to investigate their structural and thermal transport properties. Furthermore, we propose a two-channel lattice-dynamics model based on Grüneisen-derived phonon lifetimes and compare it to an approach using machine-learned interatomic potentials. Both approaches are able to predict thermal conductivity in agreement with experimental lattice thermal conductivities along the whole temperature range, highlighting their potential suitability for future high-throughput predictions. Our findings also reveal a relationship between bond heterogeneity arising from weakly bonded Ag⁺ ions and occupied antibonding states in Ag–S and Ag–Ag interactions and strong anharmonicity, including large Grüneisen parameters, and low sound velocities, which are responsible for the low lattice thermal conductivity of Ag₈SnS₆, Ag₈GeS₆, and Ag₈SiS₆. We furthermore show that thermal and ionic conductivities in all three compounds are independent of each other and can likely be tuned individually.

---

Introduction

To reduce the enormous waste of heat in energy generation, thermoelectric (TE) materials offer a promising solution for energy saving and environmental protection. They can convert heat into electricity or vice versa. The thermal conductivity of a material is crucial for its thermoelectric efficiency, and a lower thermal conductivity results in higher efficiency. For example, several argyrodites (e.g., Ag₈GeSe₆, Cu₇PSe₆, Ag₈SiSe₆ and Ag₈SnSe₆) are known for their high ionic conductivity, and many others have been investigated as potential thermoelectrics.^1–8 Halogen-free argyrodites have a general chemical formula of A_(12−n)/m^m+Tⁿ⁺Q₆²⁻(A = Ag, Cu; T = Si, Ge, Sn; and Q = S, Se and Te).^2,8–11 While high ionic conductivity could be problematic for the stability of a thermoelectric device, we and others have demonstrated that thermal and ionic conductivity of some Ag⁺ and Cu⁺ based argyrodites (e.g., Ag₈GeSe₆, Ag_8−xCu_xGeS₆, and Cu₇PSe₆) are not directly correlated with each other and can also be tuned independently.^8,12–14 A similar situation might be expected for the canfieldite Ag₈SnS₆ and Ag₈SiS₆, which are argyrodite family members and isovalently substituted variants of Ag₈GeS₆.

The canfieldite (Ag₈SnS₆) shows promising thermoelectric properties. Shen and co-workers evaluated its lattice thermal conductivity and its crystal structure in detail, finding an orthorhombic Pna2₁ crystal structure at room temperature.¹⁵ Slade's study reported an additional orthorhombic Pmn2₁ phase at 120 K.¹⁶ Additionally, a previous study pointed out the importance of thermal transport via a diffusive transport mechanism (i.e., random-walk-like rather than wave-like).¹⁷ All previous studies suggested that the weakly bonded Ag⁺ ions contribute to the low lattice thermal conductivity in the canfieldites Ag₈TS₆ (T = Si, Ge, Sn).^13,15,16 However, a complete understanding of the connection between lattice thermal conductivity, ionic conductivity, and their correlation with bonding, anharmonicity, and elastic properties remains unexplored for all three compounds.

Several models have been developed to estimate lattice thermal conductivity with limited computational resources, each offering varying degrees of mathematical complexity and accuracy. However, no existing model is both computationally efficient enough for high-throughput studies and reliably accurate across the entire temperature range. Traditional models such as the one by Slack^18–20 take into account the importance of acoustic phonons and elastic properties, often providing a temperature-dependent lattice thermal conductivity (κ_L), but occasionally with overestimated values. The Cahill²¹ and Agne²² models are alternative approaches, particularly for disordered or amorphous materials, by estimating the minimum thermal conductivity based on random-walk theory. As these models capture the diffusive heat transport limit, they cannot predict the correct temperature behaviour of thermal transport over the whole temperature range when thermal transport via phonons is important.²³

Machine learning (ML) approaches have gained popularity due to their ability to predict κ_L for certain compounds at a reasonable computational cost.^24–26 The accuracy of ML models depends on the quality of the data used to train the models, which can also limit their application. Accurate phonon properties require well-converged quantum chemical calculations. Training models from scratch for each composition makes the high-throughput use of such models unfeasible. However, cheaper, pre-trained alternatives, so-called foundation machine learned interatomic potentials (MLIP), have recently emerged. Foundation models already offer a cost-effective alternative to ab initio calculations of harmonic phonons.²⁷ It has also recently been shown that they can reproduce thermal conductivity acceptably for simple binary systems and are compatible within a factor of 2 with ab initio results. With the help of a few additional data points, they can sometimes be fine-tuned for an accurate reproduction of thermal conductivity.²⁸ Once accurately trained, MLIPs can be used to predict and investigate the lattice thermal conductivity without performing expensive full ab initio calculations of phonon lifetimes based on the relaxation time approach, as implemented in phono3py,^29,30 or without using the costly ab initio Green-Kubo approach.^31–33

Motivated by the interesting ionic and thermal transport properties and the open questions concerning thermal conductivity models, in this work, we go beyond the simple investigation of the three systems Ag₈TS₆ (T = Si, Ge, Sn) with ab initio and experimental approaches and attempt to validate a comparably low-cost, fully ab initio model for thermal conductivity that might be suitable for high-throughput investigations. We build on the recently introduced two-channel model by Xia³⁴ that incorporates both phonon and diffuson contributions from harmonic phonons and assumes that each phonon lifetime is half of its vibration period. The Xia model simplifies the full lattice-dynamics approach introduced by Simoncelli et al.³⁵ and is also connected to the analytical two-channel model by Bernges et al.¹⁴ which can be used to fit experimental data. One drawback of the model by Xia is that it has a simplified estimation of phonon lifetimes. To improve the description of the phonon–phonon scattering of each phonon mode, we combine Bjerg's³⁶ model for computing phonon lifetimes (τ) based on the ideas of Slack, and Xia's two-channel model. This offers a more versatile framework for predicting and analyzing heat conduction, particularly in materials with significant Grüneisen parameters or large unit cells. We also compare it with lattice dynamics calculations based on a foundational machine-learned interatomic potential.

This study fulfils three purposes. First, we demonstrate the clear connection between the bonding properties and lattice thermal conductivity for all three compounds. Essentially, analysing the bonding properties is sufficient to conclude that all three compounds have similar lattice thermal conductivity. Furthermore, we demonstrate that cheap ab initio methods, partly combined with machine learning, can analyse and predict lattice thermal conductivity with high accuracy, potentially enabling high-throughput predictions in the future. Lastly, we investigate the relationship between thermal and ionic conductivity.

Therefore, we begin with a detailed quantum-chemical analysis of the bonding in Ag₈TS₆ (T = Si, Ge, Sn). Next, we analyse the harmonic phonon properties, including sound velocity (v), the Debye temperature (Θ), and the volume-dependent Grüneisen parameters (γ), using both experimental and theoretical methods, and connect these to the bonding analysis. Furthermore, we use the two aforementioned approaches (Grüneisen-based lifetime estimation and foundation model) to predict lattice thermal conductivity and reproduce experimental results. Based on an accurate model of experimental thermal conductivity results and ionic conductivity measurements, we demonstrate that ionic and thermal conductivity are independent in Ag₈TS₆ (T = Si, Ge, Sn). By doing so, we also demonstrate the importance of the diffuson channel for these compounds. By integrating bonding analysis, phonon property prediction, and advanced modelling techniques, we aim to establish a robust framework for predicting the thermal conductivity of inorganic materials, with implications for high-throughput material discovery.

Results and discussion

Structural description and X-ray diffraction

Single crystal X-ray diffraction reported by Slade et al.¹⁶ revealed that the canfieldite Ag₈SnS₆ presents two phase transitions: at low-temperature (120 K), the structure adopts an orthorhombic phase (space group Pmn2₁) (Fig. 1a and b). Upon heating, it transforms to another orthorhombic phase (space group Pna2₁) (Fig. 1c–e), which remains stable at room-temperature (RT). At high-temperature, a second transition takes place around 460 K, where the structure transforms from orthorhombic Pna2₁ to the cubic phase with space group F4̄3m. The first transition from Pmn2₁ to Pna2₁ is rather unusual, as the symmetry is lowered upon heating. In the case of the powder sample, however, no change in diffraction patterns was observed below 120 K in their study.¹⁶ For the related compounds Ag₈GeS₆ and Ag₈SiS₆ (Fig. S1b and c), only the orthorhombic Pna2₁ structure has been reported at room temperature.^9,37,38 Fig. 1e and S1d, show that the Pna2₁ structures (RT) contain chains of alternating, corner-sharing Ag⁺ and T tetrahedra (T = Si, Ge, Sn) along b. Along a, chains of corner-sharing, alternating T tetrahedra and trigonal-planar-coordinated Ag⁺ sites are found. Between these chains, linearly coordinated and threefold-coordinated Ag⁺ sites are present, forming a complex network. A detailed report of the coordination environments of all the argyrodites studied here is presented in Section S2 in the SI. Key results will be discussed as part of the bonding analysis.

Fig. 1 Crystal structure of the Ag₈SnS₆ compound in the (a and b) orthorhombic Pmn2₁ space group at 120 K (low-temperature) and (c–e) orthorhombic Pna2₁ space group at room-temperature. These crystal structures form the basis for the computational analysis discussed in the following sections. The crystal structures for the Ag₈SiS₆ and Ag₈GeS₆ are presented in Fig. S1 of the SI.

In this study, we synthesized Ag₈TS₆ (T = Si, Ge, Sn) via a solid-state synthesis approach, and Rietveld refinements of their powder X-ray diffraction patterns at room temperature confirm the formation of single-phase materials (Fig. S2). Subsequently, temperature-dependent powder X-ray diffraction studies were conducted to investigate the presence of any phase transitions within the temperature range of 100 K to 400 K (Fig. S3). The Rietveld refinements of all diffraction patterns indicate that the orthorhombic phase with space group Pna2₁, remains stable for both Ag₈GeS₆ and Ag₈SiS₆ throughout the examined temperature range. Similarly, for Ag₈SnS₆, no clear change in the diffraction patterns was observed around 120 K. This may be because the structural variations are too small to detect the low-temperature structural change reported in the literature.¹⁶ Nevertheless, the refined unit cell volume of Ag₈SnS₆ below 120 K deviates slightly from a linear increase, suggesting that some structural change may occur at low temperature (Fig. S4). In contrast, a linear increase in unit cell volume for the other compositions was observed with increasing temperature (Fig. S4). We shortly note here that we predicted a potential additional phase of Ag₈SnS₆, so far not known from experimental work, via ab initio calculations (see Section S7 in the SI). This prediction might be an artifact of the density functional theory (DFT)-based methodology.

Bonding analysis

Based on the composition alone, we might naively expect Ag⁺, Si⁴⁺, Ge⁴⁺, Sn⁴⁺, and S²⁻ ions. A closer inspection of the structure, however, was already done by Krebs et al.,³⁷ and suggests the following ionic formula Ag₈(SiS₄)S₂ indicating SiS₄⁴⁻ polyanions isovalent to SiO₄⁴⁻. From previous bonding analysis results, we also expect very weakly bonded Ag atoms and highly covalent bonds from Si, Ge, and Sn to S.^{13,15,16,38,39}

Typically, the bonding situation in a material is used to estimate the sound velocities and to obtain information about the anharmonic nature of the heat transport.⁴⁰ For example, bond heterogeneity is typically made responsible for high phonon–phonon scattering rates and, consequently, low thermal conductivities.¹⁵ Specifically, in the case of Ag₈SnS₆, the rattler-like behavior of Ag⁺ is expected due to the very weak Ag–S bonds.¹⁵ Therefore, we provide a detailed analysis of the bonding situation in all three Ag₈TS₆ compounds (T = Si, Ge, Sn) by means of Crystal Orbital Hamilton Populations (COHP)⁴¹ and Crystal Orbital Bond Indices (COBI).⁴² The integration of the COHP up to the Fermi level (ICOHP) provides a quantitative measure of the bond strength. In the ICOHP analysis, the total electronic band energy is partitioned into pairwise atomic interactions, with more negative values indicating a stronger bond for a given atomic pair. Meanwhile, the integration of the COBI up to the Fermi level (ICOBI) corresponds to a bond order, serving as an indicator of covalency. The relationships between ICOHP, ICOBI, and bond lengths in inorganic materials have been explored as part of ref. 43. Beyond this, we also provide an analysis of metal–metal and multi-center interactions in these compounds, as they might be connected to the overall weak Ag–S bonds.¹⁵

The bonding situation in all three Ag₈TS₆ (T = Si, Ge, Sn) compounds is very similar. The T–S bonds are by far stronger and more covalent than the Ag–S bonds, indicated by both the ICOHP and ICOBI values (Fig. 2a and b). They also confirm the polyanionic nature of the TS₄⁴⁻ units, i.e., strong covalent bonds between T and S. The very covalent Sn–S bonds in Ag₈SnS₆ show an average ICOHP value of −4.58 eV and an average ICOBI value of 0.84 (close to the ideal ICOBI of 1 of a single bond). In contrast, the Ag–S interactions are much weaker, and the ICOHPs range from −0.66 to −1.61 eV (ICOBIs from 0.12 to 0.34). In the case of the COHPs, occupied antibonding states below the Fermi energy level weaken the Ag–S bonding interactions (see Fig. 2c). Specifically, Ag (4d) and S (3p) interactions contribute to the antibonding states. Likely due to weak Ag–S bonds, a large number of different, very distorted Ag⁺ environments exist. We found linear, trigonal planar, trigonal non-planar, and tetrahedral coordination environments for Ag⁺ in Ag₈SnS₆ (see Fig. S8 in SI). This analysis again confirms the expectation of the mobile nature of the Ag⁺ ions based on the bonding situation. In contrast, Sn only shows a nearly perfect tetrahedral environment. Besides cation–anion bonds, we also found Ag–Ag interactions in Ag₈SnS₆, with ICOHPs ranging from −0.24 to −0.32 eV (ICOBIs from 0.05 to 0.07), likely leading to additional distortions of the Ag environments and weakening of the Ag–S bonds. Ag–Ag interactions in metallic silver (fcc, Materials Project ID mp-124) have an ICOHP around −0.22 eV (ICOBI around 0.05), indicating the Ag–Ag interaction in Ag₈SnS₆ as a significant interaction. Typically, we consider ICOHPs larger −0.1 eV to be significant. The exact bond strengths and environments for all Ag₈TS₆ (T = Si, Ge, Sn) can be found in Fig. S6–S9 and Table S8–S11 in the SI.

Fig. 2 (a) and (b) show the distribution of ICOHP and ICOBI for the RT Ag₈SnS₆ structure, respectively. (c) and (d) depict the weakly bonded Ag–S COHP and COBI interactions at distinct Ag sites. Bonding interactions mainly involve Ag(5s/4d) and S(3p) orbitals, while the antibonding interactions below the Fermi level are dominated by Ag(4d) and S(3p) orbitals.

Plotting all two-center ICOBI^(2c) against each compound's bond length unveils an interesting pattern. The ICOBI vs. interatomic distance curve would fall monotonously in a regular compound without a unique bonding situation. Instead, we see unusually strong outliers for bond lengths beyond 3.5 Å (Fig. 2b). As previously shown in the literature,⁴⁴ these outliers indicate potential (hypervalent) multi-center interactions (a detailed discussion can be found in Section S2 in the SI). This is further investigated, i.e., the three-center (3c) bonds of consecutive atoms with stronger two-center ICOBI (ICOBI^(2c) ≥ 0.25) have been taken into account (see Section S2 in the SI for a detailed discussion). While we do find hypervalency/multicenter interactions (e.g., for S–Ag–S), they can be considered to be weak. Overall, it can be assumed that the weak Ag–S bonds, the Ag–Ag interactions, and the S–Ag–S multi-center interactions are closely related and therefore responsible for the anharmonicity of the compounds.

From the bonding analysis results, we see the overall bonding character stays the same when changing the tetrel species. As is known from simple binary compounds, bond strength and sound velocities are typically correlated.⁴⁵ Additionally, bond heterogeneity because of rattler-like atoms typically leads to high phonon–phonon scattering and anharmonicity. These similar results for all three compounds therefore suggest that all materials will present very similar sound velocities and high anharmonic transport behavior. Consequently, they are expected to exhibit similarly low lattice thermal conductivities and comparable features in their phonon band structures.

(Quasi-)harmonic phonon band structures

Checking the thermal stability of the thermoelectric materials is essential. Commonly, a lack of imaginary modes in the phonon band structure indicates dynamic stability of the structure. For all the argyrodites Ag₈TS₆ (T = Si, Ge, and Sn), the phonon frequencies along high-symmetry directions of their Brillouin zone and phonon density of states (PDOS) do not exhibit imaginary modes, which means that they are dynamically stable (Fig. 3 and S5). At Г, a significant splitting of the Longitudinal Optical-Transverse Optical (LO-TO) is visible, a direct result of the applied non-analytical term correction^46,47 that accounts for long-range dipole–dipole interactions. The low-temperature (LT) canfieldite Ag₈SnS₆ phase, using a 30-atom unit cell, has 90 phonon modes in total; while the Pna2₁ phases of Ag₈TS₆ (T = Si, Ge, and Sn), with 60 atoms per unit cell, have 180 modes. The phonon dispersion curve also shows considerable overlap between bands, indicating a possible high anharmonicity and a possible diffuson-dominated thermal transport.^12,17,22 Here, the PDOSs also show that the Si/Ge/Sn atoms make a small contribution across the entire region, while the S atoms mainly dominate the optical frequencies. The acoustic modes produce a dominant peak in the frequency range of 1.6 and 1.8 THz, which corresponds, due to their quantity and atomic mass, to the Ag⁺ vibrations. Similarly, the sound and group velocities computed by harmonic phonon and elastic tensor calculations are very similar between all three compounds (see a more detailed discussion in the SI Section S6). They also compare well with results from experimental sound-velocity measurements. Overall, no significant difference was found for the three compounds Ag₈TS₆ (T = Si, Ge, Sn) sharing the same crystal structure type.

Fig. 3 Computed phonon band structure along the partial phonon density of states for the room-temperature phase of the Ag₈SnS₆ canfieldite. Here, the dotted red line at 1.62 THz corresponds to the acoustic Debye frequency (ω_AC).

Grüneisen parameter

In general, the lattice thermal conductivity in a solid depends mainly on the heat capacity, speed of sound, and phonon relaxation time. Materials with low heat capacity, low group velocity, and short phonon lifetime have low lattice thermal conductivity. Both group velocity and phonon lifetime may depend on the bonding situation in the crystal. So far, we have found that Ag₈SnS₆, Ag₈GeS₆, and Ag₈SiS₆ all have weak Ag–S bonds and associated low sound velocities corresponding to Ag⁺ vibrations. Furthermore, we also expect high anharmonicity of Ag⁺ vibrations from the bonding analysis.

To quantify and evaluate anharmonicity as a function of phase and composition, we also compute the variation of phonon frequencies with respect to volume change as mode-dependent Grüneisen parameters and derive average quantities (γ). Given our previous results, we expected larger Grüneisen parameters for all three compounds, but no considerable differences between them. Fig. 4 shows strong anharmonicity represented by a large Grüneisen parameter for the low-energy vibrational modes (highlighted with grey), which are mainly dominated by Ag⁺ ions. This agrees with the expected mobile/rattler-like nature of the Ag⁺ ions. The averaged Grüneisen parameter was computed across all modes, showing good agreement with our experimental Grüneisen parameter derived from sound velocity measurements and the one reported previously in the literature (Ag₈GeS₆).¹³ Despite the change in composition, no significant differences were observed among the experimental average Grüneisen parameters derived from the experimental sound velocities.

Fig. 4 Computed mode Grüneisen parameter as a function of frequency for Ag₈SiS₆, Ag₈GeS₆, and Ag₈SnS₆ at room- and low-temperature. Here, we highlight the acoustic modes (grey color) where the anharmonicity is larger. Computed average Grüneisen parameter (γ) and experimental Grüneisen (γ_exp) are also shown for all investigated structures.

Although we observed comparable experimental and theoretical average Grüneisen values among the compounds, some differences are evident in our theoretical results at lower frequency modes. For instance, Ag₈SiS₆ mostly shows negative Grüneisen parameters for the lower frequencies. The calculation of the average Grüneisen parameter, shown in Fig. 4, was performed over all modes. Nevertheless, the average Grüneisen parameters used in the lattice thermal conductivity calculation were calculated with the acoustic modes only, as we expect them to be most important for thermal transport. A comparison of Grüneisen parameter computed over all modes, acoustic modes and up to the Debye frequency are presented in Fig. S20c in the SI.

Lattice thermal conductivity

Various models to predict lattice thermal conductivity have been developed. These range from simple empirical relationships to complex quantum mechanical calculations. These models vary strongly in required computational resources and also in how they model the heat transport in complex solids – either via phonons or diffusons.

Cahill²¹ and Agne²² have developed two alternative models that can be used cost-effectively with ab initio data. When combined with elastic properties obtained from DFT calculations, these models predict minimum lattice thermal conductivity. In both models, the amorphous solid has been used as a model system for the minimum thermal conductivity of crystalline materials; they both rely on random-walk theory, indicating heat transport in amorphous materials via diffusons. Because of this, these models can only be used in the high-temperature limit of crystalline materials.

Slack,^18–20 on the other hand, provides a lattice thermal conductivity model based on heat transport via phonons and as a function of temperature. This model emphasizes the role of acoustic phonon modes in thermal transport processes. According to Slack, the lattice thermal conductivity is influenced by factors such as the Debye temperature, sound velocity, and the Grüneisen parameter, which accounts for the anharmonicity of the lattice vibrations. The model is particularly useful for estimating the upper limit of thermal conductivity in crystalline materials with strong atomic bonding. In this approach, the lattice thermal conductivity can be computed as:

where M̄ is the average atomic mass, δ is the average volume per atom, n is the number of atoms per unit cell, Θ is the acoustic Debye temperature, T is the absolute temperature. A is the Slack coefficient, which is dependent on the anharmonicity of the structure, represented by the average Grüneisen parameter that we have already discussed above:

The Slack model is valuable as it provides a more thorough temperature-dependent analysis of thermal conductivity, offering insights that other models may not capture, especially in materials where acoustic phonons play a dominant role. While it yields important information about lattice thermal conductivity, the model also has limitations. Many studies point out that lattice thermal conductivity is generally overestimated when compared with experimental data. This discrepancy can be related to the A coefficient. Qin and coworkers⁴⁸ address this problem by scaling the A coefficient or by fitting the A parameter.

A recent model for high-throughput screening and analysis of thermal conductivity was introduced by Xia et al.,³⁴ where the lattice thermal conductivity can be estimated through harmonic phonon calculations. This model provides a complementary perspective to the previous methods since it builds upon the so-called two-channel model (phonon-gas channel and diffuson channel) where the total thermal conductivity κ_l is calculated from the sum of the phonon (s = s′) and diffuson contributions (s ≠ s′) (eqn (3)). s and s′ refer to the phonon branches.

C_ss′(q) is a heat capacity matrix element, v_ss′(q) is a velocity matrix element and τ_ss′(q) is a phonon lifetime matrix element for two phonons at the branches s and s′ in reciprocal space at q. τ_ss′(q) can be computed based on Γ_s(q) – the scattering rate or the inverse of phonon lifetime of the phonon at branch s at point q – and ω_s(q), its frequency:

This two-channel model, introduced by Simoncelli et al.,^35,49 has been very useful when disordered materials or crystals with large unit cells, such as Yb₁₄Mn₁Sb₁₁, have been investigated.⁵⁰ Simoncelli's model relies on the ab initio computation of phonon lifetimes, which can be computationally very demanding, making it extremely expensive for a large-scale screening approach. In contrast to this, the model by Xia³⁴ purely relies on harmonic phonon calculations, making it significantly more affordable. It additionally assumes that each phonon lifetime (1/Γ_s(q)) is half of its vibration period.

A comparison of lattice thermal conductivity using the models mentioned above is shown in Table S16. All models consistently predict the material's low lattice thermal conductivity, which can be attributed to its weak bonding, low sound velocities, and the high anharmonicity of the low-energy vibrational modes dominated by Ag⁺ ions. Furthermore, the diffusion-mediated minimum conductivity (κ^min_Agne) and two-channel (κ^min_Xia), both indicate that heat conduction is primarily dominated by diffusons, as it was also shown in previous studies including sulfide- and selenide-argyrodites.^12,17

Although, it is clear that all these models predict low minimal lattice thermal conductivities for Ag₈SnS₆, Ag₈GeS₆, and Ag₈SiS₆, a full ab initio model that can provide a detailed insight into the thermal properties is missing. Thus, to incorporate the anharmonicity in the prediction of the lattice thermal conductivity and to reduce the overestimation that the Slack model often shows, we start from the two-channel approach proposed by Xia et al. but go beyond the minimum lattice thermal conductivity approximation. We model the phonon lifetimes (1/Γ_s(q)) using the method proposed by Bjerg and co-workers, which is based on Slack's approach.^18,36 By incorporating inverse phonon lifetimes through the Grüneisen parameter, we effectively account for phonon–phonon scattering, which constitutes the dominant process limiting the lattice thermal conductivity in these materials. Then, the inverse phonon lifetimes are calculated as follows:

ω_s(q) is the phonon frequency, and p is a fitting function that is dependent on the average Grüneisen parameter γ, which can be determined by:Here v is the speed of sound and is determined from the Debye frequency ω_D, the number of atoms, and the volume of the cell, via the following equation:

Following the proposed model in this study, we computed the two-channel temperature-dependent lattice thermal conductivity, where the diagonal components of the heat flux matrix correspond to the phonon contribution, the off-diagonal components correspond to the diffuson contribution, and the total lattice thermal conductivity is obtained by summing both contributions. With this, in Fig. 5b and S21, we show the ultra-low total lattice thermal conductivity for the sulfide-argyrodite materials, with a very good agreement with the experimental measurements in the high temperature range.

Fig. 5 (a) Fit of the low-temperature measured thermal conductivity data (Experimental**, PPMS) of Ag₈GeS₆. High-temperature measurement using a second method (LFA, this study Experimental*), together with a measurement from literature¹⁷ are also provided in the plot for comparison. The fit is performed with the help of the analytical model as proposed in ref. 14 (scattering coefficients are presented in Table S19 in the SI. (b) Comparison of the lattice thermal conductivity following our proposed Grüneisen Model (GM) with experimental measurements for Ag₈TS₆ (T = Si, Ge, Sn). (c) Two-channel lattice thermal conductivity using the foundation model MACE-MP-03b. The model was not finetuned here. (d) Contribution of scattering process, phonon–phonon scattering (C₁), point-defect (C₂), and boundary scattering (A) on the phonon channel as obtained from the analytical model, compared with the Grüneisen model and the foundation model MACE-MP-03b (not finetuned). Although the two proposed approaches show minor differences, they remain consistent with the experimental results within a three-fold standard deviation, showing especially strong agreement for temperatures above 200 K.

For Ag₈GeS₆, where additional low-temperature experimental data are available, deviations are observed between 0 and 50 K compared to other measurements. We attribute this deviation to the presence of point-defect scattering, which can be caused when imperfections, such as atomic-scale substitutions, vacancies, or interstitials, disrupt the periodicity of the crystal lattice. This disruption could create a barrier to phonon propagation, significantly reducing lattice thermal conductivity.^12,14,51,52 Furthermore, the presence of microstructure features, such as grain boundaries, phase segregation, as well as different grain sizes in the experimental samples, can scatter phonons and decrease thermal conductivities and contribute to discrepancies with computational approaches that do not correct for these effects.^14,53

To estimate the influence of the point defects and the microstructure, we fitted the analytical model described in ref. 14 to the experimental data of Ag₈GeS₆. This fitting was previously presented and discussed in our earlier work,¹³ however, it is also included here to provide a complete comparison between our proposed models. This analytical model also accounts for phonon and diffuson channels and starts from harmonic phonon data as computed by DFT. Below the frequencies of the Ioffe-Regel limit, the Callaway model is used to describe the heat transport, while above this limit, the model introduced by Agne is used to describe the diffuson channel. However, estimations of lifetimes within the Callaway model, including effects from point defects and microstructure in the phonon lifetimes, are now fitted to the experimental data. Fig. 5a shows this fitting, demonstrating that heat transport can be accurately described based on this analytical model. It also indicates that the suppression of the phonon peak is predominantly driven by point-defect scattering and boundary scattering from microstructural characteristics, such as grain size. This analysis can be seen in Fig. 5d, which highlights the role of point-defect and boundary scattering within the phonon channel. For the Grüneisen-based approach, we observe a divergence in the 0–20 K range. However, the overall features of the temperature-dependent thermal conductivity agree very well with the experiment and especially our analytical model in which effects from point-defect and boundary scattering have been subtracted.

Given the experimental uncertainty, the foundation machine-learned interatomic potential (MACE-MP-03b, medium) reached good results in comparison with experiments as well, as shown in Fig. 5c and d. We did not fine-tune the model for the presented results. To compute the thermal conductivity, we use the full two-channel lattice dynamics approach implemented by Simoncelli et al.^35,49 Fig. 5c and S22, respectively, compare the phonon-channel and total lattice thermal conductivity obtained from the ML model and the Grüneisen parameter-based estimation. The result from the ML model agrees very well with the analytical model over the whole temperature range when point-defect and boundary scattering are subtracted. This again highlights the importance of point defects and boundary scattering for an accurate description of the thermal conductivity. Overall, the ML potential yields results consistent with the Grüneisen model, demonstrating that both approaches reliably capture this system's thermal transport behaviour – even in the low-temperature region.

Overall, following our proposed models, the results align with the findings of Ouyang and coworkers,¹⁷ confirming that heat transport in the argyrodites (Ag₈GeS₆ and Ag₈SnS₆ (RT)) is dominated by the diffuson-channel. We note that we neglected the influence of four-phonon scattering processes and additional temperature renormalizations of the harmonic phonons that slightly influence the results, in contrast to the simulations by Ouyang and coworkers.¹⁷ Theoretical predictions and experimental results reveal no significant differences among the three compositions, Ag₈SiS₆, Ag₈GeS₆, and Ag₈SnS₆. The Grüneisen-based model, however, results in a slight difference between the thermal conductivity of the Ag₈SiS₆ and Ag₈SnS₆ compounds, which also corresponds to the differences observed in the computed Grüneisen parameters.

As the Grüneisen-based model is computationally comparably cheap and a foundation MLIP model even requires less computational cost, they would both be suited for a high-throughput approach for screening thermal conductivity. However, it is currently unclear for which composition spaces foundational ML potentials might fail and how cheap finetuning for complex systems could look like. First finetuning tests with additional ab initio data from rattled supercells with an average displacement of 0.1 Å worsened the description of the phonon channel in our case, while the harmonic phonon results improved. We hope that automated MLIP training and finetuning capabilities will support establishing efficient training and finetuning procedures.^54,55 Despite these challenges, MLIPs are very promising as they allow for a full ab initio calculation of the lifetimes and include temperature renormalization effects in the phonons or four-phonon processes comparatively easily. Both the Grüneisen and MLIP approaches could also be combined to spot systematic failures of the foundation model within a high-throughput approach, or the Grüneisen model might be used together with a foundation MLIP. For heat capacity simulations, we have previously seen that even ML models with comparably poor predictions can yield good heat capacity estimates when sufficiently constrained by a physical model.⁵⁶

As in the previous studies on the argyrodite such as Ag₈GeSe₆, Cu₇PSe₆, and Ag_8–xCu_xGeS₆,^8,12–14 the thermal conductivity can also be modelled without considering changes in ionic conductivity. In those studies, the low thermal conductivity and the high ionic conductivities were shown to be independent of each other, as the ionic conductivities vary drastically and the thermal conductivities stay nearly constant within the same temperature range. This suggests that ionic conductivity does not directly control thermal transport; therefore, the thermal conductivity was modelled purely based on the lattice dynamics simulations. We furthermore experimentally investigate ionic transport properties to shed further light on the situation in Ag₈TS₆ (T = Si, Ge, Sn). We find that across all compositions, the ionic conductivity increases from an average of ∼0.003 mS cm⁻¹ at 233 K to ∼0.1 mS cm⁻¹ at 303 K, indicating an enhancement of more than one order of magnitude with rising temperature. In contrast, the experimental thermal conductivity remains nearly constant over the same temperature range, varying only marginally between ≈0.27 W mK⁻¹ and ≈0.28 W mK⁻¹. A more detailed discussion can be found in the SI Section S5. This observation suggests that ion transport has no direct influence on the observed low thermal conductivity in these materials, corroborating the previously reported findings for Ag⁺ and Cu⁺-based selenide and sulfur argyrodites^.8,12,13

Conclusions

Our results demonstrate a strong relationship between chemical bonding and lattice thermal conductivity in Ag-based sulfide argyrodites. The similar bonding strengths in all compounds lead to very similar sound velocities, while the weakly bonded Ag⁺ atoms result in high anharmonicity of vibrations, associated with high Grüneisen parameters. This weakness likely originates from occupied antibonding states in Ag–S bonds, Ag–Ag bonds, and the multi-center interactions. This further supports that the bonding situation might be predictive of a compound's overall thermal conductivity.

By applying both the Grüneisen-based and the MLIP-based models, we achieve good agreement with the experimental thermal conductivity data, especially in the medium-to high-temperature range. To capture the characteristic low-temperature peak, it is essential to include point-defect scattering, which effectively suppresses the phonon peak. In addition, microstructural features, most notably grain boundaries, introduce further boundary scattering, with grain size emerging as a key design parameter for tailoring thermal transport. Both effects are shown based on a fit of experimental data with an analytical model.

Overall, these results again demonstrate that accurately modelling heat transport in structurally complex materials over a large temperature range requires capturing the combined influence of bonding-driven anharmonicity, sound velocity, point-defect scattering, and microstructural effects. Furthermore, we identify two approaches that might be suitable for comparably cheap high-throughput screening of lattice thermal conductivity over wide temperature ranges.

Experimental and theoretical work

Methodology

Atomistic simulations. Electronic-structure computations were performed using density functional theory (DFT) as implemented in the Vienna Ab initio Simulation Package (VASP).^57–59 The exchange-correlation functional was treated in the semi-local approximation of Perdew, Burke, and Ernzerhof (PBE) with generalized gradient approximation (GGA).^60,61 The plane wave cut-off was set to 520 eV. To sample the Brillouin zone, we employed Γ-centre grid with a maximum separation of 0.12 Å⁻¹, which corresponds to a 7 × 7 × 5 and 3 × 7 × 5 k-points mesh for the orthorhombic (Pmn2₁) and orthorhombic (Pna2₁) structures, respectively. Structure optimizations were carried out in terms of volume, cell shape, and ionic positions.

The vibrational properties were computed using the supercell approach with the finite displacement method implemented in phonopy with displacements of 0.01 Å.^62,63 To obtain the dynamical matrix D(q), we used a supercell model of (3 × 3 × 2) and (1 × 3 × 2) for the LT and RT structures, respectively. The supercell calculations for the LT structure were performed at the Γ-point, while for the RT structure a 3 × 2 × 2 Γ-centred k-point grid was needed. In order to correct the dipole interaction, we also employed non-analytical term correction using Born charges as computed with VASP.^46,47

To compute volume-dependent thermal properties, we employed the Quasi-Harmonic Approximation (QHA),⁶⁴ implemented in phonopy.^29,65 To do so, we applied the harmonic approximation at expanded and contracted volumes. We start with the fully optimized structure at the ground state (V₀), and then we compute the constant volume energy of 13 different volumes from 0.943 × V₀ to 1.063 × V₀ in steps of 0.013 × V₀. The lattice parameters and atom positions were optimized by minimizing the electronic energy (ISIF = 4).⁶⁶ Additionally, to compute the anharmonicity of the structures, we compute the Grüneisen parameter. Here, two additional structural optimizations were performed at constant volume, 1% × V₀, and −1% × V₀.

To obtain Cahill's minimum thermal conductivity, we performed elastic constant calculation using an automated workflow implemented in atomate2, where elastic tensors are computed from stress-strain relationships.^67–69 More details can be found in the SI Section S6.

To compute the lattice thermal conductivity based on the foundation ML potential (here MACE-MP-03b medium model⁷⁰) together with the two-channel model introduced by Simoncelli et al.,^35,49 we solved the Wigner transport equation model as implemented in phono3py.^29,65 For this purpose, the third-order force constants were obtained with a supercell of 1 × 2 × 2, and the reciprocal space was sampled with a 6 × 14 × 10 mesh. Due to very demanding memory requirements, we only used the relaxation time approximation.

To get chemical insight into these compounds, bonding analysis was performed. To do so, we used our recently developed automatic bonding analysis workflow.⁷¹ The fully optimized structure for phonon computations is used as the input structure to start this workflow. The workflow then performs the bonding analysis with the LOBSTER^72–75 program by adding all necessary computational steps to the pipeline. This pipeline consists of a static DFT computation using the GGA functional parameterized by PBE^60,61 within the PAW framework^.76,77 A grid density of 6000 k-points per reciprocal atom is set for the DFT run. The electronic structure's convergence criterion and the plane-wave energy cutoff are set to 10⁻⁶ and 520 eV, respectively. The number of grid points (NEDOS) on which the density of states is evaluated is set to 10 000. The Brillouin zone is integrated using the tetrahedron method with Blöchl⁷⁸ correction (i.e., ISMEAR = −5). In all DFT computations, spin polarization is switched on, even though this is not required for these compounds. The workflow also performs LOBSTER computations with the available basis for projecting the wavefunctions. Here, we report the results on the minimal basis.

For bonding analysis runs via LOBSTER, COHPs and COBIs are computed for the entire energy range of VASP static runs, and the COHP/COBI energy interval step is set to 10 000 points (equal to NEDOS set in the VASP static run). The increased number of points assigned for the COHP/COBI computation poses a very good estimate of bonding and anti-bonding contribution in bonds during post-processing the results via LobsterPy. Three-center interactions to calculate three-center COBI and ICOBI were chosen according to stronger two-center ICOBIs (cutoff ICOBI⁽²⁾ = 0.2) of three consecutive atoms and automatically analyzed using a new implementation by one of the current authors in pymatgen (as of v2023.10.11).⁷⁹ Other multi-center bonds have been checked as well, but did not yield significant values (cutoff ICOBI(n) = ±0.05).

Solid-state synthesis of Ag₈TS₆ (T = Si, Ge, Sn). The synthesis of Ag₈TS₆ (T = Si, Ge, Sn) utilized reactants including silver powder (99.9%, sigma aldrich), silicon (99.999%, Thermo Scientific)), germanium (99.999%, sigma aldrich), tin (99.85%, Thermo Scientific), and sulfur powder (99.98%, sigma aldrich). A high-temperature solid-state synthesis method was conducted under static vacuum conditions to produce polycrystalline samples of Ag₈TS₆. Initially, stoichiometric amounts of the reactants were weighed inside an argon-filled glovebox and placed into carbon-coated quartz ampoules, which had been pre-dried at 1073 K for 2 hours under dynamic vacuum. These ampoules were then sealed under vacuum and heated in a horizontal tube furnace. The heating process involved ramping the temperature to 523 K at a rate of 50 K per hour, followed by a 24-hour hold. Subsequently, the temperature was increased to 1250 K at the same rate, maintained for 60 hours, and then cooled down to room temperature.

X-ray diffraction. X-ray diffraction patterns of Ag₈TS₆ (T = Si, Ge, Sn) were collected using a STOE STADIP diffractometer. The setup utilized Mo Kα1 radiation (λ = 0.7093 Å) equipped with curved Ge (111) monochromator and a Mythen2 1 K detector. Measurements were performed in the Debye–Scherrer geometry over a 2θ range from 4° to 44°, at a scan rate of 1° per minute. The temperature range during these measurements was between 100 K and 400 K, maintained using a Cryostream 1000 cooler from Oxford Cryosystems Ltd for low-temperature conditions (<298 K). Samples were prepared in borosilicate glass capillaries with a 0.5 mm diameter, and they were equilibrated for 20 minutes at each temperature step prior to the measurement. Details of structural phase analysis and Rietveld refinements are discussed in Section S1 in the SI.

Ultrasonic speed of sound measurement. An Olympus Epoch 600 with 5 MHz transducers was employed to measure speed of sound using the pulse-echo method. Variations in signal measurements and the geometrical factors (such as density and thickness) were accounted for to determine the uncertainty of the speed of sound measurement.

Thermal transport properties measurement. A Netzsch LFA-467 instrument was used to measure thermal diffusivity of all three compositions, using 10 mm diameter, disc-shaped samples with a bulk density of approximately more than 95% of the theoretical density. Measurements were conducted over a temperature range of 173 K to 500 K. An MCT detector with a ZnS furnace window was used for the measurements below room temperature; while for measurements from room temperature to high-temperature, an InSb detector with a sapphire furnace window was employed. The detection time and signal amplification were optimized automatically for each measurement. At every temperature point, three measurements were taken, with five measurements conducted at 173 K to ensure accuracy. The detector signal was analyzed using an improved Cape–Lehman model.^80–82 All samples were spray-coated with graphite to enhance the infrared light absorption and emission during the laser-flash measurements. The equations used to calculate the thermal conductivity from the measured thermal diffusivity are provided in Section S3 in the SI. Additionally, low-temperature thermal conductivity measurement for Ag₈GeS₆ was performed using a Physical Property Measurement System (PPMS) with the TTO option, under high vacuum (<10⁻⁴ Torr) and with a temperature gradient of approximately 3% between the hot and cold sides. A disc-shaped sample (4 mm × 2 mm) was used for the measurement.

Direct current (DC) polarization measurements. DC polarization measurements were performed using a press cell with a 10 mm inner diameter. The samples were filled into the press cell, and stainless-steel stamps were used as ion blocking electrode on both sides. The cells were then closed and subjected to uniaxial pressing at 3 tons for 3 minutes. A VMP-300 potentiostat (Biologic) was used to carry out DC polarization, applying a voltage ranging from 5 mV to 50 mV in 5 mV steps. To ensure equilibrium at each step, the applied voltage was kept constant for 2 hours before proceeding to the next step.

Electrochemical impedance spectroscopy. Electrochemical impedance spectroscopy measurements were carried out using a cell set up comprising two stainless steel stamps that served both as current collectors and as a means to press the sample during the measurements. The samples were placed in an insulating PEEK housing with an inner diameter of 10 mm. First, the argyrodite materials were loaded in the PEEK housing and pressed under 3 tons of uniaxial pressure for 3 minutes. Secondly, the cells were opened in a glovebox, and a thin layer of RbAg₄I₅ was pressed onto both sides of the sample, followed by an additional 5 minutes of compression using a manual screw press. Finally, AC impedance spectroscopy was performed over the temperature range of 233–303 K using an SP300 impedance analyzer (Biologic). The measurements employed an excitation amplitude of 10 mV and covered a frequency range of 5 MHz to 1 Hz. The analysis of impedance results is shown in Section S5 in the SI.

Author contributions

Conceptualization: JB, JG; data curation: JB, AN, CE, KU, AG, JG; formal analysis: JB, JG, AN, CE, KU, AG, WZ; funding acquisition: WZ, JG; investigation: JB, AN, CE, AG, WZ, JG; methodology: JB, JG, AN, CE, KU, AG, WZ; project administration: WZ, JG; resources: WZ, JG; software: JB, JG, AN, CE, KU, AG; supervision: WZ, JG; validation: JB, JG; visualization: JB, AN, CE, KU, AG; writing – original draft: JB, JG, AG, WZ with contributions from all authors; writing – review & editing: All authors.

Conflicts of interest

There are no conflicts to declare.

Data availability

The raw computational data required to reproduce the results of this study are available through our ZENODO repositories https://doi.org/10.5281/zenodo.17397457 and https://doi.org/10.5281/zenodo.17399976. Post-processing and plotting scripts used to generate figures presented in this work are openly accessible on our GitHub repository https://github.com/DigiMatChem/paper-grueneisen-model-for-argyrodites/(https://doi.org/10.5281/zenodo.19017707). Detailed instructions for reproducing each dataset and figure are provided in the accompanying README file.

Supplementary information (SI): further discussion of the computational and experimental results is available in the SI. The SI also contains citations to ref. 83–90. See DOI: https://doi.org/10.1039/d5ta08709k.

Acknowledgements

The research was supported by ERC Grant MultiBonds (grant agreement no. 101161771; funded by the European Union. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.) The research was additionally supported by the Deutsche Forschungsgemeinschaft (DFG) under grant number ZE 1010/15-1 and project number 459785385. Furthermore, we would like to acknowledge the Gauss Centre for Supercomputing e.V. (https://www.gausscentre.eu) for funding this project by providing generous computing time on the GCS Supercomputer SuperMUC-NG at Leibniz Supercomputing Centre (https://www.lrz.de) (Project No. pn73da).

Notes and references

1. S. Lin, W. Li and Y. Pei, Thermally insulative thermoelectric argyrodites, Mater. Today, 2021, 48, 198–213.
2. W. Li, S. Lin, B. Ge, J. Yang, W. Zhang and Y. Pei, Low Sound Velocity Contributing to the High Thermoelectric Performance of Ag₈SnSe₆, Adv. Sci., 2016, 3, 1600196.
3. X. Shen, High-Temperature Structural and Thermoelectric Study of Argyrodite Ag₈GeSe₈, ACS Appl. Mater. Interfaces, 2018, 11, 2168–2176.
4. Q. Ren, M. K. Gupta, M. Jin, J. Ding, J. Wu, Z. Chen, S. Lin, O. Fabelo, J. A. Rodríguez-Velamazán, M. Kofu, K. Nakajima, M. Wolf, F. Zhu, J. Wang, Z. Cheng, G. Wang, X. Tong, Y. Pei, O. Delaire and J. Ma, Extreme phonon anharmonicity underpins superionic diffusion and ultralow thermal conductivity in argyrodite Ag₈SnSe₆, Nat. Mater., 2023, 1–8.
5. M. K. Gupta, J. Ding, D. Bansal, D. L. Abernathy, G. Ehlers, N. C. Osti, W. G. Zeier and O. Delaire, Strongly Anharmonic Phonons and Their Role in Superionic Diffusion and Ultralow Thermal Conductivity of Cu₇PSe₆, Adv. Energy Mater., 2022, 12, 2200596.
6. K. S. Weldert, W. G. Zeier, T. W. Day, M. Panthöfer, G. J. Snyder and W. Tremel, Thermoelectric Transport in Cu7PSe6 with High Copper Ionic Mobility, J. Am. Chem. Soc., 2014, 136, 12035–12040.
7. B. K. Heep, K. S. Weldert, Y. Krysiak, T. W. Day, W. G. Zeier, U. Kolb, G. J. Snyder and W. Tremel, High Electron Mobility and Disorder Induced by Silver Ion Migration Lead to Good Thermoelectric Performance in the Argyrodite Ag₈SiSe₆, Chem. Mater., 2017, 29, 4833–4839.
8. A. Ghata, T. Bernges, O. Maus, B. Wankmiller, A. A. Naik, J. Bustamante, M. W. Gaultois, O. Delaire, M. R. Hansen, J. George and W. G. Zeier, Exploring the Thermal and Ionic Transport of Cu+ Conducting Argyrodite Cu₇PSe₆, Adv. Energy Mater., 2024, 14, 2402039.
9. W. F. Kuhs, R. Nitsche and K. Scheunemann, The argyrodites — A new family of tetrahedrally close-packed structures, Mater. Res. Bull., 1979, 14, 241–248.
10. S. Lin, W. Li and Y. Pei, Thermally insulative thermoelectric argyrodites, Mater. Today, 2021, 48, 198–213.
11. J. Dong, L. Hu, J. Liu, Y. Liu, Y. Jiang, Z. Yu, X. Y. Tan, A. Suwardi, Q. Zheng, Q. Li, J.-F. Li, V. P. Dravid, Q. Yan and M. G. Kanatzidis, Off-Centering of Ge Atoms in GeBi₂Te₄ and Impact on Thermoelectric Performance, Adv. Funct. Mater., 2024, 34, 2314499.
12. T. Bernges, R. Hanus, B. Wankmiller, K. Imasato, S. Lin, M. Ghidiu, M. Gerlitz, M. Peterlechner, S. Graham, G. Hautier, Y. Pei, M. R. Hansen, G. Wilde, G. J. Snyder, J. George, M. T. Agne and W. G. Zeier, Considering the Role of Ion Transport in Diffuson-Dominated Thermal Conductivity, Adv. Energy Mater., 2022, 12, 2200717.
13. A. Ghata, P. L. Eckert, T. Böger, P. Garg and W. G. Zeier, Influence of Cu+ Substitution on the Structural, Ionic, and Thermal Transport Properties of Ag_8–xCu_xGeS₆ Argyrodites, Chem. Mater., 2025, 37, 6900–6911.
14. T. Bernges, M. Peterlechner, G. Wilde, M. T. Agne and W. G. Zeier, Analytical model for two-channel phonon transport engineering, Mater. Today Phys., 2023, 35, 101107.
15. X. Shen, Y. Xia, C. Yang, Z. Zhang, S. Li, Y.-H. Tung, A. Benton, X. Zhang, X. Lu, G. Wang, J. He and X. Zhou, High Thermoelectric Performance in Sulfde-Type Argyrodites Compound Ag₈Sn(S_1-xSe_x)₆ Enabled by Ultralow Lattice Thermal Conductivity and Extended Cubic Phase Regime, Adv. Funct. Mater., 2020, 30, 2000526,  DOI:10.1002/adfm.202000526.
16. T. J. Slade, V. Gvozdetskyi, J. M. Wilde, A. Kreyssig, E. Gati, L.-L. Wang, Y. Mudryk, R. A. Ribeiro, V. K. Pecharsky, J. V. Zaikina, S. L. Bud’ko and P. C. Canfield, A Low-Temperature Structural Transition in Canfieldite, Ag₈SnS₆, Single Crystals, Inorg. Chem., 2021, 60, 19345–19355.
17. N. Ouyang, D. Shen, C. Wang, R. Cheng, Q. Wang and Y. Chen, Positive temperature-dependent thermal conductivity induced by wavelike phonons in complex Ag-based argyrodites, Phys. Rev. B, 2025, 111, 064307.
18. G. A. Slack and S. Galginaitis, Thermal Conductivity and Phonon Scattering by Magnetic Impurities in CdTe, Phys. Rev., 1964, 133, A253–A268.
19. G. A. Slack, in Solid State Physics, H. Ehrenreich, F. Seitz and D. Turnbull, Academic Press, 1979, pp. 1–71.
20. G. A. Slack, Nonmetallic crystals with high thermal conductivity, J. Phys. Chem. Solids, 1973, 34, 321–335.
21. D. G. Cahill, S. K. Watson and R. O. Pohl, Lower limit to the thermal conductivity of disordered crystals, Phys. Rev. B, 1992, 46, 6131–6140.
22. M. T. Agne, R. Hanus and G. J. Snyder, Minimum thermal conductivity in the context of diffuson-mediated thermal transport, Energy Environ. Sci., 2018, 11, 609–616.
23. R. Hanus, R. Gurunathan, L. Lindsay, M. T. Agne, J. Shi, S. Graham and G. Jeffrey Snyder, Thermal transport in defective and disordered materials, Appl. Phys. Rev., 2021, 8, 031311.
24. F. Eriksson, E. Fransson and P. Erhart, The Hiphive Package for the Extraction of High-Order Force Constants by Machine Learning, Adv. Theory Simul., 2019, 2, 1800184.
25. A. Rosado-Miranda, V. Posligua, J. Fdez Sanz, A. Márquez, P. Nath and J. Plata, Design Principles Guided by DFT Calculations and High-Throughput Frameworks for the Discovery of New Diamond-like Chalcogenide Thermoelectric Materials, ACS Appl. Mater. Interfaces, 2024, 16, 28590–28598.
26. Y. Luo, M. Li, H. Yuan, H. Liu and Y. Fang, Predicting lattice thermal conductivity via machine learning: a mini review, npj Comput. Mater., 2023, 9, 1–11.
27. A. Loew, D. Sun, H.-C. Wang, S. Botti and M. A. L. Marques, Universal machine learning interatomic potentials are ready for phonons, npj Comput. Mater., 2025, 11, 178.
28. B. Póta, P. Ahlawat, G. Csányi and M. Simoncelli, arXiv, 2025, preprint, arXiv:2408.00755,  DOI:10.48550/arXiv.2408.00755.
29. A. Togo, L. Chaput, T. Tadano and I. Tanaka, Implementation strategies in phonopy and phono3py, J. Phys.: Condens. Matter, 2023, 35, 353001.
30. A. Togo, L. Chaput and I. Tanaka, Distributions of phonon lifetimes in Brillouin zones, Phys. Rev. B, 2015, 91, 094306.
31. R. Kubo, The fluctuation-dissipation theorem, Rep. Prog. Phys., 1966, 29, 255.
32. M. S. Green, Markoff Random Processes and the Statistical Mechanics of Time-Dependent Phenomena. II. Irreversible Processes in Fluids, J. Chem. Phys., 1954, 22, 398–413.
33. L. Chaput, Direct Solution to the Linearized Phonon Boltzmann Equation, Phys. Rev. Lett., 2013, 110, 265506.
34. Y. Xia, D. Gaines, J. He, K. Pal, Z. Li, M. G. Kanatzidis, V. Ozoliņš and C. Wolverton, A unified understanding of minimum lattice thermal conductivity, Proc. Natl. Acad. Sci. U. S. A., 2023, 120, e2302541120.
35. M. Simoncelli, N. Marzari and F. Mauri, Unified theory of thermal transport in crystals and glasses, Nat. Phys., 2019, 15, 809–813.
36. L. Bjerg, B. B. Iversen and G. K. H. Madsen, Modeling the thermal conductivities of the zinc antimonides ZnSb and Zn 4 Sb 3, Phys. Rev. B, 2014, 89, 024304.
37. B. Krebs and J. Mandt, Zur Kenntnis des Argyrodit-Strukturtyps: Die Kristallstruktur von Ag₈SiS₆, Z. Naturforsch. B, 1977, 32, 373–379.
38. G. Eulenberger, Die Kristallstruktur der Tieftemperaturmodifikation von Ag₈GeS₆, Monatsh. Chem., 1977, 108, 901–913.
39. C. Sturm, N. Boccalon, D. Ramirez and H. Kleinke, Stability and Thermoelectric Properties of the Canfieldite Ag₈SnS₆, ACS Appl. Energy Mater., 2021, 4, 10244–10251.
40. W. G. Zeier, A. Zevalkink, Z. M. Gibbs, G. Hautier, M. G. Kanatzidis and G. J. Snyder, Thinking Like a Chemist: Intuition in Thermoelectric Materials, Angew. Chem., Int. Ed., 2016, 55, 6826.
41. R. Dronskowski and P. E. Blöchl, Crystal orbital Hamilton populations (COHP): energy-resolved visualization of chemical bonding in solids based on density-functional calculations, J. Phys. Chem., 1993, 97, 8617–8624.
42. P. C. Müller, C. Ertural, J. Hempelmann and R. Dronskowski, Crystal Orbital Bond Index: Covalent Bond Orders in Solids, J. Phys. Chem. C, 2021, 125, 7959–7970.
43. A. A. Naik, C. Ertural, N. Dhamrait, P. Benner and J. George, A Quantum-Chemical Bonding Database for Solid-State Materials, Sci. Data, 2023, 10, 610.
44. C. Ertural, Über die elektronische Struktur funktioneller Festkörpermaterialien und ihre Beschreibung mittels lokaler Bindungsindikatoren, RWTH Aachen University, 2022.
45. W. G. Zeier, A. Zevalkink, Z. M. Gibbs, G. Hautier, M. G. Kanatzidis and G. J. Snyder, Thinking Like a Chemist: Intuition in Thermoelectric Materials, Angew. Chem., Int. Ed., 2016, 55, 6826–6841.
46. X. Gonze, J.-C. Charlier, D. C. Allan and M. P. Teter, Interatomic force constants from first principles: The case of \ensuremath{\alpha}-quartz, Phys. Rev. B:Condens. Matter Mater. Phys., 1994, 50, 13035–13038.
47. X. Gonze and C. Lee, Dynamical matrices, Born effective charges, dielectric permittivity tensors, and interatomic force constants from density-functional perturbation theory, Phys. Rev. B:Condens. Matter Mater. Phys., 1997, 55, 10355–10368.
48. G. Qin, A. Huang, Y. Liu, H. Wang, Z. Qin, X. Jiang, J. Zhao, J. Hu and M. Hu, High-throughput computational evaluation of lattice thermal conductivity using an optimized Slack model, Mater. Adv., 2022, 3, 6826–6830.
49. M. Simoncelli, N. Marzari and F. Mauri, Wigner Formulation of Thermal Transport in Solids, Phys. Rev. X, 2022, 12, 041011.
50. R. Hanus, J. George, M. Wood, A. Bonkowski, Y. Cheng, D. L. Abernathy, M. E. Manley, G. Hautier, G. J. Snyder and R. P. Hermann, Uncovering design principles for amorphous-like heat conduction using two-channel lattice dynamics, Mater. Today Phys., 2021, 18, 100344.
51. Y. Sun, Y. Liu, R. Li, Y. Li and S. Bai, Strategies to Improve the Thermoelectric Figure of Merit in Thermoelectric Functional Materials, Front. Chem., 2022, 10, 865281.
52. R. Gurunathan, R. Hanus, M. Dylla, A. Katre and G. J. Snyder, Analytical Models of Phonon-Point-Defect Scattering, Phys. Rev. Appl., 2020, 13, 034011.
53. T. Böger, T. Bernges, M. T. Agne, P. Canepa, F. Tietz and W. G. Zeier, On the Thermal Conductivity and Local Lattice Dynamical Properties of NASICON Solid Electrolytes, J. Am. Chem. Soc., 2024, 146, 32678–32688.
54. Y. Liu, J. D. Morrow, C. Ertural, N. L. Fragapane, J. L. A. Gardner, A. A. Naik, Y. Zhou, J. George and V. L. Deringer, An automated framework for exploring and learning potential-energy surfaces, Nat. Commun., 2025, 16, 7666.
55. S. Menon, Y. Lysogorskiy, A. L. M. Knoll, N. Leimeroth, M. Poul, M. Qamar, J. Janssen, M. Mrovec, J. Rohrer, K. Albe, J. Behler, R. Drautz and J. Neugebauer, From electrons to phase diagrams with machine learning potentials using pyiron based automated workflows, npj Comput. Mater., 2024, 10, 261.
56. E. J. Cope, J. Bustamante, Z. M. Johnson, A. Lancaster, R. Gurunathan, J. George and M. T. Agne, Heat capacity estimation of complex materials for energy technologies, Joule, 2025.
57. G. Kresse and J. Hafner, Ab initio molecular dynamics for liquid metals, Phys. Rev. B:Condens. Matter Mater. Phys., 1993, 47, 558–561.
58. G. Kresse, J. Furthmüller and J. Hafner, Theory of the crystal structures of selenium and tellurium: The effect of generalized-gradient corrections to the local-density approximation, Phys. Rev. B:Condens. Matter Mater. Phys., 1994, 50, 13181–13185.
59. G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B:Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
60. J. P. Perdew, K. Burke and M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett., 1996, 77, 3865–3868.
61. J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh and C. Fiolhais, Atoms, molecules, solids, and surfaces: Applications of the generalized gradient approximation for exchange and correlation, Phys. Rev. B:Condens. Matter Mater. Phys., 1992, 46, 6671–6687.
62. A. Togo and I. Tanaka, First principles phonon calculations in materials science, Scr. Mater., 2015, 108, 1–5.
63. A. Togo, F. Oba and I. Tanaka, First-principles calculations of the ferroelastic transition between rutile-type and CaCl₂ -type SiO ₂ at high pressures, Phys. Rev. B:Condens. Matter Mater. Phys., 2008, 78, 134106.
64. R. P. Stoffel, C. Wessel, M.-W. Lumey and R. Dronskowski, Ab Initio Thermochemistry of Solid-State Materials, Angew. Chem., Int. Ed., 2010, 49, 5242–5266.
65. A. Togo, First-principles Phonon Calculations with Phonopy and Phono3py, J. Phys. Soc. Jpn., 2023, 92, 012001.
66. K. Latimer, S. Dwaraknath, K. Mathew, D. Winston and K. A. Persson, Evaluation of thermodynamic equations of state across chemistry and structure in the materials project, npj Comput. Mater., 2018, 4, 40.
67. M. de Jong, W. Chen, T. Angsten, A. Jain, R. Notestine, A. Gamst, M. Sluiter, C. Krishna Ande, S. van der Zwaag, J. J. Plata, C. Toher, S. Curtarolo, G. Ceder, K. A. Persson and M. Asta, Charting the complete elastic properties of inorganic crystalline compounds, Sci. Data, 2015, 2, 150009.
68. A. Togo, K. Shinohara and I. Tanaka, Spglib: a software library for crystal symmetry search, Sci. Technol. Adv. Mater., 2024, 4, 2384822.
69. A. M. Ganose, H. Sahasrabuddhe, M. Asta, K. Beck, T. Biswas, A. Bonkowski, J. Bustamante, X. Chen, Y. Chiang, D. C. Chrzan, J. Clary, O. A. Cohen, C. Ertural, M. C. Gallant, J. George, S. Gerits, R. E. A. Goodall, R. D. Guha, G. Hautier, M. Horton, T. J. Inizan, A. D. Kaplan, R. S. Kingsbury, M. C. Kuner, B. Li, X. Linn, M. J. McDermott, R. S. Mohanakrishnan, A. N. Naik, J. B. Neaton, S. M. Parmar, K. A. Persson, G. Petretto, T. A. R. Purcell, F. Ricci, B. Rich, J. Riebesell, G.-M. Rignanese, A. S. Rosen, M. Scheffler, J. Schmidt, J.-X. Shen, A. Sobolev, R. Sundararaman, C. Tezak, V. Trinquet, J. B. Varley, D. Vigil-Fowler, D. Wang, D. Waroquiers, M. Wen, H. Yang, H. Zheng, J. Zheng, Z. Zhu and A. Jain, Atomate2: modular workflows for materials science, Digital Discovery, 2025, 4, 1944–1973.
70. I. Batatia, P. Benner, Y. Chiang, A. M. Elena, D. P. Kovács, J. Riebesell, X. R. Advincula, M. Asta, M. Avaylon, W. J. Baldwin, F. Berger, N. Bernstein, A. Bhowmik, F. Bigi, S. M. Blau, V. Cărare, M. Ceriotti, S. Chong, J. P. Darby, S. De, F. Della Pia, V. L. Deringer, R. Elijošius, Z. El-Machachi, E. Fako, F. Falcioni, A. C. Ferrari, J. L. A. Gardner, M. J. Gawkowski, A. Genreith-Schriever, J. George, R. E. A. Goodall, J. Grandel, C. P. Grey, P. Grigorev, S. Han, W. Handley, H. H. Heenen, K. Hermansson, C. H. Ho, S. Hofmann, C. Holm, J. Jaafar, K. S. Jakob, H. Jung, V. Kapil, A. D. Kaplan, N. Karimitari, J. R. Kermode, P. Kourtis, N. Kroupa, J. Kullgren, M. C. Kuner, D. Kuryla, G. Liepuoniute, C. Lin, J. T. Margraf, I.-B. Magdău, A. Michaelides, J. H. Moore, A. A. Naik, S. P. Niblett, S. W. Norwood, N. O'Neill, C. Ortner, K. A. Persson, K. Reuter, A. S. Rosen, L. A. M. Rosset, L. L. Schaaf, C. Schran, B. X. Shi, E. Sivonxay, T. K. Stenczel, C. Sutton, V. Svahn, T. D. Swinburne, J. Tilly, C. van der Oord, S. Vargas, E. Varga-Umbrich, T. Vegge, M. Vondrák, Y. Wang, W. C. Witt, T. Wolf, F. Zills and G. Csányi, A foundation model for atomistic materials chemistry, J. Chem. Phys., 2025, 163, 184110.
71. J. George, G. Petretto, A. Naik, M. Esters, A. J. Jackson, R. Nelson, R. Dronskowski, G. Rignanese and G. Hautier, Automated Bonding Analysis with Crystal Orbital Hamilton Populations, ChemPlusChem, 2022, 87, e202200246.
72. V. L. Deringer, A. L. Tchougréeff and R. Dronskowski, Crystal Orbital Hamilton Population (COHP) Analysis As Projected from Plane-Wave Basis Sets, J. Phys. Chem. A, 2011, 115, 5461–5466.
73. S. Maintz, V. L. Deringer, A. L. Tchougréeff and R. Dronskowski, Analytic projection from plane-wave and PAW wavefunctions and application to chemical-bonding analysis in solids, J. Comput. Chem., 2013, 34, 2557–2567.
74. S. Maintz, V. L. Deringer, A. L. Tchougréeff and R. Dronskowski, LOBSTER: A tool to extract chemical bonding from plane-wave based DFT: Tool to Extract Chemical Bonding, J. Comput. Chem., 2016, 37, 1030–1035.
75. R. Nelson, C. Ertural, J. George, V. L. Deringer, G. Hautier and R. Dronskowski, LOBSTER : Local orbital projections, atomic charges, and chemical-bonding analysis from projector-augmented-wave-based density-functional theory, J. Comput. Chem., 2020, 41, 1931–1940.
76. G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B:Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
77. P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B:Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
78. P. E. Blöchl, O. Jepsen and O. K. Andersen, Improved tetrahedron method for Brillouin-zone integrations, Phys. Rev. B:Condens. Matter Mater. Phys., 1994, 49, 16223–16233.
79. S. P. Ong, W. D. Richards, A. Jain, G. Hautier, M. Kocher, S. Cholia, D. Gunter, V. L. Chevrier, K. A. Persson and G. Ceder, Python Materials Genomics (pymatgen): A robust, open-source python library for materials analysis, Comput. Mater. Sci., 2013, 68, 314–319.
80. J. A. Cape and G. W. Lehman, Temperature and Finite Pulse-Time Effects in the Flash Method for Measuring Thermal Diffusivity, J. Appl. Phys., 1963, 34, 1909–1913.
81. R. E. Taylor and J. A. Cape, Finite Pulse-Time Effects In The Flash Diffusivity Technique, Appl. Phys. Lett., 1964, 5, 212–213.
82. J. Blumm and J. Opfermann, Improvement of the mathematical modeling of flash measurements, High. Temp.–High. Press., 2002, 34, 515–521.
83. A. A. Coelho, TOPAS and TOPAS-Academic: an optimization program integrating computer algebra and crystallographic objects written in C++, J. Appl. Crystallogr., 2018, 51, 210–218.
84. P. Thompson, D. E. Cox and J. B. Hastings, Rietveld refinement of Debye–Scherrer synchrotron X-ray data from Al2O3, J. Appl. Crystallogr., 1987, 20, 79–83.
85. J. Hempelmann, P. C. Müller, C. Ertural and R. Dronskowski, The Orbital Origins of Chemical Bonding in Ge−Sb−Te Phase-Change Materials**, Angew. Chem., Int. Ed., 2022, 61, e202115778.
86. A. Yadav, P. Deshmukh, K. Roberts, N. Jisrawi and S. Valluri, An analytic study of the Wiedemann–Franz law and the thermoelectric figure of merit, J. Phys. Commun., 2019, 3, 105001.
87. K. Funke, R. D. Banhatti, D. Wilmer, R. Dinnebier, A. Fitch and M. Jansen, Low-Temperature Phases of Rubidium Silver Iodide: Crystal Structures and Dynamics of the Mobile Silver Ions, J. Phys. Chem. A, 2006, 110, 3010–3016.
88. A. I. Pogodin, M. J. Filep, V. Yu. Izai, O. P. Kokhan and P. Kúš, Crystal growth and electrical conductivity of Ag₇PS₆ and Ag₈GeS₆argyrodites, J. Phys. Chem. Solids, 2022, 168, 110828.
89. J. T. S. Irvine, D. C. Sinclair and A. R. West, Electroceramics: Characterization by Impedance Spectroscopy, Adv. Mater., 1990, 2, 132–138.
90. W. Chen, J.-H. Pöhls, G. Hautier, D. Broberg, S. Bajaj, U. Aydemir, Z. M. Gibbs, H. Zhu, M. Asta, G. J. Snyder, B. Meredig, M. A. White, K. Persson and A. Jain, Understanding thermoelectric properties from high-throughput calculations: trends, insights, and comparisons with experiment, J. Mater. Chem. C, 2016, 4, 4414–4426.

---