Accelerated Discovery of a Large Family of Quaternary Chalcogenides with very Low Lattice Thermal Conductivity

The development of efficient thermal energy management devices such as thermoelectrics, barrier coatings, and thermal data-storage disks often relies on compounds that possess very low lattice thermal conductivity ($\kappa_l$). Here, we present the computational prediction of a large family of 628 thermodynamically stable quaternary chalcogenides, AMM'Q$_3$ (A = alkali/alkaline earth/post-transition metals; M/M' = transition metals, lanthanides; Q = chalcogens) using high-throughput density functional theory (DFT) calculations. We validate the presence of low-$\kappa_l$ in this family of materials by calculating $\kappa_l$ of several predicted stable compounds using the Peierls-Boltzmann transport equation within a first-principles DFT framework. Our analysis reveals that the low-$\kappa_l$ in the AMM'Q$_3$ compounds originates from the presence of either a strong lattice anharmonicity that enhances the phonon scatterings or rattlers cations that lead to multiple scattering channels in their crystal structures. Our predictions suggest new experimental research opportunities in the synthesis and characterization of these stable, low-$\kappa_l$ compounds.


INTRODUCTION
An important focus in materials science research has been to discover novel materials with properties that might hold the keys to solving the most pressing problems in renewable energy, energy harvesting, or semiconductor power electronics. The augmentation of new materials discovery and the prediction of their properties have been accelerated by the advent of advanced computer algorithms coupled with high-throughput (HT) screening methods [1][2][3][4][5][6][7][8][9][10] using accurate quantum mechanical calculations based on density functional theory (DFT). In the recent past, several computational predictions have led to the successful synthesis of new solid-state compounds in a variety of chemistry and structure types in the family of half-Heuslers [9,10], double half-Heuslers [11], electrides [12], AB 2 X 4 based chalcogendes [13], and rocksalt based compounds [14].
Crystalline solids with extreme thermal transport properties are technologically important for the efficient management of thermal energy [15]. While materials with high lattice thermal conductivity (κ l ) are used in microelectronic devices for heat dissipation, materials with low κ l are used in thermal barrier coatings [16], thermal data-storage devices [17], and high-performance thermoelectrics (TEs) [18,19] which can convert heat into electrical energy. The conversion efficiency of the TEs is defined by the figure-of-merit, ZT=S 2 σ/(κ l + κ e ), where S, σ, κ e and κ l are the Seebeck coefficient, electrical conductivity, electronic thermal conductivity, and lattice thermal conductivity, respectively. Engineering the electronic band structures of crystalline compounds that already possess * koushik.pal.physics@gmail.com † c-wolverton@northwestern.edu low κ l has emerged to be a very popular strategy to increase the ZT. Therefore, crystalline semiconductors with intrinsically low-κ l are highly sought after in thermoelectrics and other thermal energy management devices.
In pursuit of finding new low-κ l materials, different classes of crystalline compounds like the perovskites [20], half-Heuslers [21], full-Heuslers [22], and double half-Heuslers [11], have been explored using HT computational methods, and subsequently, some of the predicted compounds were experimentally synthesized successfully [9][10][11]. Despite the existence of a large number of crystalline compounds in various materials databases e.g., the Inorganic Crystal Structure Database (ICSD) [23], Open Quantum Materials Database (OQMD) [3,4], Materials Project [5], and Aflowlib [24], it is quite important to look for novel stable and metastable compounds which might exhibit exciting physical and chemical properties. In this work, we present the computational discovery of a large number of stable (628) and low-energy metastable (852) quaternary chalcogenides AMM'Q 3 (A = alkali, alkaline earth, post-transition metals; M/M' = transition metals, lanthanides; Q = chalcogens) that span a huge chemical space across the periodic table. Our results are based on reliable, accurate and robust HT-DFT calculations where (a) we generated initial AMM'Q 3 compositions following the experimentally known AMM'Q 3 compounds formation criteria, (b) calculated their energetics in all known seven crystallographic prototypes that are found in this materials family, and (c) performed thermodynamic phase stability analysis of these compositions against all possible competing phases that are present in the OQMD.
About 192 quaternary chalcogenides (see Supplementary Information (SI) for a complete list) with the generic formula AMM'Q 3 have been synthesized experimentally [25][26][27][28][29][30][31][32][33] which reveal that these compounds possess rich chemistries and structure types like arXiv:2012.10778v1 [cond-mat.mtrl-sci] 19 Dec 2020 the perovskites and Heusler compounds. Koscielski et al. [25] noted that these known AMM'Q 3 compounds contain no Q-Q bonds and the elements (A, M, M', Q) balance the charge in their crystal structures with their expected formal oxidation states, making them charge balanced. The AMM'Q 3 compounds are further classified into three categories depending on the nominal oxidation states of the three cations A, M, and M', namely: where the oxidation states of the cations are indicated with the superscripts. In all cases, we assume that the chalcogen atoms Q (S, Se, and Te) have a nominal 2charge. Examining the experimentally known AMM'Q 3 compounds (see SI Figure S1), we observe the following: (a) The A site in these compounds is always occupied by alkali, alkaline earth, or post-transition metals with the only exception of Eu (in a 2+ charge state) which occupies the A-site of some of the Type-II compounds. (b) Whereas only transition metals occupy the M site, the M' site can be filled either by the transition metals, lanthanides, or actinides. (c) No observed AMM'Q 3 compound contains more than one alkali, alkaline earth, or post-transition metals. As detailed later, we will use these criteria in designing our HT workflow for generating the initial crystal structures through prototype decoration. Although the crystal chemistries of these compounds have been characterized somewhat in detail [25][26][27][28][29][30][31][32][33], their properties have remained largely unexplored. Recently, it was shown experimentally [34] and theoretically [34][35][36][37] that many known semiconducting compounds in this crystal family exhibit ultralow κ l . In addition, some compounds are shown to possess electronic bands favorable to support high thermoelectric performance [34,36,37].
From the distribution of the elements forming the known AMM'Q 3 compounds, we see that all the cations (A, M, M') are in their common oxidation states of 1+, 2+, 3+ or 4+ coming from the chemical groups (i.e., alkali, alkaline earth, transition, post-transition metals, lanthanides, actinides) that span a large part of the periodic table. Yet, although a large number (192) of these compounds have been reported experimentally, this number is small compared to the vast number of possible compounds that can be obtained based on chargebalanced combinatorial substitutions of the elements in a prototype crystal structure of AMM'Q 3 . Performing this combinatorial exercise experimentally would require a massive amount of resources and time to discover new AMM'Q 3 compounds. However, computational screening can be very helpful in narrowing down the search space of the target compounds that would have higher chances of synthesizability in the laboratory [1,[8][9][10][11][38][39][40][41][42][43]. Here, we have performed HT-DFT calculations followed by accurate ground-state phase stability analysis, and suggest (T = 0 K) thermodynamically stable and metastable AMM'Q 3 compounds for experimental synthesis and exploration of their properties. Our calculations of the thermal transport properties of some of the predicted stable compounds using the Peierls-Boltzmann transport equation (PBTE) show that these compounds exhibit innate low-κ l due to the presence of strong lattice anharmonicity or rattler cations.

Structural prototypes
The experimentally known AMM'Q 3 compounds crystallize in seven structure types [25][26][27][28][29][30][31][32][33]: KCuZrSe 3 (space group (SG): Cmcm, #63), Eu 2 CuS 3 (SG: Pnma, #62), BaCuLaS 3 (SG: Pnma, #62), Ba 2 MnS 3 (SG: Pnma, #62), NaCuTiS 3 (SG: Pnma, #62), BaAgErS 3 (SG: C2/m, #12) and TlCuTiTe 3 (SG: P2 1 /m, #11). All these structure types are visualized in their extended unit cells in Figure 1, where the conventional unit cell is outlined with black lines. Among these, five structure types (KCuZrSe 3 , Eu 2 CuS 3 , Ba 2 MnS 3 , NaCuTiS 3 , TlCuTiTe 3 ) are layered where the rows of A m+ cations stack alternatively with the [MM'Q 3 ] m− layers and interact through electrostatic interactions [35,37]. The strength of the interactions vary with the charges on the cations (i.e., m+) and induce modifications in the structures as well as in their properties [36]. When the interactions between the layers increase significantly, the atoms from the neighboring [MM'Q 3 ] m− layers interact, giving rise to the three-dimensional channel structures (BaCuLaS 3 and BaAgErS 3 ). In BaCuLaS 3 (Figure 1c) and BaAgErS 3 (Figure 1f), two rows of the A-site cations occupy the empty spaces inside the channels formed by the M, M', and Q atoms. Figure 1h shows that 71% of the known AMM'Q 3 compounds crystallize in the KCuZrSe 3 structure followed by 16% of the compounds crystallizing in the Eu 2 CuS 3 structure type. The rest of the known AMM'Q 3 compounds (13 %) crystallize in the other five structure types. It is worth noting that in Eu 2 CuS 3 , the Eu atoms have mixed oxidation states (Eu 2+ CuEu 3+ S 3 ) and sit in two different sites in the crystal structure. Similarly, the Ba atoms in Ba 2 MnS 3 also occupy two different sites. We have used all these structure types in our HT-DFT design and discovery of new AMM'Q 3 compounds. We note that the KCuZrSe3 and TlCuTiTe 3 structure types have 12 atoms in their primitive unit cells and the rest of the structure types have 24 atoms.

Materials design strategy
The discovery of new compounds through HT-DFT method often starts with the decoration of prototype crystal structures with chemically similar elements from the periodic table to generate their initial crystal structures. DFT calculations are then performed on these newly decorated compounds followed by rigorous thermodynamic phase stability analysis to screen for stable and metastable compounds. Rather than generating the input crystal structures in a brute-force manner by substituting every element of the periodic table at all atomic sites in the prototype structures, in this work we restricted ourselves by the following screening (c) We only consider compound compositions that are charge-balanced. (d) We exclude any radioactive elements during the prototype decorations although some of the known AMM'Q 3 compounds contain them. Adhering to these pre-conditions helps us narrow down our search space of compound exploration, reduces the computational cost, and most importantly it increases the success rates of stable compound prediction through HT-DFT calculations which will be evident later.
First, we generate the crystal structures of 4659 AMM'Q 3 compositions (see Table I Elemental distributions of the experimentally known and theoretically decorated AMMQ3 compounds. The AMM'Q3 compounds are categorized into three types depending on the nominal oxidation states of the cations: The color-coded periodic table indicates the group of elements (whose oxidation states are written in blue color) that are used for prototype decorations in our HT-DFT calculations. Elements that are denoted with the check marks constitute the experimentally known AMM'Q3 compounds. Salmon, cyan and yellow colors represent elements that occupy A, M/M' and Q sites, respectively. The grey color represents the elements that were not used for prototype decoration. Eu is the only element that occupies both A and M' site. All experimental and decorated compounds are charge-balanced.
this family and (ii) all experimentally known AMM'Q 3 compounds have low energies (within 50 meV/atom above the convex hull) in this structure type. After performing DFT calculations for all these compositions followed by T=0 K phase stability analysis, we kept only those compounds ( ∼ 1700) that have an energy within 50 meV/atom of the ground state convex hull and discarded the rest from our search space. The DFT relaxed structures of all these 1700 compounds retain the KCuZrSe 3 structure type. In the next step, we take these 1700 compositions and regenerate their crystal structures in each of six additional structure types to perform DFT calculations of 6 × 1,700 = 10,200 AMM'Q 3 compounds. Next, we perform T = 0K thermodynamic phase stability analysis of those 1700 compositions considering all seven structure types and their competing phases that are available in the OQMD. In the final step, we obtain 628 thermodynamically stable and 852 metastable hitherto unknown AMM'Q 3 compounds after performing a total number of 4,659 + 10,200 = 14,859 DFT calculations. The stable 628 compounds include 69 Type-I, 231 Type-II, and 328 Type-III compounds, and among them, a total number of 570 compounds possess finite band gaps. A schematic of the HT-DFT flowchart is shown in Figure 3a, and a summary is given in Table I.

Phase stability analysis
We now present a detailed analysis of the T= 0 K ground state phase stability of all (known and predicted) AMM'Q 3 compounds. We begin our assessment with the phase stability analysis of the experimentally known AMM'Q 3 . Out of 192 known compounds, we find only 119 compounds (Type-I: 39, Type-II: 30, and Type-III: 50) in OQMD before we performed any new calculations from this work. We designate these 119 compounds as Set-I. As the initial crystal structures of the experimentally known compounds in OQMD mostly come from the ICSD, the DFT calculations of Set-I compounds in OQMD were performed based on their experimental crystal structures taken from the ICSD. We designate the rest of the known 192-119 = 73 compounds (Type-I: 1, Type-II: 60, and Type-III: 12) as Set-II, which were not present in OQMD due to (a) the absence of their structures in the ICSD and (b) no previous HT-DFT calculations were performed based on the prototype decorations in this AMM'Q 3 family. However, our HT-DFT calculations of all the decorated AMM'Q 3 compositions include the experimentally known AMM'Q 3 compounds in Set-II. Hence, we will first analyze the phase stability of the Set-I compounds to see if our DFT calculations are able to correctly capture the energetics of the known compounds in Set-I and then we will utilize the phase stability data of Set-II compounds to validate the reliability of our approach for the discovery of new stable AMM'Q 3 Total no. of AMM'Q 3 compositions (Type-I + type-II + type-III) = 4659. All compounds are generated in the KCuZrSe 3 prototype structure 1700 compounds within 50 meV/atom above the convex hull  compounds based on prototype decoration through HT-DFT calculations.

Set-I compounds
As detailed in the Methods section and as well as in other references [1, 3-6, 24, 41, 44] the hull distance (hd) is a metric of the thermodynamic stability of a compound. If the formation energy of the compound breaks the convex hull, then it is considered to be thermodynamically stable with hd = 0, indicating the likelihood of its synthesizability. On the other hand, compounds with a small positive hull distance (typically within a few tens of meV/atom) are called metastable and may also be in some cases experimentally synthesized [39,40]. According to this criterion, all experimentally known AMM'Q 3 compounds should possess zero or small positive hull distances. Our analysis reveals that in Set-I, all (39) Type-I compounds and all but one (29) Type-II compounds have hd = 0, which is in line with our expectations. The one Type-II compound that has a small positive hull distance is Eu 2+ CuEu 3 +S 3 (hd = 37 meV/atom). Also, all but three (47) Type-III compounds of Set-I have hd = 0. These three Type-III compounds are CsCoYbS 3 (hd =192 meV/atom), CsCoYbSe3 (hd =151 meV/atom), CsZnYbSe 3 (hd=76 meV/atom). So, 115 of 119 experimentally synthesized compounds in Set-I are thermodynamically stable in the OQMD. Hence, stability is an excellent metric for the synthesizability of the predicted compounds. From this analysis, we also find that one of the Type-II compounds in Set-I, BaAgErS 3 , possesses a small positive hd of 12 meV/atom when its calculation is performed on the experimental crystal structure having C2/m SG (#12), indicating that it is metastable at 0 K in this structure. Interestingly, our HT-DFT calculations reveal that BaAgErS 3 is stable T = 0 K in the KCuZrSe 3 structure type (see SI Figure S2).

Set-II compounds
Since the 73 compounds in Set-II did not exist in OQMD before our HT-DFT calculations, we generate their crystal structures through prototype decorations as mentioned before to perform DFT calculations and T = 0 K phase stability analysis. As these compounds have already been synthesized experimentally, our DFT calculations and phase stability analysis provide a key test of our methodology. It also gives us an opportunity to examine how reliably our calculations can predict hitherto unknown stable AMM'Q 3 compounds. After performing the thermodynamic stability analysis, we found that among the 73 compounds in Set-II (type-I: 1, type-II: 60, type-III: 12), 64 compounds have hd = 0, and only 4 compounds (SrCuYbS3: hd = 84 meV/atom, BaCuYbTe 3 : hd=62 meV/atom, EuCuYbS 3 : hd=72 meV/atom, PbCuYbS 3 : hd =104 meV/atom) of Type-II and 5 compounds (CsMnYbSe 3 : hd=8 meV/atom, CsZnYbS 3 : hd=163 meV/atom, CsZnYbTe 3 : hd=71 meV/atom, RbZnYbSe 3 : hd=102 meV/atom, RbZnYbTe 3 : hd=82 meV/atom) of Type-III have positive hull distances.
Thus, our HT-DFT calculations based on the decorated structures successfully capture the stability of all but 13 (out of 192) of the experimentally known AMM'Q 3 compounds. The 13 compounds which are experimentally observed but with hd > 0 all contain Eu and Yb, and the calculated energetics of these compounds may originate from the choice of incorrect pseudopotentials as discussed later. These results and analysis give us strong confidence in designing and discovering new AMM'Q 3 compounds using the HT-DFT calculations and thermodynamic T=0 K stability as a metric for synthesizable compounds. We have provided the phase stability data of all 192 experimentally known AMM'Q 3 compounds in the SI.

Novel AMM'Q3 compounds
After performing the T= 0 K stability analysis of all newly decorated AMM'Q 3 compounds for which HT-DFT calculations are performed in all 7 structure types, we discover a large number of 628 (type-I: 69, type-II: 231, type-III:328) thermodynamically stable compounds that exclude the experimentally known 192 compounds. To put this number into perspective, the OQMD (containing more than 900,000 entries as of September 2020) has stable (hd = 0) 1161 full-Heuslers (SG: Fm3m, #225), 618 half-Heuslers (SG: F43m, #216), 353 cubic (SG: Pm3m, #221) perovskites, 242 orthorhombic (SG: Pnma, #62) perovskites. Similar to the experimentally known AMM'Q 3 compounds, the KCuZrSe 3 and Eu 2 CuS 3 structure types are the most common, which constitute 41% and 25% of the predicted stable compounds, respectively (Figure 3b). Among the other structure types, NaCuTiSe 3 and TlCuTiTe 3 are quite common, which constitute 19% and 9% of the predicted stable compounds. The rest 6% compounds have the BaCuLaS 3 and BaAgErS 3 structure types. However, we found no new stable compound in the Ba 2 MnS 3 structure type. Our analysis shows that 570 out of 628 compounds possess finite bands that range from 0.2 eV to 5.34 eV, among which most of the compounds have band gaps within 1.5 eV (Figure 3c). This is not surprising since all the decorated compositions are charge-balanced. In ad-dition, we found 852 potentially synthesizable metastable compounds (Type-I: 59, Type-II: 282, Type-III: 511) with small positive hull distances (i.e., 0 < hd ≤ 50 meV/atom). A summary of the HT-DFT calculations is shown in Table-I, and the lists of all predicted stable and metastable compounds are given in the SI.
We further examine the predicted stable compounds within each type in terms of their chemistries (sulfides, selenides, and tellurides) and the cations (A, M, M'). The results are displayed as bar charts in Figures 4, where the bar corresponding to an element represents the number of stable compounds that contain it. We see that: (a) the elements that form the stable compounds constitute almost the entire periodic table. (b) There are more Cu-containing compounds in Type-I and Type-II categories compared to Ag or Au. In Type-III, there are more Mn compounds than other M elements with a 2+ oxidation state. (c) The number of stable compounds increases from Type-I (69) to Type-II (231) to Type-III (328). This is not surprising given that the number of elements that can occupy the M and M' sites (satisfying the charge-neutrality criteria) increase from Type-I to Type-III compounds ( Table-I). A general trend that is noticeable across three types is that as we go from sulfides to selenides to tellurides, the number of compounds with smaller A cations decrease whereas the number increases with larger A cations. For example, there are two, two, and no Li-containing compounds in sulfides, selenides, and tellurides of Type-I compounds, respectively. Similarly, the number of compounds that have Sr increase from sulfides (7) to selenides (20) to tellurides (34) in Type-II.

Lattice thermal transport properties
We now focus on exploring the lattice thermal transport properties of the predicted stable AMM'Q 3 compounds. An accurate estimation of κ l of a compound within a first-principles DFT framework is computationally very expensive [45]. Hence, the calculations of κ l for all the predicted stable compounds would require a massive amount of computational resources. However, to demonstrate the thermal transport properties of our newly predicted stable compounds, we randomly select a handful of compounds with some criteria. We first screen for non-magnetic and semiconducting compounds, where the lattice contribution dominates the total thermal conductivity. Next, we search for those compounds which have the KCuZrSe 3 structure type as it possesses the highest crystal symmetry (SG: Cmcm, #63) and the smallest unit cell (12 atoms). Finally, we randomly select 10 compounds for κ l calculations. The selected compounds, which include sulfides, selenides and tellurides, are: CsCuZrS 3 , BaCuScSe 3 , BaCuScTe 3 , BaCuTbSe 3 , BaAgGdSe 3 , CsZnYS 3 , CsZnGdS 3 , CsZnScSe 3 , CsZnScTe 3 , CsCdYTe 3 . We note that this list also includes Type-I (the first in the list), Type-II (the next 4) as well as Type-III (the last 5) compounds. The electronic structures and phonon dispersions of these compounds are given in the SI.
We calculate the κ l of these 10 compounds using the PBTE (see Methods section) and present the results in Figure 5. We see that all these compounds exhibit very low κ l where the in-plane (κ ⊥ l ) and the cross-plane (κ l ) components are lower than 3 Wm −1 K −1 and 1.2 Wm −1 K −1 , respectively, for T ≥ 300 K. Here, κ ⊥ l is perpendicular to the stacking direction of the layers in the crystal structure of the AMM'Q 3 compounds and κ l is parallel to it. As a reference, we compare our calculated κ l with that of a prototypical thermoelectric material SnSe, which was experimentally shown to possesses low-κ l that leads to a very high thermoelectric figure-of-merit [46]. The measured [47] κ ⊥ l and κ l for single-crystalline stoichiometric samples of SnSe are 1.9 Wm −1 K −1 and 0.9 Wm −1 K −1 , respectively, at 300 K, which become much lower in the off-stoichiometric polycrystalline samples [46]. Further examination of the results in Figure 5 reveals that in terms of the anisotropy of the κ ⊥ l and κ l components, Type-I and Type-III compounds are quite similar, but Type-II compounds are different from the rest i.e., κ ⊥ l /κ l (Type-I/III) > κ ⊥ l /κ l (Type-II). The difference in the anisotropy of the properties arises from the fact that in Type-II compound the electrostatic attractions between A 2+ and [MM'Q 3 ] 2− layers are stronger than that of between A 1+ and [MM'Q 3 ] 1− layers in Type-I/III compounds. The stronger interlayer interactions give rise to a shorter interlayer distance in Type-II, which make κ ⊥ l and κ l less anisotropic.
To gain a deeper understanding, we examine the lattice dynamics and thermal transport properties of BaCuScTe 3 (Type-II) and CsCdYTe 3 (Type-III) in detail. Our analysis reveals that the underlying physical principles governing the low-κ l in Type-II compounds are different from the Type-I/III compounds. We start our analysis with the phonon dispersion of BaCuScTe 3 (Figure 6a) which shows the presence of very low-frequency acoustic (< 45 cm −1 ), and optical (∼ 15 cm −1 along X-S-R directions) phonon modes, which give rise to low phonon group velocities. The phonon dispersion and the phonon density of states (Figure 6b) of BaScCuTe 3 show that a strong hybridization exists between the phonon branches up to 100 cm −1 , where Ba and Te atoms have large contributions. Soft acoustic phonon branches give rise to very low sound velocities and a strong hybridization between the phonons at low energies enhances the phonon scattering phase space. Both these factors help suitably to give rise to a very low-κ l (κ ⊥ l = 1.7 Wm −1 K −1 and κ l = 0.76 Wm −1 K −1 at 300 K). On the other hand, the phonon dispersion of CsCdYTe 3 (Figure 6d) features nearly dispersion-less optical phonon branches along the X-S-R-Z directions in the Brillouin zone at low energies, which are the characteristics of rattler vibrations in the crystal structure. In addition, it also has soft acoustic phonon branches (< 35 cm −1 ). The calculated κ l of CsCdYTe 3 becomes ultralow with κ ⊥ l and κ l being 0.82 Wm −1 K −1 0.25 Wm −1 K −1 , respectively at 300 K.
Rattler phonon modes are highly localized, which strongly inhibit the transport of phonons, giving rise to ultralow-κ l in many crystalline solids such as TlInTe 2 [48], CsAg 5 Te 3 [49], etc. It was shown that the filler atoms in clathrates [50] and skutterudites [51] act as ideal rattlers, which give rise to dispersion-less phonon branches where the phonon frequencies remain highly localized having very small participation ratio (PR) values ∼ 0.2 (see Methods section). The PRs of the phonon modes of BaCuScTe 3 and CsCdYTe 3 are color-coded in Figures 6a and 6d, respectively. We see that while most of the low-energy phonon modes (< 100 cm −1 ) of BaCuScTe 3 have PR values close to 1, signifying the absence of phonon localization, the low-energy dispersionless phonon branches of CsCdYTe 3 have small PRs ( < 0.2), indicating their highly localized nature of the phonon modes. The phonon density of states (Figure 6e) also reveals that these localized phonons primarily arise from the Cs atoms (confined in 25-55 cm −1 ) that act as rattlers. Thus, our analysis shows that the ultralow-κ l in  CsCdYTe 3 is primarily caused by the localized vibrations of the rattling phonons.
We further investigate the origin of such poor thermal transport properties in BaCuScTe 3 in terms of a more fundamental material quantity, the lattice anharmonic-ity. Strong lattice anharmonicity is one of the important factors that induces very low-κ l in compounds like SnSe [52], NaSbTe 2 [53], etc. To estimate the strength of the intrinsic anharmonicity, we calculate the macroscopic average Gruneisen parameter, where γ i and C v,i are the Gruneisen parameter and specific heat capacity at constant volume for the i-th phonon mode. The calculated γ of BaCuScTe 3 (1.5) is larger than that of CsCdYTe 3 (1.2), signifying the presence of stronger anharmonicity in the former. These γ values are comparable to that of NaSbSe 2 (1.7), NaSbTe 2 (1.6), NaBiTe 2 (1.5) etc. compounds which are experimentally [53] shown to possess ultralow-κ l . Thus, we see that while the low-κ l in Type-II compounds is caused by the low-sound velocities and a stronger lattice anharmonicity, the presence of rattler cations in Type-III as well Type-I compounds (see SI Figure S4) are primarily responsible for inducing low-κ l in them.

Rattler Rattler
To examine which phonons are primarily responsible for the conduction of heat in BaCuScTe 3 and CsCdYTe 3 , we plot the cumulative-κ l and their first-order derivatives (κ l ) with respect to the frequency at T = 300 K. We see from Figure 6c that while the acoustic and lowenergy optical phonons up to 90 cm −1 mainly contribute to κ ⊥ l , κ l is primarily contributed by the phonons up to 45 cm −1 in BaCuScTe 3 . On the other hand, Figure  6f shows that while the acoustic as well as the optical phonons up to 100 cm −1 have large contributions towards κ ⊥ l , only the acoustic phonons (up to 35 cm −1 ) primarily carry heat for κ l in CsCdYTe 3 . We also notice that the anisotropy (κ ⊥ l /κ l = 3.3 at T = 300 K) in CsCdYTe 3 is much larger that of BaCuScTe 3 (κ ⊥ l /κ l = 2.2 at T = 300 K). The origin of this anisotropy can be attributed to the contrasting interlayer and intralayer interactions in BaCuScTe 3 and CsCdYTe 3 . For example, the analysis of the interatomic force-constants (IFCs) reveals that interlayer interactions in CsCdYTe 3 are much weaker (IFC (Cs−T e) = -0.333 eV/Å 2 ) compared to BaCuScTe 3 (IFC (Ba−T e) =-1.204 eV/Å 2 ), which makes the transport of optical phonons (above 25 cm −1 ) along the stacking direction (i.e., κ l ) of CsCdYTe 3 less effective. On the other hand, the intralayer interactions in CsCdYTe 3 are much stronger (IFC (Cd−T e) = -3.923 eV/Å 2 , IFC (Y −T e) : -2.216 eV/Å 2 ) than those of BaCuScTe 3 (IFC (Cu−T e) = -2.396 eV/Å 2 , IFC (Sc−T e) = -1.890 eV/Å 2 ). As a result, while phonons up to 90 cm −1 mainly carry the heat for κ ⊥ l in BaCuScTe 3 , in CsCdYTe 3 they are carried by phonons with frequencies up to 100 cm −1 very effectively.

DISCUSSION
We notice that in Set-I and Set-II, all the known compounds that have positive hull distances have either Yb 3+ or Eu 3+ cations in them. Given the fact OQMD used Yb 2 and Eu 2 PPs which are intended for compounds having Yb 2+ and Eu 2+ cations, those energetic results are somewhat suspect. However, we note that our HT-DFT calculations predict the stability of other rare-earth elements containing known compounds in this family correctly. Hence, the newly predicted compounds containing Eu 3+ and Yb 3+ should be taken with caution. We note that none of our predicted stable compounds contain these cations. However, 11 and 7 of the predicted metastable compounds contain Eu 3+ and Yb 3+ cations, respectively.
Concerning the experimental validation of our prediction, we note that the experimental synthesis of a large number of compounds is a daunting task. While automated synthesis of a large number of compounds is now possible through high-throughput experimental facilities [54], here, we suggest only four compounds (Table II) that can be immediately picked up for experimental verification of our prediction. These compounds contain no toxic element and their calculated lattice thermal conductivity is very low. In Table II, we provide their DFT calculated band gaps, γ, and average κ l calculated at 300 K. These calculated quantities along with κ l can be compared with the experimentally measured values of these materials. Finally, we note that κ l for each compound has been calculated using only three-phonon scattering processes. The inclusion of additional four-phonon scattering rates [55,56] and grain-boundary [57] limited phonon scattering mechanisms could further lower the calculated κ l in this family of compounds. Also, as the electronic structure of the compounds features nearly flat bands and multiple peaks near the valence/conduction band extrema ( Figure S3, SI), some of these compounds are expected to exhibit good potential for thermoelectric applications as well.
In summary, we use HT-DFT calculations to discover a large number of 628 thermodynamically stable quaternary chalcogenides (AMM'Q 3 ). As all compositions in this family are charge-balanced, our analysis shows that 570 of 628 compounds possess finite band gaps which vary between 0.2 and 5.34 eV. Our calculations of the thermal transport properties show that AMM'Q 3 compounds exhibit intrinsically very low κ l , and the anisotropy in κ l is much smaller in Type-II compounds compared to Type-I/III compounds. Our analysis further reveals that lowκ l in this family originates either due to the presence of rattling cations (in Type-I/III compounds) or stronger lattice anharmonicity (Type-II compound). While the rattler cations give rise to localized phonon modes that inhibit the propagation of phonons, a stronger lattice anharmonicity enhances the phonon scattering phase space, leading to a low-κ l . In addition, there exists a strong coupling between the acoustic and low-energy optical phonon modes in Type-II compounds which increases the phonon scattering rates of the heat-carrying phonons. Our work is thus interesting not just from the perspective of new materials discovery but also for finding the presence of low κ l in them, which hold promises for further research and possible applications in energy materials and related devices.

DFT calculations
We performed DFT calculations using the Vienna Ab-initio Simulation Package (VASP) [58] using the projector-augmented wave (PAW) [59] potentials with the Perdew-Burke-Ernzerhof (PBE) [60] generalized gradient approximation (GGA) to the exchange-correlation functional. The atomic positions and other cell degrees of freedom of the compounds were fully relaxed and spinpolarized calculations were performed for compounds that contain partially filled d or f-shells elements with a ferromagnetic arrangement of spins in accordance with the high-throughput framework as laid out in the qmpy suite of tools [3,4]. For more details on the calculation parameters, we refer to Ref. [3,4]. T=0 K phase stability analysis often serves as an excellent indicator for the possibility of synthesizability of a predicted compound in the laboratory [9][10][11]. To assess the T = 0 K thermodynamic stability of the compounds, we calculate their formation energies (∆H f ) utilizing the DFT total energy (ground state) of each compound using the formula, ∆H f (σ) = E(σ)i n i µ i , where E is the DFT total energy (at 0 K) of an AMM'Q 3 compound in a crystal structure denoted by σ, µ i is the chemical potential of element i with its fraction n i in that compound. For each composition, we used a number of prototype crystal structures, σ, based on structural prototypes of known AMM'Q 3 compounds. To determine the thermo-dynamic stability of a compound, we need to compare its formation energy against all its competing phases, not only against other compounds at the same composition. To this end, we generated the quaternary phase diagram (i.e., the T = 0 K convex hull of the A-M-M'-Q phase space) for every AMM'Q 3 compound considering all elemental, binary, ternary, and quaternary phases present in the OQMD, which (as of September 2020) corresponded to nearly 900,000 entries of DFT-calculated energies. The calculated convex hull distance (hd, explained in the next section) then serves as a metric to determine whether a compound is stable (i.e., hd = 0), metastable (small positive hd), or unstable (large positive hd).

Convex hull construction
To construct the convex hull of an AMM'Q 3 compound, it is necessary to identify the set of phases (elemental, binary, ternary as well as quaternary) in the fourdimensional composition space of A-M-M'-Q that have the lowest formation energies at their compositions. In OQMD, we present the convex hull of a quaternary compound through a four-dimensional phase diagram which is represented as an isometric shot of the Gibbs tetrahedron (Figure 7a). Each face of the tetrahedron represents a three-dimensional phase diagram of a ternary composition that is represented as the Gibbs triangle ( Figure  7b). The vertices in Figures 7a & 7b represent the elements constituting the quaternary and ternary phase space, respectively, and the edges connecting any two vertices represent the binary composition axis between those two elements. Any node within Figures 7a & 7b represents a stable compound, which is denoted by a cyan disk. The metastable and unstable compounds are not shown in these figures for clarity which fall off the nodes. The stability of an AMM'Q 3 compound is given by the difference (which is defined as the hull distance, hd = ∆H f -∆H e ) between the calculated formation energy (∆H f ) of an AMM'Q 3 compound under consideration and its hull energy (∆H e ) at that composition. The hull energy is defined as the energy at the convex hull at that AMM'Q 3 composition. By definition, the hd of a stable compound is zero, whereas for metastable and unstable compounds they are real positive numbers. In keeping with the heuristic conventions used in literature [10,61,62] we term those AMM'Q 3 compounds to be metastable whose hull distances lie within 50 meV/atom above the convex hull (i.e., 0 < hd ≤ 50 meV/atom). The metastable compounds are also potentially synthesizable in the laboratory [39,40].

Phonon participation ratio
Phonon dispersions have been calculated using 2 × 2 × 1 supercell of the primitive unit cell using Phonopy [63]. The high symmetry paths in the Brillouin zones were adopted following the conventions used by Setyawan et al. [64]. To examine the extent of localization of the phonon modes, we calculate their phonon participation ratio using the formula [50,65]: , where e i (ω q ) is the eigenvector of the phonon mode at q with frequency ω, M i is the mass of the i-th atom in the unit cell containing a total number of N atoms. The value of P (ω q ) ranges between 0 and 1. In an ordered crystal, when P (ω q ) becomes close to 1, it indicates that the phonon mode is propagative where all the atoms in the unit cell participate. On the other hand, very low values of PR (∼ 0.2) [65,66] indicate the strong localization of the phonon modes (e.g., rattling phonons) where only a few atoms in the unit cell participate in the vibrations. Examples of rattler atoms containing compounds include filled clathrates (Ba 8 Si 46 and Ba 8 Ga 16 Ge 30 ) [50,65], where the filler atoms act as ideal rattlers that induce ultralow-κ l in them.

Thermal conductivity calculations
We calculate the κ l utilizing the phonon lifetimes obtained from the third-order interatomic force constants (IFCs) [67][68][69], which was shown to reproduce κ l within 5 % of the experimentally measured κ l in this AMM'Q 3 family of compounds [34,35]. We constructed the thirdorder IFCs of each compound based on DFT calculations of displaced supercell configurations by limiting the cut-off distance (r c ) up to the third nearest neighbor. We used 2 × 2 × 1 supercell (containing 48 atoms) of the primitive unit cell (with 12 atoms) using thirdorder.py [45] utility for the calculation of IFCs. Using the second and third-order IFCs in the Sheng-BTE code [45], we calculate the temperature-dependent phonon scattering rates and κ l utilizing a full iterative solution to the PBTE for phonons using a 12 × 12 × 12 q-point mesh. The calculated κ l generally depends on the r c which accounts for the maximum range of interaction in the third-order IFCs [52]. It was shown that good convergence of κ l was obtained by limiting r c to the third nearest neighbor within the crystal structure in this family of compounds [35].

DATA AVAILABILITY
The data that supports the findings of the work are in the manuscript and Supplementary Information. The The four dimensional phase diagram (T= 0 K) of the Li-Cu-Zr-S quaternary system, which is the isometric shot of a Gibbs tetrahedron. (b) One of the faces of the tetrahedron represents the three-dimensional phase diagram of the Cu-Zr-S ternary system which is presented as the Gibbs triangle. Each cyan node in (a) and (b) represents a stable compound. For clarity, we do not show any metastable/unstable compounds and mark only few stable compounds in this figure. LiCuZrS3 is one of the predicted stable compounds in this work, which and its competing phases are denoted with dashed arrows in (a).
structures and energetics of the predicted compounds would be made available through the Open Quantum Materials database (OQMD) in a future release. Additional data will be available upon reasonable request.

CODE AVAILABILITY
Open-source codes are used throughout this work.

MATERIALS & CORRESPONDENCE
Request for materials and correspondence should be addressed to K.P. or C.W.