Structural Dynamics, Phonon Spectra and Thermal Transport in the Silicon Clathrates

Abstract

The potential of thermoelectric power to reduce energy waste and mitigate climate change has led to renewed interest in “phonon-glass electron-crystal” materials, of which the inorganic clathrates are an archetypal example. In this work we present a detailed first-principles modelling study of the structural dynamics and thermal transport in bulk diamond Si and five framework structures, including the reported Si Clathrate I and II structures and the recently-synthesised oC24 phase, with a view to understanding the relationship between the structure, lattice dynamics, energetic stability and thermal transport. We predict the IR and Raman spectra, including ab initio linewidths, and identify spectral signatures that could be used to confirm the presence of the different phases in material samples. Comparison of the energetics, including the contribution of the phonons to the finite-temperature Helmholtz free energy, shows that the framework structures are metastable, with the energy differences to bulk Si dominated by differences in the lattice energy. Thermal-conductivity calculations within the single-mode relaxation-time approximation show that the framework structures have significantly lower ${\kappa }_{\mathrm{latt}}$ than bulk Si, which we attribute quantitatively to differences in the phonon group velocities and lifetimes. The lifetimes vary considerably between systems, which can be largely accounted for by differences in the three-phonon interaction strengths. Notably, we predict a very low ${\kappa }_{\mathrm{latt}}$ for the Clathrate-II structure, in line with previous experiments but contrary to other recent modelling studies, which motivates further exploration of this system.

1. Introduction

Among the most important challenges to contemporary science is meeting the ever-increasing worldwide demand for energy whilst reducing anthropogenic greenhouse gas emissions. Achieving this goal will require a suite of technologies including new and improved sources of clean energy together with improvements to the efficiency of energy-intensive processes.

An estimated 60% of the energy used globally is currently wasted as heat, which is only estimated to improve to around 50% by 2030, and this has led to intensive research interest in thermoelectric power [1,2]. Thermoelectric generators (TEGs) harness the Seebeck effect in a thermoelectric material to directly convert waste heat into electrical power [3]. TEGs are solid-state devices with no moving parts, and as such have a broad range of potential applications from powering remote sensing devices, to extracting waste heat from combustion engine exhaust trains, to repurposing decommissioned oil rigs as geothermal power plants [3,4].

The performance of a thermoelectric material is determined by the dimensionless figure of merit, $ZT$: [1].

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

where S is the Seebeck coefficient, $\sigma$ is the electrical conductivity and $\kappa$ is the thermal conductivity, which is the sum of the electronic and lattice (phonon) components ${\kappa }_{\mathrm{elec}}$ and ${\kappa }_{\mathrm{latt}}$. S, $\sigma$ and ${\kappa }_{\mathrm{elec}}$ all depend on the electronic structure of the material and are linked through the carrier concentration, n, such that the best balance is typically found in heavily-doped semiconductors [1]. In particular, increasing n increases the conductivity but generally decreases the Seebeck coefficient and increases the electronic contribution to the thermal conductivity, which in heavily-doped semiconductors is often proportional to $\sigma$ through the Wiedemann–Franz law. The ${\kappa }_{\mathrm{latt}}$ depends on the structure and chemical bonding of the material, and low ${\kappa }_{\mathrm{latt}}$ is typically found in materials composed of heavy elements with weak chemical bonding [5] and/or strongly anharmonic lattice dynamics associated with structural features such as phase transitions [6,7,8] and active lone pairs [9].

The requirement to optimise both the electronic transport and the lattice thermal conductivity has historically led to a focus on “heavy” chalcogenide materials, which include the current industry standards for room- and high-temperature applications Bi₂Te₃ ($ZT$ ≃ 1 between 350–450 K) and PbTe ($ZT$ ≃ 2.2 at 915 K with endotaxial nanostructuring) [10]. However, an issue with these types of materials is the toxicity of heavy metals such as Pb and the rarity of Te, both of which mean these current “champion” materials are unsuitable for mass production and are therefore restricted to niche applications. This has led to wide-ranging research to identify more abundant and less toxic candidate thermoelectrics, including other chalcogenides, oxides and oxyclalcogenides alongside skutterudite and clathrate framework structures, intermetallics and half-Heuslers, among others [3,11,12,13].

Approaches to improving $ZT$ fall into two general classes, viz. optimising the power factor $S^2\sigma$ by tuning the electronic structure e.g., using band engineering [1,14,15], and reducing the lattice thermal conductivity e.g., through alloying, “discordant atom” doping, and hierarchical nanostructuring [10,16,17,18,19]. The interlinked nature of the transport properties makes optimisation challenging: for example, a high electrical conductivity requires a low carrier effective mass $m\ast$ (high carrier mobility $\mu$), whereas a large Seebeck coefficient is favoured by a large $m\ast$, and the conductivity and electronic thermal conductivity are typically proportional such that increasing $\sigma$ also increases ${\kappa }_{\mathrm{elec}}$ [1]. On the other hand, the lattice thermal conductivity is largely independent and can in principle be minimised without degrading the electrical transport [1]. However, the comparatively limited microscopic understanding of heat transport through phonons and how it can be controlled means that at present this is largely done through empirical strategies and trial and error.

The dual requirement for good electrical transport characteristic of crystalline semiconductors and poor thermal transport characteristic of amorphous materials (glasses) is encapsulated in the “phonon glass electron crystal” (PGEC) concept initially proposed by Slack [20]. A well-known subclass of PGEC compounds are materials based on a crystalline framework into which loosely-bound guest species can be incorporated. The “rattling” of the guest species in the cages couples to the framework modes and suppress the thermal transport through the lattice [21], while the electronic transport through the framework is unaffected, and the guest species can also be used to optimise the electrical properties e.g., by reducing or oxidising framework atoms to control the carrier concentrations. This idea has been used to good effect in CoSb₃, where the $ZT$ of 0.05 at 723 K shown by the pristine material [22] can be increased by filling with single ions (e.g., Na_0.48Co₄Sb₁₂ − $ZT$ = 1.32 at 800 K [23]) or multiple species (e.g., Ba_0.08La_0.05Yb_0.04Co₄Sb₁₂ − $ZT$ = 1.7 at 850 K [24]).

A second subclass of potential PGEC compounds are the inorganic clathrates typically formed by Group 14 elements. Clathrates are “open-structured” 3D frameworks composed of polyhedral cages, and the wide variety of possible structural motifs lead to a complex relationship between the cage structure, structural dynamics, and the thermal conductivity. Clathrates show potential for thermometric applications because of their transport properties, wide band gap [25,26,27] and reported low thermal conductivity [25,26,27]. Elemental Si adopts a bulk diamond structure, which can be easily doped to obtain a competitive power factor of around 50 $\mathsf{\mu }$W K${}^{-2}$ between ∼800–1200 K [28,29] but has a high ${\kappa }_{\mathrm{latt}}$ on the order of ∼25–45 W m${}^{-1}$ K${}^{-1}$ over the same temperature range [30,31]. However, Si can also form a number of clathrate phases [32,33,34] and other open-framework structures [35,36], at least one of which has been measuered to have a much lower thermal conductivity than diamond Si [37]. Ab inito calculations are playing an increasingly important role in building our microscopic understanding of lattice thermal conductivity, providing insight into both the origin of the intrinsically-low ${\kappa }_{\mathrm{latt}}$ of high-performance thermoelectrics such as SnSe [7,19,38] and also how material modifications such as alloying and filling can be used to optimise it [18,39,40]. Given the potential of clathrate structures in general, it is of considerable interest to use these techniques to unravel the link between the structure and thermal conductivity in the allotropes of Si.

In this work we utilise periodic density-functional theory (DFT) calculations to investigate the structural dynamics of bulk diamond Si and five framework structures drawn from experiments and previous theoretical studies [32,35,36]. We model and analyse the phonon spectra and compute reference infrared (IR) and Raman spectra for confirming the identity of the materials in experimental syntheses. We compare the energetic stabilities using both the athermal lattice energies and the constant-volume Helmholtz free energies including the phonon contributions to the internal energy and entropy. Finally, we calculate and analyse the lattice thermal conductivities of the six materials and elucidate the microscopic origin of the significant differences in the heat transport of the bulk and framework structures. From these calculations we obtain a detailed understanding of how moving from a dense bulk phase to open framework structures affects the structural dynamics and thermal transport, and we hope the improved fundamental understanding of the factors that determine the thermal transport in clathrates and related systems will inform the future design of high-performance thermoelectrics with ultra-low ${\kappa }_{\mathrm{latt}}$.

2. Computational Modelling

Periodic density-functional theory calculations were performed using the Vienna Ab initio Simulation Package (VASP) code [41].

As described in Section 3.1, initial structures of bulk diamond Si and five framework structures were obtained from the Materials Project database [42] and the Inorganic Crystal Structure Database (ICSD).

The PBEsol functional [43] was used to describe electron exchange and correlation, as we have previously found that this gives good results for lattice dynamics and thermal conductivity [44]. The Si atoms were described using projector augmented-wave (PAW) pseudopotentials [45,46] treating the Si 3s and 3p electrons as valence states. The Kohn-Sham wavefunctions of the valence electrons were represented in a plane-wave basis with a kinetic-energy cutoff of 500 eV. Brillouin-zone integrations were performed using the $\mathrm{\Gamma }$-centered Monkhorst-Pack $k$-point meshes [47] listed in Table 1 and a Gaussian smearing with a width $\sigma$ of 0.01 eV. The cutoff and $k$-point meshes were both chosen to converge the absolute total energy and cell pressure to <1 meV atom${}^{-1}$ and <1 kbar (0.1 GPa) respectively. The size of the charge-density grids was determined automatically to avoid aliasing errors (PREC = Accurate in VASP), and non-spherical contributions to the gradient corrections inside the PAW spheres were accounted for (LASPH = .TRUE.).

Each structure was symmetrised using spglib [48] and fully optimised to tolerances of 10${}^{-8}$ eV on the electronic total energy and 10${}^{-2}$ eV Å${}^{-1}$ on the forces.

Lattice-dynamics and thermal-conductivity calculations were set up and post-processed using the Phonopy and Phono3py software [49,50]. The second- and third-order force constants were determined using supercell finite-displacement calculations with the supercells shown in Table 1 and the default displacement step sizes of 10${}^{-2}$ and 3 × 10${}^{-2}$ Å respectively. Where required conversions between primitive and conventional cells were performed using spglib and/or by applying appropriate appropriate transformation matrices during post-processing. During the single-point force calculations additional support grids with 8× the number of points as the standard charge-density grids were employed to reduce numerical noise (ADDGRID = .TRUE.).

Phonon density of states (DoS) curves $g\left(\nu \right)$ were computed by interpolating the phonon frequencies onto the regular $\mathrm{\Gamma }$-centered $\mathit{q}$-point grids listed in Table 1. The DoS curves were broadened using a Gaussian smearing with a width $\sigma$ = 0.064 THz, corresponding to a full-width at half-maximim (FWHM) of 5 cm${}^{-1}$. These meshes were also used to compute the phonon contributions to the constant-volume (Helmholtz) free energy $A_{\mathrm{vib}}$. Phonon dispersion curves $\nu \left(\mathit{q}\right)$ were computed by interpolating frequencies along strings of $\mathit{q}$ passing between high-symmetry wavevectors in the respective Brillouin zones. The lattice thermal conductivities were computed using the single-mode relaxation-time approximation (SM-RTA) [50] from modal properties evaluated on the regular $\mathit{q}$-point grids listed in Table 1.

Finally, the infrared (IR) and Raman activities of the $\mathrm{\Gamma }$-point phonon modes were computed using the procedure outlined in Refs [51,52]. The Born effective-charge tensors ${\mathit{Z}}^{*}$ required for the IR activities were calculated using the density-functional perturbation theory (DFPT) routines in VASP [53]. The polarisability derivatives $\partial \mathit{\alpha }/\partial Q\equiv \partial {\epsilon }^{\infty }/\partial Q$ were computed using numerical differentiation with a step size $\Delta Q$ chosen such that the largest Cartesian displacement was 10${}^{-2}$ Å. The ${\epsilon }^{\infty }$ were computed using DFPT, with the exception of two modes of the C-II structure for which the DFPT calculation failed to converge and we therefore used instead the finite-field method outlined in Ref. [54] with a field strength of 10${}^{-3}$ eV Å${}^{-1}$.

3. Results and Discussion

3.1. Selection of Framework Structures

A search of the Materials Project (MP) database [42] yielded 39 entries for Si. 25 entries have Inorganic Crystal Structure Database (ICSD) codes and associated publications, of which 13 are labelled as being experimentally-determined and 12 as theoretical structures. The remaining 14 structures are user submissions with little other identifying information.

From the ICSD structures marked as experimental we selected the ambient diamond structure of bulk Si (mp-149), the Type-II clathrate structure published with Ref. [32] (mp-16220) and the “oC24” framework structure reported in Ref. [36] (mp-1095269). The oC24 structure is not formally designated as a clathrate but has a similar cage-like structure and is formed in a similar manner to the Type-I and Type-II clathrates [33,36,55]. The other structures include the hexagonal diamond structure that can be prepared in nanowire form [56,57] and a number of dense phases formed under pressure [58,59,60,61,62], neither of which are of interest to the present study.

The theoretical study in Ref. [35] proposes 11 so-called “Kelvin” framework structures. From these, we selected the three lowest-energy structures, viz. the ${\mathrm{Si}}_{\mathrm{K}}^{\mathrm{VII}}$ structure reported as the Type-V clathrate in Ref. [34] (mp-1203790), the ${\mathrm{Si}}_{\mathrm{K}}^{\mathrm{II}}$ structure identified as corresponding to the empty Type-I clathrate (mp-971662), and the ${\mathrm{Si}}_{\mathrm{K}}^{\mathrm{V}}$ structure similar to the Type-VI structure reported in Ref. [34]. The MP entry for the latter (mp-1196961) has a different spacegroup and cell volume to the linked ICSD structure, so we instead downloaded this structure from the ICSD (ICSD: 189396).

Table 2 lists the six structures selected for calculations together with their spacegroups, optimised lattice parameters, and the labels used to identify them. Corresponding images of the optimised structures are shown in Figure 1. The optimised lattice parmeters are generally in good agreement with experimental measurements and other theoretical literature. The room-temperature lattice constant of bulk Si is 5.431 Å [63], to which our calculated value of 5.436 Å is a good match. Our predicted lattice constants for the oC24 structure are a = 3.826, b = 10.70 and c = 12.66, which similarly compare well to the measured values of a = 3.82, b = 10.70 and c = 12.63 Å [36]. For the C-II structure we predict a lattice constant of 14.63 Å, which is very close to the measured 14.62 Å reported in the original structural characterisation in Ref. [64]. Our predicted volume of 23.02 Å${}^3$ atom${}^{-1}$ also lies between the values of 22.65 and 23.49 Å${}^3$ atom${}^{-1}$ obtained from the calculations in Ref. [32]. Likewise, our predicted volumes for the K-II/C-I, K-V/C-VI and K-VII/C-V structures are within 2% of the volumes calculated in Ref. [35].

3.2. Phonon Spectra

The calculated phonon dispersion $\nu \left(\mathit{q}\right)$ and density of states (DoS) curves $g\left(\nu \right)$ of the six structures are compared in Figure 2. The primitive cell of bulk diamond Si has a simple unit cell containing $n_{\mathrm{a}}$ = 2 atoms, resulting in 3$n_{\mathrm{a}}$ = 6 branches at each wavevector $\mathit{q}$ in the dispersion, and the high symmetry of the $Fd\overline{3}m$ spacegroup is such that two pairs of the branches are degenerate along the path shown in Figure 2a. In comparison, the more complex structures of the five framework materials result in larger primitive cells and thus more complex phonon dispersion curves, with 3$n_{\mathrm{a}}$ = 102 (C-II), 36 (oC24), 138 (K-II/C-I), 120 (K-V/C-VI) and 204 branches (K-VII/C-V).

The three acoustic modes in d-Si show a wide dispersion, reaching up to around 400 cm${}^{-1}$ at the zone-boundary wavevectors, while the optic modes span a relatively narrower range of ∼150 cm${}^{-1}$. In the five framework structures, the acoustic modes overlap with low-lying optic modes towards the zone boundaries, making it difficult to determine the frequency range precisely, but the acoustic modes do not appear to reach frequencies above ∼200 cm${}^{-1}$ and therefore have a considerably narrower dispersion than in d-Si. Each of the framework structures shows groups of optic modes both with relatively flat and with wider dispersions that are in some cases comparable to the optic modes in d-Si.

Despite the variations in complexity, the phonon frequencies in all six structures span a similar range of ∼500 cm${}^{-1}$. To assess the differences in the distribution of phonon frequencies over this range, we compared the cumulative integral of the phonon density of states as a function of frequency (Figure 3). Interestingly, we find that the phonon spectrum of the C-II structure shows a marked shift to lower frequencies compared to the other five systems, which otherwise have similar frequency distributions. This difference indicates that the C-II structure has significantly weaker chemical bonding compared to the other structures, a point to which we return in the following subsections.

While the dispersion and DoS can be measured using inelastic X-ray or neutron scattering (IXS/INS), these are relatively complex experiments that usually require access to synchrotron X-ray or neutron sources. Since the wavelengths of modes close to the Brillouin zone centre $\mathit{q}$ = $\mathrm{\Gamma }$ are comparable to those of visible and infrared (IR) light, these modes can also be probed using more routine IR and Raman spectroscopy. Simulated spectra can be generated from these calculations as a weighted sum of Lorentzian peak functions according to: [52]

(2)$I\left(\nu \right)=\sum _{\mathit{q}=\mathrm{\Gamma },j}\frac{I_j}{\pi }\frac{\frac{1}{2}{\mathrm{\Gamma }}_j}{{\left(\nu -{\nu }_j\right)}^2+{\left(\frac{1}{2}{\mathrm{\Gamma }}_j\right)}^2}$

where the sum runs over the 3$n_{\mathrm{a}}$ modes j at $\mathit{q}$ = $\mathrm{\Gamma }$ with ordinal frequencies ${\nu }_j$, spectroscopic activities $I_j$ and linewidths ${\mathrm{\Gamma }}_j$, and we have dropped the $\mathit{q}$ subscript from the terms in the summand for brevity.

The IR activity depends on the change in polarisation as the atoms are displaced along the mode. Given a $\mathrm{\Gamma }$-point mode eigenvector ${\mathit{W}}_j$ the Cartesian displacement ${\mathit{X}}_{j,k}$ of the kth atom is given by:

(3)${\mathit{X}}_{j,k}=\frac{1}{\sqrt{m_k}}{\mathit{W}}_{j,k}$

where $m_k$ is the atomic mass. The infrared activities $I^{\mathrm{IR}}$ are given by the change in polarisation along the mode eigenvector and are calculated from the Born effective-charge tensors of the atoms ${\mathit{Z}}_k^{*}$ using the expression: [51]

(4)$I_j^{\mathrm{IR}}=\sum _{\alpha }{\left|\sum _k\sum _{\beta }{Z}_{k,\alpha \beta }{X}_{j,k,\beta }\right|}^2$

where $\alpha$ and $\beta$ are the Cartesian directions x, y and z.

The Raman activity $I^{\mathrm{Raman}}$ depends on the change in polarisability $\mathit{\alpha }$, which can be computed from the derivative of the high-frequency dielectric constant ${\epsilon }^{\infty }$ with respect to the mode amplitude Q: [51]   

(5)$I_j^{\mathrm{Raman}}=\frac{V}{4\pi }\frac{\partial \mathit{\alpha }}{\partial Q_j}\equiv \frac{V}{4\pi }\frac{\partial {\epsilon }^{\infty }}{\partial Q_j}$

where V is the unit-cell volume. The Raman activities are tensor quantities and the $I_j$ therefore depend on the measurement geometry. To account for this, we perform an orientational average to obtain the scalar intensities ${\overline{I}}^{\mathrm{Raman}}$ that would be measured from a powder sample: [51]

(6)$\begin{array}{cc}\hfill {\overline{I}}^{\mathrm{Raman}}=& 45{\left[\frac{1}{3}\left(I_{11}+I_{22}+I_{33}\right)\right]}^2+\hfill \\ \hfill & \frac{7}{2}\left[{\left(I_{11}-I_{22}\right)}^2+{\left(I_{11}-I_{33}\right)}^2+{\left(I_{22}-I_{33}\right)}^2+6\left(I_{12}^2+I_{13}^2+I_{23}^2\right)\right]\hfill \end{array}$

where $I_{\alpha \beta }$ denote the components of the Raman tensor.

In this work we compute the linewidths ${\mathrm{\Gamma }}_j$ from the imaginary part of the phonon self-energy from many-body perturbation theory using the method documented in Ref. [50]. These are the intrinsic linewidths from energy- and momentum-conserving three-phonon scattering processes and do not include any instrument broadening.

The simulated IR and Raman spectra of the six materials are compared in Figure 4. The irreducible representations ${\mathrm{\Gamma }}_{3N}$ spanned by the $\mathrm{\Gamma }$-point modes in the six systems, based on which we identified the expected number of IR- and Raman-active modes, are summarised in Table 3.

The cubic $Fm\overline{3}m$ and $Pm\overline{3}m$ spacegroups adopted by d-Si, C-II and K-II/C-I both have the ${\mathrm{O}}_{\mathrm{h}}$ point group in which the T1u modes are IR active and the A1g, Eg and T2g modes are Raman active. In all three systems, the three acoustic modes form a T1u set. Although potentially IR active, since these modes correspond to rigid translations that cannot induce a change in polarisation they are constrained to have zero intensity. d-Si therefore shows no IR activity and a single peak in the Raman spectrum. We compute a central frequency of 508 cm${}^{-1}$ and a room-temperature (300 K) linewidth of 3.3 cm${}^{-1}$ for the Raman peak, which compare well to the experimental measurements of 520 cm${}^{-1}$ and 4.6 cm${}^{-1}$ respectively. [66] The much larger numbers of modes of the C-II and K-II/C-I structures result in more featured spectra with distinctive IR and Raman bands. The peak linewidths in the C-II spectrum are notably broader than in the other five spectra, which is indicative of stronger phonon anharmonicity and shorter lifetimes. We return to this point in Section 3.4.

The orthorhombic $Cmcm$ and $Cmmm$ spacegroups adopted by the oC24 and K-V/C-VI structures have the ${\mathrm{D}}_{2\mathrm{h}}$ point group where the B1u, B2u and B3u modes are IR active and the Ag, B1g, B2g and B3g modes are Raman active. The acoustic modes span a representation of B1u + B2u + B3u, so again a subset of the IR-active modes are restricted to have zero intensity. The oC24 structure has a considerably smaller unit cell than K-V/C-VI and as such has fewer modes, resulting in a notably less featured IR spectrum. On the other hand, although the K-V/C-VI structure potentially has ∼$3\times$ the number of Raman-active modes as the oC24 phase, differences in the relative intensities mean the spectrum of the latter shows more prominent features.

The $P{6}_3/mmm$ spacegroup of the K-VII/C-V structure has the ${\mathrm{D}}_{6\mathrm{h}}$ point group in which the A2u and E1u modes are IR active, the A1g, E1g and E2g modes are Raman active, and the acoustic modes span a representation of A2u + E1u. Of the six structures examined in this work K-III/C-V has the largest unit cell and the largest number of modes, although symmetry restrictions mean that there are fewer IR and Raman active modes than in the orthorhombic K-V/C-VI structure. The net result is that the spectra of these two structures are of comparable complexity. The Raman spectrum again has fewer prominent features than the oC24 structure due to differences in the relative intensities.

Overall, our simulations suggest that the differences in the six spectra would be sufficient to identify the different phases experimentally, especially using data from both IR and Raman spectroscopy. To aid with experimental characterisation and spectrum assignment, we provide a set of peak tables listing the central frequencies, irreducible representations, IR/Raman activities and linewidths calculated for the six structures as Supporting Information.

3.3. Thermodynamics

Within the harmonic approximation the temperature-dependent constant-volume (Helmholtz) free energy $A\left(T\right)$ is given by:

(7)$A\left(T\right)=U_{\mathrm{latt}}+A_{\mathrm{vib}}\left(T\right)=U_{\mathrm{latt}}+U_{\mathrm{vib}}\left(T\right)-T{S}_{\mathrm{vib}}\left(T\right)$

where $U_{\mathrm{latt}}$ is the lattice internal energy, equal to the DFT total energy obtained after structural optimisation, and $A_{\mathrm{vib}}=U_{\mathrm{vib}}-T{S}_{\mathrm{vib}}$ is the vibrational free energy from contributions of the phonons to the free energy and entropy.

$A_{\mathrm{vib}}$ can be computed directly from the vibrational partition function $Z_{\mathrm{vib}}$ using the bridge relation:

(8)$A_{\mathrm{vib}}\left(T\right)=-\frac{1}{N}{k}_{\mathrm{B}}T\ln \left[Z_{\mathrm{vib}}\left(T\right)\right]$

where N is the number of wavevectors $\mathit{q}$ included in the partition function product-equivalent to the number of unit cells in the crystal-so the normalisation yields the $A_{\mathrm{vib}}$ per unit cell. The phonon partition function $Z_{\mathrm{vib}}$ is computed from the (angular) phonon frequencies ${\omega }_{\mathit{q}j}$ according to:

(9)$Z_{\mathrm{vib}}\left(T\right)=\prod _{\mathit{q}j}\frac{\exp \left(-ħ{\omega }_{\mathit{q}j}/2k_{\mathrm{B}}T\right)}{1-\exp \left(-ħ{\omega }_{\mathit{q}j}/k_{\mathrm{B}}T\right)}$

where the product runs over phonon modes with wavevector $\mathit{q}$ and band index j. The vibrational internal energy $U_{\mathrm{vib}}$ can be computed according to:

(10)$U_{\mathrm{vib}}\left(T\right)=\sum _{\mathit{q}j}\left[n_{\mathit{q}j}\left(T\right)+\frac{1}{2}\right]ħ{\omega }_{\mathit{q}j}=\sum _{\mathit{q}j}\left[\frac{1}{\exp \left(ħ{\omega }_{\mathit{q}j}/k_{\mathrm{B}}T\right)-1}+\frac{1}{2}\right]ħ{\omega }_{\mathit{q}j}$

where $n_{\mathit{q}j}$ is the phonon occupation number given by the Bose-Einstein distribution. The phonon entropy $S_{\mathrm{vib}}$ can then be obtained from $A_{\mathrm{vib}}$ and $U_{\mathrm{vib}}$.

Table 4 compares the calculated lattice energies and the 300 K $A_{\mathrm{vib}}$ of the five framework structures to bulk diamond Si, while Figure 5 compares the free energy differences $\Delta A$ of the five structures to d-Si as a function of temperature from 0–1000 K. The calculated lattice energies of the five framework structures are all higher than that of d-Si, indicating them to be metastable. For four of the five framework structures, viz. oC24, K-II/C-I, K-V/C-VI and K-VII/C-V, the differences in lattice energy $\Delta U_{\mathrm{latt}}$ fall within a relatively narrow range of 6.82–9.14 kJ mol${}^{-1}$ atom${}^{-1}$, which suggests these phases are energetically similar. On the other hand, C-II is predicted to have a ∼4 × larger $U_{\mathrm{latt}}$ than the other framework structures, which can be taken as a reflection of the weaker chemical bonding evident in the phonon spectra. Comparing the 300 K free-energy differences $\Delta A$ gives similar results. The differences in internal energy $\Delta U$, including both lattice and phonon contributions, are comparable to the $\Delta U_{\mathrm{latt}}$, indicating that differences in vibrational internal energy are much smaller than differences in $U_{\mathrm{latt}}$. The differences in the vibrational entropy $-T\Delta S$ are generally negative, i.e., the framework structures are stabilised by a higher entropy, but the stabilisation of 0.12–0.92 kJ mol${}^{-1}$ atom${}^{-1}$ is again much smaller than the differences in the lattice energies.

To further investigate the impact of phonon contributions to the free energy on the relative stabilities of the framework structures, we compare the varition in $\Delta A$ with temperture to the differences in internal energy $\Delta U$ and entropy $-T\Delta S$ (Figure 6).

The $\Delta U$ all show slight increases with temperature, indicative of a larger $U_{\mathrm{vib}}$ that destabilises the framework materials relative to d-Si, but this effect is clearly secondary to the differences in the lattice internal energy. All five structures show a similar variation in internal energy with temperature, such that the $\Delta U$ between them remain roughly constant and the vibrational internal energy $U_{\mathrm{vib}}$ does not change the stability ordering predicted by the differences in the lattice energies. This is even the case for the orthorhombic oC24 and K-C/C-VI structures whose free energies are relatively close.

In contrast, differences in vibrational entropy favour the five framework structures relative to d-Si. For four of the five framework structures the stabilisation is <1 kJ mol${}^{-1}$ atom${}^{-1}$ at 1000 K, which is ∼5–10× smaller than the differences in the internal energy. For the C-II structure, on the other hand, the entropic stabilisation is much larger at around 3.5 kJ mol${}^{-1}$ atom${}^{-1}$. This can be attributed to the shift in the phonon spectrum to lower frequencies identified in Section 3.2 (c.f. Figure 3). However this does not compensate for the weaker bonding in this structure.

The study in Ref. [35] predicted the K-II/C-I, K-V/C-VI and K-VII/C-V structures to be 5.98, 7.62 and 5.31 kJ mol${}^{-1}$ atom${}^{-1}$ higher in energy than d-Si, respectively, which are comparable to our values of 7.47, 9.14 and 6.82 kJ mol${}^{-1}$ atom${}^{-1}$. On the other hand, Ref. [35] predicted the C-II structure to be 5.02 kJ mol${}^{-1}$ atom${}^{-1}$ higher in energy than bulk Si, which is a much smaller difference than predicted by our calculations and which would make this framework structure energetically competitive with the others. Similarly, our result is at odds with the computational study in Ref. [32], which predicted energy differences of 7.42 and 5.31 kJ mol${}^{-1}$ atom${}^{-1}$ with two different DFT functionals.

To investigate this further, we performed a series of additional calculations on d-Si and C-II to attempt to identify the source of the discrepancy (see Supplementary Materials). These calculations indicated that the higher energy of the C-II structure predicted by our calculations is not due to our choice of functional or technical parameters, and indeed we obtained very similar results to those in the Materials Project entry (mp-16220) using the PBE functional including a calculated $\Delta U_{\mathrm{latt}}$ of 32.8 kJ mol${}^{-1}$ atom${}^{-1}$. In our view, the most likely explanation is that the ICSD structure reported with Ref. [32] (ICSD: 56721), which is listed as the source of mp-16220, is different to the structure(s) used in previous modelling studies. However, this is the only C-II structure available in either of the ICSD and MP database.

Despite this potential issue, our phonon calculations confirm that the C-II structure we investigate here is dynamically stable by the absence of imaginary modes in the calculated phonon spectrum (c.f. Figure 2b).

The empty Type-II structure has reportedly been synthesised by heating NaSi to form Na₈Si₁₃₆, followed by repeated acid washing and heating under vacuum to remove the Na guest atoms [37,67]. Heating a mixture of elemental Na and Si under pressure is reported to yield the filled Type-I structure, Na₈Si₄₆, and NaSi₆, which can be converted to the empty oC24 by “thermal degassing” to remove the Na [36,55]. We would expect these preparation conditions potentially to yield metastable phases, and we therefore believe it feasible that the C-II structure examined here could in principle be formed. Indeed, as discussed in the following section, we predict a very low lattice thermal conductivity for this structure, which appears to be more in line with experimental measurements than predictions from previous computational studies [37,67,68,69], and which suggests this structure may be more representative of the phase obtained in experiments.

Given the K-V/C-VI and K-VII/C-V framework structures examined here are similar in energy to or lower in energy than the C-II, oC24 and K-II/C-I structures, which have (at least in principle) been prepared experimentally, this implies they should also be synthetically accessible given suitable synthesis conditions.

3.4. Thermal Transport

Using the single-mode relaxation-time approximation (SM-RTA) model, the macroscopic thermal conductivity ${\kappa }_{\mathrm{latt}}$ is computed as a sum of contributions from individual phonon modes as:

(11)${\kappa }_{\mathrm{latt}}=\frac{1}{N}\sum _{\mathit{q}j}{\kappa }_{\mathit{q}j}\left(T\right)=\frac{1}{NV}\sum _{\mathit{q}j}{C}_{\mathit{q}j}\left(T\right){\nu }_{\mathit{q}j}\otimes {\nu }_{\mathit{q}j}{\tau }_{\mathit{q}j}\left(T\right)$

Here N is the number of wavevectors in the summation and V is the unit-cell volume. The $C_{\mathit{q}j}$ are the modal heat capacities given by:

(12)$C_{\mathit{q}j}\left(T\right)=\sum _{\mathit{q}j}{\left[\frac{ħ{\omega }_{\mathit{q}j}}{k_{\mathrm{B}}T}\right]}^2\frac{\exp \left(ħ{\omega }_{\mathit{q}j}/k_{\mathrm{B}}T\right)}{{[\exp \left(ħ{\omega }_{\mathit{q}j}/k_{\mathrm{B}}T\right)-1]}^2}$

The ${\nu }_{\mathit{q}j}$ are the mode group velocities given by the derivative of the frequencies ${\omega }_{\mathit{q}j}$ with respect to the wavevector:

(13)${\nu }_{\mathit{q}j}=\frac{\partial {\omega }_{\mathit{q}j}}{\partial \mathit{q}}$

The ${\tau }_{\mathit{q}j}$ are the phonon lifetimes, which are related to the phonon linewidths ${\mathrm{\Gamma }}_{\mathit{q}j}$ by:

(14)${\tau }_{\mathit{q}j}\left(T\right)=\frac{1}{{\mathrm{\Gamma }}_{\mathit{q}j}\left(T\right)}$

As some of the clathrates have non-cubic spacegroups and therefore show anisotropic thermal transport, we also compute a scalar average ${\kappa }_{\mathrm{ave}}$ according to:

(15)${\kappa }_{\mathrm{ave}}=\frac{1}{3}\mathrm{Tr}\left({\kappa }_{\mathrm{latt}}\right)=\frac{1}{3}\left({\kappa }_{xx}+{\kappa }_{yy}+{\kappa }_{zz}\right)$

where ${\kappa }_{xx}$, ${\kappa }_{yy}$ and ${\kappa }_{zz}$ are the three diagonal elements of the ${\kappa }_{\mathrm{latt}}$ tensor, corresponding respectively to transport along the principal x, y and z directions, and we have omitted the explicit temperature dependence for brevity. This average is comparable to what would be measured in experiments on e.g., thin films or pressed pellets consisting of randomly-oriented crystal grains.

The calculated ${\kappa }_{\mathrm{ave}}$ of bulk diamond Si is shown as a function of temperature in Figure 7 and compared to experimental measurements from Refs [30,31]. The SM-RTA model gives a near-quantitative reproduction of the experiments over a range of ∼150–1000 K with a mean absolute relative error of 8.46%. Over this temperature range the largest difference occurs at 1000 K, where the calculated value of 35.6 W m${}^{-1}$ K${}^{-1}$ overestimates the experimental values of 29.8 and 31 W m${}^{-1}$ K${}^{-1}$ reported in Refs [30,31] by 19.6 and 14.9% respectively. Given the high temperature, this could be due to limitations of the model, including not accounting for changes in volume due to thermal expansion at finite temperature and/or the omission of higher- (e.g., fourth-) order phonon scattering processes when computing the linewidths. However, we also note that there are significant differences in the two sets of experimental measurements at some temperatures—e.g., at 300 K Refs [30,31] report values of 142.2 and 156 W m${}^{-1}$ K${}^{-1}$, which differ from our calculated value of 136.2 W m${}^{-1}$ K${}^{-1}$ by 4.2 and 12.7%, respectively. These discrepancies highlight the fact that measuring thermal conductivity experimentally can be somewhat challenging.

Comparison of the calculated ${\kappa }_{\mathrm{ave}}$ of the six systems (Figure 8) shows a clear grouping whereby d-Si >oC24 ≃ K-II/C-I ≃ K-V/C-VI ≃ K-VII/C-V > C-II. At room temperature, the calculated ${\kappa }_{\mathrm{ave}}$ of d-Si and C-II are 136.2 and 6.33 W m${}^{-1}$ K${}^{-1}$, respectively, while the conductivities of the other four framework structures range from 30.45 W m${}^{-1}$ K${}^{-1}$ (K-VII/C-V) to 57.86 W m${}^{-1}$ K${}^{-1}$ (oC24) (Table 5). With reference to the phonon spectra (Figure 2 and Figure 3), the high thermal conductivity of d-Si can be explained by the wide dispersion of the acoustic modes and consequent large group velocities ${\nu }_{\mathit{q}j}$, and the low thermal conductivity of the C-II structure could in principle be explained by its generally lower phonon frequencies, weaker chemical bonding, and resulting smaller group velocities.

Among the other four framework systems the ${\kappa }_{\mathrm{ave}}$ fall in the order of oC24 > K-II/C-I > K-V/C-VI > K-VII/C-V. This can be largely explained by the fact that high thermal conductivity is generally associated with less complex structures (i.e., with smaller primitive cells) and high symmetry: the oC24 and K-VII/C-V structures have the smallest and largest cells, respectively, while the K-II/C-I and K-V/C-VI structures have similar-sized primitive cells but K-II/C-I is cubic whereas K-V/C-VI is othorhombic (c.f. Table 2).

To better understand the differences in the thermal transport of the six systems we use the constant relaxation-time approximation (CRTA) model developed in Refs [18,39], whereby the ${\kappa }_{\mathrm{latt}}$ is written as the product of a harmonic term and a weighted-average lifetime ${\tau }^{\mathrm{CRTA}}$:

(16)${\kappa }_{\mathrm{latt}}\simeq {\tau }^{\mathrm{CRTA}}\times \frac{1}{N}\sum _{\mathit{q}j}\frac{{\kappa }_{\mathit{q}j}}{{\tau }_{\mathit{q}j}}={\tau }^{\mathrm{CRTA}}\times \frac{1}{NV}\sum _{\mathit{q}j}{C}_{\mathit{q}j}{\nu }_{\mathit{q}j}\otimes {\nu }_{\mathit{q}j}$

where we have again omitted the temperature dependence for brevity. Comparing the harmonic term in the summand and the ${\tau }^{\mathrm{CRTA}}$ allows the difference in ${\kappa }_{\mathrm{latt}}$ between systems to be attributed quantitatively to differences in the group velocities and lifetimes. We note that both the harmonic and lifetime terms are implicitly temperature dependent. While in this analysis we follow Ref. [39] and treat ${\tau }^{\mathrm{CRTA}}$ as a scalar quantity, the harmonic term is a tensor, and we therefore compute an equivalent scalar average to that in Equation (15) in order to analyse the ${\kappa }_{\mathrm{ave}}$.

The harmonoic and lifetime terms of the six systems as a function of temperature are compared in Figure 8b and Figure 8c respectively, and values at 300 K are included in Table 5. This analysis clearly shows that the differences in both the harmonic and lifetime terms are responsible for the variation in the ${\kappa }_{\mathrm{ave}}$. At 300 K the harmonic term for bulk Si is 5 W m${}^{-1}$ K${}^{-1}$ ps${}^{-1}$, which is the largest of the six compounds and can, as noted above, be explained by the wide dispersion of the acoustic modes shown in Figure 2. The oC24 structure has the second-highest $\kappa /{\tau }^{\mathrm{CRTA}}$ at 2.3 W m${}^{-1}$ K${}^{-1}$ ps${}^{-1}$ (54% smaller), and the remaining framework structures all have $\kappa /{\tau }^{\mathrm{CRTA}}$ < 1 W m${}^{-1}$ K${}^{-1}$ ps${}^{-1}$ (83–91% smaller). This ordering shows that the harmonic term is primarily governed by the size of the primitive cells such that larger cells result in lower $\kappa /{\tau }^{\mathrm{CRTA}}$.

The variation in the lifetime term is more intricate. At 300 K, d-Si has a ${\tau }^{\mathrm{CRTA}}$ of 27.2 ps. The oC24 framework has a 35% smaller averaged lifetime of 17.8 ps which, in addition to its smaller harmonic term, results in its 70% lower ${\kappa }_{\mathrm{ave}}$. On the other hand, the K-II/C-I, K-V/C-VI and K-VII/C-V structures have averaged lifetimes that are ∼1.5–2× longer than that of d-Si, which partially compensates for their smaller $\kappa /{\tau }^{\mathrm{CRTA}}$.

The low thermal conductivity of the C-II structure, on the other hand, appears to be characterised by both a low $\kappa /{\tau }^{\mathrm{CRTA}}$ and a short ${\tau }^{\mathrm{CRTA}}$. The former is the smallest among the six systems, at ∼10% of the value for d-Si, which can be explained by the impact of the weaker chemical bonding on the group velocities. The ${\tau }^{\mathrm{CRTA}}$ is also the shortest of the six compounds and is around 50% smaller than that of d-Si.

The phonon lifetimes are the inverse of the phonon linewidths ${\mathrm{\Gamma }}_{\mathit{q}j}$ (c.f. Equation (14)). Within the SM-RTA model the ${\mathrm{\Gamma }}_{\mathit{q}j}$ are calculated as a sum of contributions from three-phonon scattering processes between triplets of modes $\mathit{q}j$, ${\mathit{q}}^{\prime }{j}^{\prime }$ and ${\mathit{q}}^{″}{j}^{″}$ according to: [50]

(17)$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_{\mathit{q}j}\left(T\right)=& \frac{18\pi }{{ħ}^2}\sum _{{\mathit{q}}^{\prime }{j}^{\prime },{\mathit{q}}^{″}{j}^{″}}{\left|{\mathrm{\Phi }}_{-\mathit{q}j,{\mathit{q}}^{\prime }{j}^{\prime },{\mathit{q}}^{″}{j}^{″}}\right|}^2\times \hfill \\ \hfill & \{\left[n_{{\mathit{q}}^{\prime }{j}^{\prime }}\left(T\right)-n_{{\mathit{q}}^{″}{j}^{″}}\left(T\right)\right]\left[\delta \left({\omega }_{\mathit{q}j}+{\omega }_{{\mathit{q}}^{\prime }{j}^{\prime }}-{\omega }_{{\mathit{q}}^{″}{j}^{″}}\right)-\delta \left({\omega }_{\mathit{q}j}-{\omega }_{{\mathit{q}}^{\prime }{j}^{\prime }}+{\omega }_{{\mathit{q}}^{″}{j}^{″}}\right)\right]+\hfill \\ \hfill & \left[n_{{\mathit{q}}^{\prime }{j}^{\prime }}\left(T\right)+n_{{\mathit{q}}^{″}{j}^{″}}\left(T\right)+1\right]\delta \left({\omega }_{\mathit{q}j}-{\omega }_{{\mathit{q}}^{\prime }{j}^{\prime }}-{\omega }_{{\mathit{q}}^{″}{j}^{″}}\right)\}\hfill \end{array}$

The contribution from each triplet is determined by the three-phonon interaction strengths ${\mathrm{\Phi }}_{-\mathit{q}j,{\mathit{q}}^{\prime }{j}^{\prime },{\mathit{q}}^{″}{j}^{″}}$, the occupation numbers $n_{\mathit{q}j}$, and selection rules enforcing the conservation of energy. The two terms in the brace in the summation describe two types of scattering process, viz. collisions (Class 1), where two phonons coalesce to form a third, and decay (Class 2), where a phonon decays and emits two others. (There is also a selection rule for the conservation of (crystal) momentum, which is incorpared into the definition of the ${\mathrm{\Phi }}_{-\mathit{q}j,{\mathit{q}}^{\prime }{j}^{\prime },{\mathit{q}}^{″}{j}^{″}}$ [50]).

The three-phonon interaction strengths are temperature independent and calculated from the harmonic frequencies and eigenvectors ${\omega }_{\mathit{q}j}$/${\mathit{W}}_{\mathit{q}j}$ and the third-order (anharmonic) force constants ${\mathbf{\Phi }}^{\left(3\right)}$, while the conservation of energy expression is calculated solely from the ${\omega }_{\mathit{q}j}$ and is temperature dependant through the $n_{\mathit{q}j}$. Based on the analysis in Ref. [50], we approximate the linewidths as:

(18)${\mathrm{\Gamma }}_{\mathit{q}j}\left(T\right)\simeq \frac{18\pi }{{ħ}^2}\tilde{P}{N}_2({\mathit{q}}_j,{\omega }_{\mathit{q}j},T)$

where $\tilde{P}$ is a weighted squared average three-phonon interaction strength and $N_2(\mathit{q},\omega ,T)$ is a weighted two-phonon joint density of states (w-JDoS) defined as:

(19)$N_2(\mathit{q},\omega ,T)=N_2^{\left(1\right)}(\mathit{q},\omega ,T)+N_2^{\left(2\right)}(\mathit{q},\omega ,T)$

with:

(20)$\begin{array}{cc}\hfill N_2^{\left(1\right)}(\mathit{q},\omega ,T)=& \frac{1}{N}\sum _{{\mathit{q}}^{\prime }{j}^{\prime },{\mathit{q}}^{″}{j}^{″}}\Delta \left(-\mathit{q}+{\mathit{q}}^{\prime }+{\mathit{q}}^{″}\right)\left[n_{{\mathit{q}}^{\prime }{j}^{\prime }}\left(T\right)-n_{{\mathit{q}}^{″}{j}^{″}}\left(T\right)\right]\times \hfill \\ \hfill & \left[\delta \left(\omega +{\omega }_{{\mathit{q}}^{\prime }{j}^{\prime }}-{\omega }_{{\mathit{q}}^{″}{j}^{″}}\right)-\delta \left(\omega -{\omega }_{{\mathit{q}}^{\prime }{j}^{\prime }}+{\omega }_{{\mathit{q}}^{″}{j}^{″}}\right)\right]\hfill \end{array}$

(21)$\begin{array}{cc}\hfill N_2^{\left(2\right)}(\mathit{q},\omega ,T)=& \frac{1}{N}\sum _{{\mathit{q}}^{\prime }{j}^{\prime },{\mathit{q}}^{″}{j}^{″}}\Delta \left(-\mathit{q}+{\mathit{q}}^{\prime }+{\mathit{q}}^{″}\right)\times \hfill \\ \hfill & \left[n_{{\mathit{q}}^{\prime }{j}^{\prime }}\left(T\right)+n_{{\mathit{q}}^{″}{j}^{″}}\left(T\right)+1\right]\delta \left(\omega -{\omega }_{{\mathit{q}}^{\prime }{j}^{\prime }}-{\omega }_{{\mathit{q}}^{″}{j}^{″}}\right)\hfill \end{array}$

For a given mode with a wavevector $\mathit{q}$ and frequency $\omega$ these functions count respectively the number of possible energy- and momentum-conserving Class 1 and Class 2 events.

As per Equation (11) the ${\kappa }_{\mathrm{latt}}$ is proportional to the ${\tau }_{\mathit{q}j}$ and therefore inversely proportional to the ${\mathrm{\Gamma }}_{\mathit{q}j}$, so that when the squared three-phonon interaction strengths are replaced by a constant $\tilde{P}$ the ${\kappa }_{\mathrm{latt}}$ is proportional to $1/\tilde{P}$. A suitable $\tilde{P}$ to reproduce the scalar average ${\kappa }_{\mathrm{latt}}$ can therefore be obtained from a linear fit of the ${\kappa }_{\mathrm{latt}}$ as a function of $1/\tilde{P}$. As noted in Ref. [39], although the ${\mathrm{\Phi }}_{-\mathit{q}j,{\mathit{q}}^{\prime }{j}^{\prime },{\mathit{q}}^{″}{j}^{″}}$ in Equation (17) are temperature independent, replacing them with a single value cannot in general capture the temperature dependence of the ${\kappa }_{\mathrm{latt}}$, so the $\tilde{P}$ are weakly temperature dependent and the fit must be performed at a chosen temperature T.

Both the $\tilde{P}$ and $N_2(\mathit{q},\omega ,T)$ can be used to compare between materials and attribute differences in the phonon linewidths/lifetimes to differences in the phonon interaction strengths, which are anharmonic by definition, and to differences in the numbers of allowed scattering pathways dictated by differences in the harmonic phonon spectra. To compare the $N_2(\mathit{q},\omega ,T)$ it is prudent to average over $\mathit{q}$ to obtain a function of frequency only, i.e.,:

(22)$\begin{array}{cc}\hfill {\overline{N}}_2(\omega ,T)& ={\overline{N}}_2^{\left(1\right)}(\omega ,T)+{\overline{N}}_2^{\left(2\right)}(\omega ,T)\hfill \\ \hfill & =\frac{1}{N}\sum _{\mathit{q}}{N}_2^{\left(1\right)}(\mathit{q},\omega ,T)+\frac{1}{N}\sum _{\mathit{q}}{N}_2^{\left(2\right)}(\mathit{q},\omega ,T)\hfill \end{array}$

We note that both the $N_2(\mathit{q},\omega ,T)$/${\overline{N}}_2(\omega ,T)$ and $\tilde{P}$ are “extensive”, i.e., they scale with the number of atoms in the primitive cell, and that to compare between materials the w-JDoS must be scaled down by ${\left(3n_{\mathrm{a}}\right)}^2$ and the $\tilde{P}$ scaled up by ${\left(3n_{\mathrm{a}}\right)}^2$.

The scaled ${\overline{N}}_2\left(\omega \right)$ of the six structures at T = 300 K are compared in Figure 9a. The w-JDoS of all of the structures apart from C-II have a very similar shape, particular at lower frequencies, suggesting there are similar numbers of allowed scattering pathways. On the other hand, the ${\overline{N}}_2\left(\omega \right)$ of C-II indicates that there are a considerably larger number of scattering pathways available across most of its frequency spectrum. The shape of the w-JDoS depends on the shape of the harmonic phonon spectrum, and this therefore reflects the different distribution of frequencies in the C-II spectrum compared to the other structures shown in Figure 3.

A value of $\tilde{P}$ was fitted for each of the six structures at T = 300 K (see Supplementary Materials), and the scaled values are listed alongside the ${\kappa }_{\mathrm{latt}}$ and CRTA components in Table 5. Excluding C-II, a plot of $1/\tilde{P}$ against the averaged lifetime ${\tau }^{\mathrm{CRTA}}$ (Figure 9b) shows a strong linear correlation across all the framework structures ($R^2$ = 0.997), providing clear evidence that differences in the ${\tau }^{\mathrm{CRTA}}$ can predominantly be attributed to differences in the phonon-phonon interaction strengths. The $\tilde{P}$ fitted for the C-II structure lie slightly below the trendline, i.e., its lifetimes are lower than would be expected considering only the strength of the phonon-phonon interactions, and this can be accounted for by its generally larger ${\overline{N}}_2\left(\omega \right)$ (c.f. Figure 9a).

We therefore attribute the variation in the phonon lifetimes in the framework structures relative to bulk diamond Si predominantly to differences in the phonon-phonon interaction strengths, and the short ${\tau }_{\mathit{q}j}$ of the C-II structure to a combination of strong phonon interactions plus a larger number of allowed scattering pathways. This conclusion is in line with the modelling study in Ref. [69], which also highlighted the importance of phonon-phonon interaction strengths in determining the thermal conductivities of d-Si and the Si clatrates.

Our calculated room-temperature ${\kappa }_{\mathrm{ave}}$ of 6.33 W m${}^{-1}$ K${}^{-1}$ is of a similar magnitude to the value of 2.5 W m${}^{-1}$ K${}^{-1}$ reported for the empty Type-II clathrate, Ref. [37] and to the values of ∼2.5 and 1 W m${}^{-1}$ K${}^{-1}$ reported for Na₁Si₁₃₆ and Na₈Si₁₃₆ respectively [67]. However, our value is considerably smaller than the values of 52 and 57.5 W m${}^{-1}$ K${}^{-1}$ obtained from calculations similar to those in the present work. [68,69] On the other hand, the study in Ref. [69] obtained values of 151 and 53 W m${}^{-1}$ K${}^{-1}$ for d-Si and the Type-I clathrate structure, which are much more comparable to, albeit slightly higher than, our predicted values. Similarly, the study in Ref. [68] predicted a ${\kappa }_{\mathrm{latt}}$ on the order of 100 W m${}^{-1}$ K${}^{-1}$ for bulk Si, which is again comparable to the results obtained in the present work.

To try to elucidate the origin of the discrepancy with previous modelling studies we performed some additional calculations of the phonon dispersion and ${\kappa }_{\mathrm{latt}}$ (see Supplementary Materials), based on which we again concluded that the most likely explanation is the C-II structure reported with Ref. [32] being different to the structure(s) used in the modelling studies in Refs [68,69]. On the face of it, our calculated ${\kappa }_{\mathrm{latt}}$ is much more in keeping with the low thermal conductivity measured experimentally [37,67]. However, the low experimental values may be due in part to high concentrations of structural defects, which are difficult to avoid given the harsh processing conditions required to remove the Na guest atoms, and indeed the low density of the Na-filled samples in the study in Ref. [67] may be indicative of similar issues in these measured thermal conductivities. Although these issues were raised in Ref. [69] as a possible reason for the large difference between the calculations and the measurements in Ref. [37], it is debatable whether they could account for the ∼20× difference.

Finally, we also consider the anisotropy in the thermal conductivities of the three framework structures with non-cubic spacegroups. Figure 10 shows the temperature dependence of the ${\kappa }_{xx}$, ${\kappa }_{yy}$ and ${\kappa }_{zz}$ components of the ${\kappa }_{\mathrm{latt}}$ tensors in the oC24, K-V/C-VI and K-VII/C-V frameworks. For the orthorhombic oC24 and K-V/C-VI structures the three ${\kappa }_{\mathrm{latt}}$ components correspond to transport along the crystallographic a, b and c axes. The oC24 structure shows relatively large anisotropy. The largest ${\kappa }_{xx}$ = 54.86 W m${}^{-1}$ K${}^{-1}$ is along the shortest crystallographic a axis, which makes sense as the short lattice constant of 3.826 Å would enforce a high degree of regularity in the bonding along this direction (c.f. Table 2). On the other hand, the b and c axes have similar lengths of b = 10.697 and 12.660 Å, respectively, but the corresponding ${\kappa }_{\mathrm{latt}}$ differ by a factor of two. The anisotropy in the K-V/C-VI structure is much smaller, with the ${\kappa }_{\mathrm{latt}}$ components ranging from ${\kappa }_{xx}$ = 38.24 W m${}^{-1}$ K${}^{-1}$ to ${\kappa }_{zz}$ = 32.96 W m${}^{-1}$ K${}^{-1}$. Although the largest transport is again associated with the shortest a axis, the transport along the longest b axis is similar and the smallest ${\kappa }_{\mathrm{latt}}$ occurs along the c axis, which has a similar length to the a axis. Finally, the K-VII/C-V system has a hexagonal spacegroup in which the ${\kappa }_{xx}$ and ${\kappa }_{yy}$ components are equal and smaller than the ${\kappa }_{zz}$ component, corresponding to transport along the longer c axis, by around 2 W m${}^{-1}$ K${}^{-1}$.

4. Discussion

In this work we have performed a detailed study of the structural dynamics, energetics and thermal transport of bulk diamond Si and five Si framework structures including four clathrates.

From the calculated phonon spectra, we find that four of the five framework structures show a similar distribution of phonon frequencies to bulk Si, despite their more complex dispersions, whereas the spectrum of the Clathrate II structure is notably skewed towards lower frequencies. We have provided reference infrared and Raman spectra for each of the six systems, comparison of which indicates that it should be possible to use a combination of both techniques to establish or confirm the identity of phases prepared in experiments.

The energy differences predicted by the calculations indicate that all five framework structures are metastable with respect to bulk Si. For the majority of the frameworks the energy differences are on the order of 5–10 kJ mol${}^{-1}$ atom${}^{-1}$, but the C-II phase is predicted to be much less stable at >30 kJ mol${}^{-1}$ atom${}^{-1}$ above d-Si. While the framework structures, and in particular C-II, are stabilised at elevated temperature due to their higher vibrational entropy, this is insufficient to compensate for differences in the lattice energies and therefore does not affect the energetic ordering. Although the predicted high energy of the C-II phase is at odds with other computational studies, the calculated phonon dispersion indicates it to be dynamically stable, and the experimental synthesis technique employed to form it would in principle be capable of yielding such high-energy metastable phases.

All of the framework structures are predicted to have substantially smaller lattice thermal conductivity than d-Si, which can be attributed quantitatively to a reduction in the phonon group velocities and an intricate variation in the phonon lifetimes. The differences in group velocity correlate to the size and, to a lesser extent, symmetry of the primitive cells. The variation in the phonon lifetimes is predominantly due to differences in the strength of the interactions between phonon modes, indicting that the different structures show different levels of phonon anharmonicity.

Notably, we predict a very low room-temperature ${\kappa }_{\mathrm{latt}}$ of 6.33 W m${}^{-1}$ K${}^{-1}$ for the C-II structure, which is an order of magnitude smaller than predicted by two previous computational studies but appears to be a good match to experimental measurements. The limited heat transport in this framework structure can be attributed to low group velocities, reflecting weak bonding evident in its phonon spectrum, and to short lifetimes arising from a combination of strong phonon interactions and an increased number of allowed three-phonon scattering processes compared to the other systems.

Overall, this study highlights the utility of first-principles calculations for studying the structural dynamics and heat transport in complex framework structures, and in particular for establishing the origin of differences in lattice thermal conductivity. The latter is of particular interest in the field of thermoelectric materials, as it may provide valuable guidance on appropriate materials engineering strategies for optimising the ${\kappa }_{\mathrm{latt}}$. In this direction, the low ${\kappa }_{\mathrm{latt}}$ of the C-II structure predicted in this study and obtained in experiments suggests that further studies on this phase are warranted and may prove fruitful.

Author Contributions

Conceptualization: J.M.S.; methodology: J.M.F. and J.M.S.; formal analysis: B.W., J.M.F. and J.M.S.; investigation: B.W., J.M.F. and J.M.S.; resources: J.M.F. and J.M.S.; data curation: B.W. and J.M.S.; writing: J.M.F. and J.M.S.; supervision: J.M.F. and J.M.S.; funding acquisition: J.M.S. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Data Availability Statement

Raw data from this study, including optimised structures, input and output files for the Phonopy and Phono3py codes, and sample input files for the Vienna Ab initio Simulation Package code, will be made available to download free of charge from an online repository at https://doi.org/10.17632/9hywbzt8zd, accessed on 20 September 2022.

Acknowledgments

A subset of the calculations in this work were performed on the University of Manchester Computational Shared Facility (CSF) HPC system, which is maintained by UoM Research IT.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Images of the six structures examined in this work: (a) diamond Si (d-Si); (b) the Type-II clathrate (C-II); (c) the “oC24” structure reported in Ref. [36]; (d) the Type-I clathrate reported as [Forumla omitted. See PDF.] in Ref. [35] (K-II/C-I); (e) the [Forumla omitted. See PDF.] structure reported in Ref. [35], which is similar to the Type-VI clathrate [34] (K-V/C-VI); and (f) the Type-V clathrate identified as [Forumla omitted. See PDF.] in Ref. [34] (K-VII/C-V). These images were prepared using the VESTA software [65].

Figure 2. Calculated phonon dispersion [Forumla omitted. See PDF.] and density of states [Forumla omitted. See PDF.] curves of the six structures examined in this work: (a) d-Si, (b) C-II, (c) oC24, (d) K-II/C-I, (e) K-V/C-VI, and (f) K-VII/C-V.

Figure 3. Cumulative integral of the phonon density of states curves [Forumla omitted. See PDF.] in Figure 2 as a function of frequency showing the distributions of modes over frequencies.

Figure 4. Simulated infrared (IR; red) and Raman (blue) spectra of the six structures examined in this work obtained using the method in Ref. [52]: (a) d-Si, (b) C-II, (c) oC24, (d) K-II/C-I, (e) K-V/C-VI, and (f) K-VII/C-V. The Raman activities used to generate these spectra are the scalar averages [Forumla omitted. See PDF.] defined in Equation (6).

Figure 5. Calculated Helmholtz free energy differences [Forumla omitted. See PDF.] as a function of temperature for the five Si framework structures relative to bulk diamond Si. Note the break in the y axis.

Figure 6. Breakdown of the Helmholtz free energy differences [Forumla omitted. See PDF.] (a), shown in Figure 5, into differences in the internal energy [Forumla omitted. See PDF.] (b) and entropy [Forumla omitted. See PDF.] (c).

Figure 7. Calculated averaged lattice thermal conductivity of bulk diamond Si as a function of temperature (black line) compared to experimental data from Refs [30,31] (blue triangles/red circles).

Figure 8. Calculated averaged lattice thermal conductivity of the six structures examined in this work: d-Si (black), C-II (blue), oC24 (red), K-II/C-I (green), K-V/C-VI (purple) and K-VII/C-V (orange). Subplot (a) shows the [Forumla omitted. See PDF.] as a function of temperature, and subplots (b,c) show its decomposition into a harmonic term [Forumla omitted. See PDF.] and averaged lifetime [Forumla omitted. See PDF.] using the constant-relaxation time approximation (CRTA) model defined in Equation (16) in the text.

Figure 9. Analysis of the phonon lifetimes of the six systems examined in this work using the model defined in Equation (18) in the text. Plot (a) compares the scaled weighted two-phonon joint density of states (w-JDoS) functions [Forumla omitted. See PDF.] at T = 300 K (Equation (22)). Plot (b) compares the averaged lifetimes [Forumla omitted. See PDF.] defined in Equation (16) to the inverse averaged squared phonon-phonon interaction strengths [Forumla omitted. See PDF.]. The line colours used in Plot (a) match the marker colours in Plot (b). In Plot (b) the dashed black line is a linear fit to the data for d-Si and the oC24, K-II/C-I, K-V/C-VI and K-VII/C-V framework structures.

Figure 10. Anisotropy in the lattice thermal conductivity of the oC24 (a), K-V/C-VI (b) and K-VII/C-V (c) structures. For the orthorhombic structures shown in (a,b) the three diagonal components of the [Forumla omitted. See PDF.] tensor, [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] are shown together with the average [Forumla omitted. See PDF.] defined in Equation (15). For the hexagonal K-VII/C-V structure shown in (c) the [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] are equivalent, and [Forumla omitted. See PDF.] is shown together with [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Table 1

Summary of the technical parameters used in calculations on the six structures examined in this work: unit-cell definition and $k$-point sampling used for geometry optimisation and to calculate the infrared (IR) and Raman activities; cell definition, supercells and $k$-point sampling used to compute the second- and third-order force constants ${\mathbf{\Phi }}^{\left(2\right)}$/${\mathbf{\Phi }}^{\left(3\right)}$; and $\mathit{q}$-point meshes used to compute the phonon density of states (DoS) and vibrational Helmholtz free energy $A_{\mathrm{vib}}$ and to evaluate the lattice thermal conductivity ${\kappa }_{\mathrm{latt}}$. Cell definitions are either the primitive (“Prim.”) or conventional (“Conv.”) cells, for systems where these differ. $k$- and $\mathit{q}$-point meshes are given as regular $\mathrm{\Gamma }$-centered meshes generated using the Monkhorst-Pack scheme [47].

		d-Si	C-II	oC24	K-II/C-I	K-V/C-VI	K-VII/C-V
Optimisation	Cell	Conv.	Prim.	Prim.	-	Prim.	-
	$k$-points	6 × 6 × 6	4 × 4 × 4	8 × 8 × 2	3 × 3 × 3	4 × 4 × 3	4 × 4 × 2
$I_{\mathrm{IR}}$/$I_{\mathrm{Raman}}$	Cell	Prim.	Prim.	Prim.	-	Prim.	-
	$k$-points	10 × 10 × 10	4 × 4 × 4	8 × 8 × 2	3 × 3 × 3	4 × 4 × 3	4 × 4 × 2
${\mathbf{\Phi }}^{\left(2\right)}$	Cell	Conv.	Conv.	Conv.	-	Conv.	-
	Supercell (# atoms) 1	3 × 3 × 3 (216)	$\left(\begin{array}{ccc}0& 1& 1\\ 1& 0& 1\\ 1& 1& 0\end{array}\right)$ (272)	$\left(\begin{array}{ccc}2& 2& 0\\ -2& 2& 0\\ 0& 0& 2\end{array}\right)$ (384)	2 × 2 × 2 (368)	$\left(\begin{array}{ccc}1& 1& 0\\ -1& 1& 0\\ 0& 0& 2\end{array}\right)$ (320)	2 × 2 × 1 (272)
	$k$-points	2 × 2 × 2	2 × 2 × 2	2 × 2 × 1	2 × 2 × 2	2 × 2 × 2	2 × 2 × 2
${\mathbf{\Phi }}^{\left(3\right)}$	Cell	Conv.	Conv.	Conv.	-	Conv.	-
	Supercell (# atoms)	2 × 2 × 2 (64)	1 × 1 × 1 (136)	3 × 1 × 1 (72)	1 × 1 × 1 (46)	1 × 1 × 1 (80)	1 × 1 × 1 (68)
	$k$-points	3 × 3 × 3	2 × 2 × 2	3 × 2 × 2	3 × 3 × 3	4 × 2 × 3	4 × 4 × 2
$\mathit{q}$-points	DoS/$A_{\mathrm{vib}}$	48 × 48 × 48	24 × 24 × 24	40 × 40 × 16	24 × 24 × 24	24 × 24 × 24	24 × 24 × 16
	${\kappa }_{\mathrm{latt}}$ 2	36 × 36 × 36	16 × 16 × 16	13 × 13 × 5	13 × 13 × 13	9 × 9 × 9	14 × 14 × 9

¹ The non-diagonal supercell matrices used for C-II, oC24 and K-V/C-VI are equivalent to 2 × 2 × 2, 4 × 4 × 2 and 2 × 2 × 2 expansions of the respective primitive cells. ² Convergence of the ${\kappa }_{\mathrm{latt}}$ of C-II with respect to the $\mathit{q}$-point sampling density is slow, such that a coarser 8 × 8 × 8 mesh yields a ${\kappa }_{\mathrm{latt}}$ within ~5% of the converged value obtained with the 24 × 24 × 24 mesh used to compute the density of states and vibrational free energy.

Table 2

List of the six structures examined in this study with the sources of the initial structures (mp—Materials Project; ICSD—Inorganic Crystal Structure Database), the crystallographic spacegroups, and the number of atoms, lattice constants (a, b, c) and volumes V of the conventional unit cells.

	Source	Spacegroup	${\mathit{n}}_{\mathbf{a}}$	a (Å)	b (Å)	c (Å)	V (Å3)
d-Si	mp-149	$Fd\overline{3}m$	8	5.436	-	-	160.6
C-II	mp-16220 [32]	$Fd\overline{3}m$	136	14.629	-	-	3130.9
oC24	mp-1095269 [36]	$Cmcm$	24	3.826	10.697	12.660	518.1
K-II/C-I	mp-971662 [35]	$Pm\overline{3}m$	46	10.227	-	-	1052.5
K-V/C-VI	ICSD: 189396 [35]	$Cmmm$	80	10.158	17.596	10.325	1845.6
K-VII/C-V	mp-1203790 [35]	$P{6}_3/mmc$	68	10.358	-	16.956	1575.3

Table 3

Analysis of the $3N$ $\mathrm{\Gamma }$-point phonon modes in the six structures examined in this work. For each structure we list the irreducible representation ${\mathrm{\Gamma }}_{3N}$ spanned by the modes together with the expected numbers of infrared (IR)- and Raman-active modes based on symmetry.

	$3\mathit{N}$	${\mathbf{\Gamma }}_{3\mathit{N}}$	# IR Active	# Raman Active
d-Si	6	T1u + T2g	1	1
C-II	102	3 A1g + A1u + A2g + 3 A2u + 4 Eg + 4 Eu + 5 T1g + 8 T1u + 8 T2g + 5 T2u	8	15
oC24	36	6 Ag + 3 Au + 3 B1g + 6 B1u + 3 B2g + 6 B2u + 6 B3g + 3 B3u	15	18
K-II/C-I	138	3 A1g + 2 A1u + 4 A2g + 2 A2u + 7 Eg + 4 Eu + 8 T1g + 10 T1u + 8 T2g + 9 T2u	10	18
K-V/C-VI	120	18 Ag + 11 Au + 15 B1g + 16 B1u + 13 B2g + 17 B2u + 14 B3g + 16 B3u	49	50
K-VII/C-V	204	11 A1g + 5 A1u + 7 A2g + 11 A2u + 5 B1g + 11 B1u + 11 B2g + 7 B2u + 16 E1g + 18 E1u + 18 E2g + 16 E2u	29	45

Table 4

Relative energies of the five framework structures examined in this study compared to bulk diamond Si. For each system we show the differences in the lattice internal energies $U_{\mathrm{latt}}$ together with the differences in the Helmholtz energies $A=U-TS$ at T = 300 K and the corresponding differences in the internal energy $U=U_{\mathrm{latt}}+U_{\mathrm{vib}}$ and entropy $-TS$.

	$\mathbf{\Delta }$ (kJ mol−1 atom−1)
	$\mathbf{\Delta }{\mathit{U}}_{\mathrm{latt}}$	$\mathbf{\Delta }\mathit{A}$	$\mathbf{\Delta }\mathit{U}$	$\mathbf{-}\mathit{T}\mathbf{\Delta }\mathit{S}$
d-Si	-	-	-	-
C-II	32.62	31.47	32.39	−0.92
oC24	8.97	8.73	8.91	−0.18
K-II/C-I	7.47	7.25	7.39	−0.13
K-V/C-VI	9.14	8.88	2.70	−0.17
K-VII/C-V	6.82	6.62	6.74	−0.12

Table 5

Calculated lattice thermal conductivity of the six structures examined in this work at T = 300 K. For each structure the unique diagonal components ${\kappa }_{xx}$, ${\kappa }_{yy}$ and ${\kappa }_{zz}$ of the ${\kappa }_{\mathrm{latt}}$ tensor, corresponding to transport along the Cartesian x, y and z directions, are shown together with the average ${\kappa }_{\mathrm{ave}}$ defined in Equation (15). Decompositions of the ${\kappa }_{\mathrm{ave}}$ into a harmonic term $\kappa /{\tau }^{\mathrm{CRTA}}$ and lifetime term ${\tau }^{\mathrm{CRTA}}$ as defined in Equation (16) are also shown, together with the scaled averaged three-phonon iteraction strengths $\tilde{P}$ defined in Equation (18).

	$\mathit{\kappa }$ (W m−1 K−1)				$\mathit{\kappa }/{\mathit{\tau }}^{\mathbf{CTRA}}$ (W m−1 K−1 ps−1)	${\mathit{\tau }}^{\mathbf{CRTA}}$ (ps)	$\tilde{\mathit{P}}\times {\left(3{\mathit{n}}_{\mathbf{a}}\right)}^2$ (10−8 eV2)
	${\mathbf{\kappa }}_{\mathbf{xx}}$	${\mathbf{\kappa }}_{\mathbf{yy}}$	${\mathbf{\kappa }}_{\mathbf{zz}}$	${\mathbf{\kappa }}_{\mathbf{ave}}$
d-Si	136.24	-	-	136.24	5.002	27.24	2.164
C-II	6.33	-	-	6.33	0.458	13.81	2.704
oC24	57.86	44.70	20.21	40.92	2.295	17.83	2.983
K-II/C-I	43.54	-	-	43.54	0.829	52.52	1.184
K-V/C-VI	38.24	37.68	32.96	36.29	0.815	44.53	1.327
K-VII/C-V	30.45	-	32.56	31.16	0.770	40.45	1.477

Supplementary Materials

The following are available online at https://www.mdpi.com/article/10.3390/molecules27196431/s1, additional calculations on the C-II structure and comparison to previous modelling studies; complete peak tables for the $\mathrm{\Gamma }$-point vibrational modes of the six systems examined in this work including frequencies, irreducible representations, intensities, and room-temperature linewidths; and analysis of the room-temperature phonon linewidths/lifetimes in each of the six structures including unscaled w-JDoS functions ${\overline{N}}_2^{\left(1\right)}(\omega ,T)$/${\overline{N}}_2^{\left(2\right)}(\omega ,T)$ and linear fits to obtain averaged three-phonon interaction strengths $\tilde{P}$. Refs [70,71,72,73,74] are cited in the supplementary~materials.

References

1. Tan, G.; Zhao, L.D.; Kanatzidis, M.G. Rationally Designing High-Performance Bulk Thermoelectric Materials. Chem. Rev.; 2016; 116, pp. 12123-12149. [DOI: https://dx.doi.org/10.1021/acs.chemrev.6b00255] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27580481]

2. Firth, A.; Zhang, B.; Yang, A. Quantification of global waste heat and its environmental effects. Appl. Energy; 2019; 235, pp. 1314-1334. [DOI: https://dx.doi.org/10.1016/j.apenergy.2018.10.102]

3. Freer, R.; Powell, A.V. Realising the potential of thermoelectric technology: A Roadmap. J. Mater. Chem. C; 2020; 8, pp. 441-463. [DOI: https://dx.doi.org/10.1039/C9TC05710B]

4. Yan, Q.; Kanatzidis, M.G. High-performance thermoelectrics and challenges for practical devices. Nat. Mater.; 2022; 21, pp. 503-513. [DOI: https://dx.doi.org/10.1038/s41563-021-01109-w] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34675376]

5. Zeier, W.G.; Zevalkink, A.; Gibbs, Z.M.; Hautier, G.; Kanatzidis, M.G.; Snyder, G.J. Thinking Like a Chemist: Intuition in Thermoelectric Materials. Angew. Chem. Int. Ed.; 2016; 55, pp. 6826-6841. [DOI: https://dx.doi.org/10.1002/anie.201508381] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27111867]

6. Li, C.W.; Hong, J.; May, A.F.; Bansal, D.; Chi, S.; Hong, T.; Ehlers, G.; Delaire, O. Orbitally driven giant phonon anharmonicity in SnSe. Nat. Phys.; 2015; 11, pp. 1063-1069. [DOI: https://dx.doi.org/10.1038/nphys3492]

7. Aseginolaza, U.; Bianco, R.; Monacelli, L.; Paulatto, L.; Calandra, M.; Mauri, F.; Bergara, A.; Errea, I. Phonon Collapse and Second-Order Phase Transition in Thermoelectric SnSe. Phys. Rev. Lett.; 2019; 122, 075901. [DOI: https://dx.doi.org/10.1103/PhysRevLett.122.075901]

8. Aseginolaza, U.; Bianco, R.; Monacelli, L.; Paulatto, L.; Calandra, M.; Mauri, F.; Bergara, A.; Errea, I. Strong anharmonicity and high thermoelectric efficiency in high-temperature SnS from first principles. Phys. Rev. B; 2019; 100, 214307. [DOI: https://dx.doi.org/10.1103/PhysRevB.100.214307]

9. Nielsen, M.D.; Ozolins, V.; Heremans, J.P. Lone pair electrons minimize lattice thermal conductivity. Energy Environ. Sci.; 2013; 6, pp. 570-578. [DOI: https://dx.doi.org/10.1039/C2EE23391F]

10. Biswas, K.; He, J.; Blum, I.D.; Wu, C.I.; Hogan, T.P.; Seidman, D.N.; Dravid, V.P.; Kanatzidis, M.G. High-performance bulk thermoelectrics with all-scale hierarchical architectures. Nature; 2012; 489, pp. 414-418. [DOI: https://dx.doi.org/10.1038/nature11439]

11. Gonçalves, A.P.; Godart, C. New promising bulk thermoelectrics: Intermetallics, pnictides and chalcogenides. Eur. Phys. J. B; 2014; 87, 42. [DOI: https://dx.doi.org/10.1140/epjb/e2014-40989-3]

12. Dolyniuk, J.A.; Owens-Baird, B.; Wang, J.; Zaikina, J.V.; Kovnir, K. Clathrate thermoelectrics. Mater. Sci. Eng. R Rep.; 2016; 108, pp. 1-46. [DOI: https://dx.doi.org/10.1016/j.mser.2016.08.001]

13. Xia, K.; Hu, C.; Fu, C.; Zhao, X.; Zhu, T. Half-Heusler thermoelectric materials. Appl. Phys. Lett.; 2021; 118, 140503. [DOI: https://dx.doi.org/10.1063/5.0043552]

14. Pei, Y.; Wang, H.; Snyder, G.J. Band Engineering of Thermoelectric Materials. Adv. Mater.; 2012; 24, pp. 6125-6135. [DOI: https://dx.doi.org/10.1002/adma.201202919]

15. Tan, X.; Wang, H.; Liu, G.; Noudem, J.G.; Hu, H.; Xu, J.; Shao, H.; Jiang, J. Designing band engineering for thermoelectrics starting from the periodic table of elements. Mater. Today Phys.; 2018; 7, pp. 35-44. [DOI: https://dx.doi.org/10.1016/j.mtphys.2018.10.004]

16. Liu, W.; Lukas, K.C.; McEnaney, K.; Lee, S.; Zhang, Q.; Opeil, C.P.; Chen, G.; Ren, Z. Studies on the Bi2Te3–Bi2Se3–Bi2S3 system for mid-temperature thermoelectric energy conversion. Energy Environ. Sci.; 2013; 6, pp. 552-560. [DOI: https://dx.doi.org/10.1039/C2EE23549H]

17. Xie, H.; Su, X.; Hao, S.; Zhang, C.; Zhang, Z.; Liu, W.; Yan, Y.; Wolverton, C.; Tang, X.; Kanatzidis, M.G. Large Thermal Conductivity Drops in the Diamondoid Lattice of CuFeS₂ by Discordant Atom Doping. J. Am. Chem. Soc.; 2019; 141, pp. 18900-18909. [DOI: https://dx.doi.org/10.1021/jacs.9b10983]

18. Skelton, J.M. Approximate models for the lattice thermal conductivity of alloy thermoelectrics. J. Mater. Chem. C; 2021; 9, pp. 11772-11787. [DOI: https://dx.doi.org/10.1039/D1TC02026A]

19. Lu, Y.; Zheng, F.w.; Yang, Y.; Zhang, P.; Zhang, D.B. Phase stabilities of Cmcm and Pnma SnSe studied by phonon quasiparticle approach. Phys. Rev. B; 2019; 100, 054304. [DOI: https://dx.doi.org/10.1103/PhysRevB.100.054304]

20. Slack, G.A. New Materials and Performance Limits for Thermoelectric Cooling. CRC Handbook of Thermoelectrics; Rowe, D.M. CRC Press: Boca Raton, FL, USA, 1985; 407.

21. Toberer, E.S.; Zevalkink, A.; Snyder, G.J. Phonon engineering through crystal chemistry. J. Mater. Chem.; 2011; 21, pp. 15843-15852. [DOI: https://dx.doi.org/10.1039/c1jm11754h]

22. Kawaharada, Y.; Kurosaki, K.; Uno, M.; Yamanaka, S. Thermoelectric properties of CoSb₃. J. Alloy. Compd.; 2001; 315, pp. 193-197. [DOI: https://dx.doi.org/10.1016/S0925-8388(00)01275-5]

23. Pei, Y.Z.; Yang, J.; Chen, L.D.; Zhang, W.; Salvador, J.R.; Yang, J. Improving thermoelectric performance of caged compounds through light-element filling. Appl. Phys. Lett.; 2009; 95, 042101. [DOI: https://dx.doi.org/10.1063/1.3182800]

24. Shi, X.; Yang, J.; Salvador, J.R.; Chi, M.; Cho, J.Y.; Wang, H.; Bai, S.; Yang, J.; Zhang, W.; Chen, L. Multiple-Filled Skutterudites: High Thermoelectric Figure of Merit through Separately Optimizing Electrical and Thermal Transports. J. Am. Chem. Soc.; 2011; 133, pp. 7837-7846. [DOI: https://dx.doi.org/10.1021/ja111199y]

25. Nolas, G.S.; Slack, G.A.; Schujman, S.B. Semiconductor clathrates: A phonon glass electron crystal material with potential for thermoelectric applications. Recent Trends in Thermoelectric Materials Research I; Tritt, T.M. Academic Press: London, UK, 2001; Volume 69, pp. 255-300. [DOI: https://dx.doi.org/10.1016/S0080-8784(01)80152-6]

26. Nolas, G.S.; Cohn, J.L.; Slack, G.A.; Schujman, S.B. Semiconducting Ge clathrates: Promising candidates for thermoelectric applications. Appl. Phys. Lett.; 1998; 73, pp. 178-180. [DOI: https://dx.doi.org/10.1063/1.121747]

27. Cohn, J.L.; Nolas, G.S.; Fessatidis, V.; Metcalf, T.H.; Slack, G.A. Glasslike Heat Conduction in High-Mobility Crystalline Semiconductors. Phys. Rev. Lett.; 1999; 82, pp. 779-782. [DOI: https://dx.doi.org/10.1103/PhysRevLett.82.779]

28. Zhu, G.H.; Lee, H.; Lan, Y.C.; Wang, X.W.; Joshi, G.; Wang, D.Z.; Yang, J.; Vashaee, D.; Guilbert, H.; Pillitteri, A. et al. Increased Phonon Scattering by Nanograins and Point Defects in Nanostructured Silicon with a Low Concentration of Germanium. Phys. Rev. Lett.; 2009; 102, 196803. [DOI: https://dx.doi.org/10.1103/PhysRevLett.102.196803] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19518985]

29. He, R.; Heyn, W.; Thiel, F.; Pérez, N.; Damm, C.; Pohl, D.; Rellinghaus, B.; Reimann, C.; Beier, M.; Friedrich, J. et al. Thermoelectric properties of silicon and recycled silicon sawing waste. J. Mater.; 2019; 5, pp. 15-33. [DOI: https://dx.doi.org/10.1016/j.jmat.2018.11.004]

30. Shanks, H.R.; Maycock, P.D.; Sidles, P.H.; Danielson, G.C. Thermal Conductivity of Silicon from 300 to 1400°K. Phys. Rev.; 1963; 130, pp. 1743-1748. [DOI: https://dx.doi.org/10.1103/PhysRev.130.1743]

31. Glassbrenner, C.J.; Slack, G.A. Thermal Conductivity of Silicon and Germanium from 3°K to the Melting Point. Phys. Rev.; 1964; 134, pp. A1058-A1069. [DOI: https://dx.doi.org/10.1103/PhysRev.134.A1058]

32. Dong, J.; Sankey, O.F.; Kern, G. Theoretical study of the vibrational modes and their pressure dependence in the pure clathrate-II silicon framework. Phys. Rev. B; 1999; 60, pp. 950-958. [DOI: https://dx.doi.org/10.1103/PhysRevB.60.950]

33. Gryko, J.; McMillan, P.F.; Marzke, R.F.; Ramachandran, G.K.; Patton, D.; Deb, S.K.; Sankey, O.F. Low-density framework form of crystalline silicon with a wide optical band gap. Phys. Rev. B; 2000; 62, pp. R7707-R7710. [DOI: https://dx.doi.org/10.1103/PhysRevB.62.R7707]

34. Karttunen, A.J.; Fässler, T.F.; Linnolahti, M.; Pakkanen, T.A. Structural Principles of Semiconducting Group 14 Clathrate Frameworks. Inorg. Chem.; 2011; 50, pp. 1733-1742. [DOI: https://dx.doi.org/10.1021/ic102178d]

35. Zhao, H.Y.; Wang, J.; Ma, Q.M.; Liu, Y. sp₃-Bonded silicon allotropes based on the Kelvin problem. Phys. Chem. Chem. Phys.; 2013; 15, pp. 17619-17625. [DOI: https://dx.doi.org/10.1039/c3cp50946j]

36. Kim, D.Y.; Stefanoski, S.; Kurakevych, O.O.; Strobel, T.A. Synthesis of an open-framework allotrope of silicon. Nat. Mater.; 2015; 14, pp. 169-173. [DOI: https://dx.doi.org/10.1038/nmat4140]

37. Nolas, G.S.; Beekman, M.; Gryko, J.; Lamberton, G.A.; Tritt, T.M.; McMillan, P.F. Thermal conductivity of elemental crystalline silicon clathrate Si₁₃₆. Appl. Phys. Lett.; 2003; 82, pp. 910-912. [DOI: https://dx.doi.org/10.1063/1.1544640]

38. Skelton, J.M.; Burton, L.A.; Parker, S.C.; Walsh, A.; Kim, C.E.; Soon, A.; Buckeridge, J.; Sokol, A.A.; Catlow, C.R.A.; Togo, A. et al. Anharmonicity in the High-Temperature Cmcm Phase of SnSe: Soft Modes and Three-Phonon Interactions. Phys. Rev. Lett.; 2016; 117, 075502. [DOI: https://dx.doi.org/10.1103/PhysRevLett.117.075502]

39. Tang, J.; Skelton, J.M. Impact of noble-gas filler atoms on the lattice thermal conductivity of CoSb₃ skutterudites: First-principles modelling. J. Phys. Condens. Matter; 2021; 33, 164002. [DOI: https://dx.doi.org/10.1088/1361-648X/abd8b8]

40. Cen, J.; Pallikara, I.; Skelton, J.M. Structural Dynamics and Thermal Transport in Bismuth Chalcogenide Alloys. Chem. Mater.; 2021; 33, pp. 8404-8417. [DOI: https://dx.doi.org/10.1021/acs.chemmater.1c02777]

41. Kresse, G.; Hafner, J. Ab Initio Mol. Dyn. Liq. Met. Phys. Rev. B; 1993; 47, pp. 558-561. [DOI: https://dx.doi.org/10.1103/PhysRevB.47.558]

42. Jain, A.; Ong, S.P.; Hautier, G.; Chen, W.; Richards, W.D.; Dacek, S.; Cholia, S.; Gunter, D.; Skinner, D.; Ceder, G. et al. The Materials Project: A materials genome approach to accelerating materials innovation. APL Mater.; 2013; 1, 011002. [DOI: https://dx.doi.org/10.1063/1.4812323]

43. Perdew, J.P.; Ruzsinszky, A.; Csonka, G.I.; Vydrov, O.A.; Scuseria, G.E.; Constantin, L.A.; Zhou, X.; Burke, K. Restoring the Density-Gradient Expansion for Exchange in Solids and Surfaces. Phys. Rev. Lett.; 2008; 100, 136406. [DOI: https://dx.doi.org/10.1103/PhysRevLett.100.136406]

44. Skelton, J.M.; Tiana, D.; Parker, S.C.; Togo, A.; Tanaka, I.; Walsh, A. Influence of the exchange-correlation functional on the quasi-harmonic lattice dynamics of II-VI semiconductors. J. Chem. Phys.; 2015; 143, 064710. [DOI: https://dx.doi.org/10.1063/1.4928058]

45. Blöchl, P.E. Projector augmented-wave method. Phys. Rev. B; 1994; 50, pp. 17953-17979. [DOI: https://dx.doi.org/10.1103/PhysRevB.50.17953]

46. Kresse, G.; Joubert, D. From ultrasoft pseudopotentials to the projector augmented-wave method. Phys. Rev. B; 1999; 59, pp. 1758-1775. [DOI: https://dx.doi.org/10.1103/PhysRevB.59.1758]

47. Monkhorst, H.J.; Pack, J.D. Special points for Brillouin-zone integrations. Phys. Rev. B; 1976; 13, pp. 5188-5192. [DOI: https://dx.doi.org/10.1103/PhysRevB.13.5188]

48. Togo, A.; Tanaka, I. Spglib: A software library for crystal symmetry search. arXiv; 2018; arXiv: cond-mat.mtrl-sci/1808.01590

49. Togo, A.; Tanaka, I. First principles phonon calculations in materials science. Scr. Mater.; 2015; 108, pp. 1-5. [DOI: https://dx.doi.org/10.1016/j.scriptamat.2015.07.021]

50. Togo, A.; Chaput, L.; Tanaka, I. Distributions of phonon lifetimes in Brillouin zones. Phys. Rev. B; 2015; 91, 094306. [DOI: https://dx.doi.org/10.1103/PhysRevB.91.094306]

51. Porezag, D.; Pederson, M.R. Infrared intensities and Raman-scattering activities within density-functional theory. Phys. Rev. B; 1996; 54, pp. 7830-7836. [DOI: https://dx.doi.org/10.1103/PhysRevB.54.7830] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/9984457]

52. Skelton, J.M.; Burton, L.A.; Jackson, A.J.; Oba, F.; Parker, S.C.; Walsh, A. Lattice dynamics of the tin sulphides SnS₂, SnS and Sn₂S₃: Vibrational spectra and thermal transport. Phys. Chem. Chem. Phys.; 2017; 19, pp. 12452-12465. [DOI: https://dx.doi.org/10.1039/C7CP01680H] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28470289]

53. Gajdoš, M.; Hummer, K.; Kresse, G.; Furthmüller, J.; Bechstedt, F. Linear optical properties in the projector-augmented wave methodology. Phys. Rev. B; 2006; 73, 045112. [DOI: https://dx.doi.org/10.1103/PhysRevB.73.045112]

54. Souza, I.; Íñiguez, J.; Vanderbilt, D. First-Principles Approach to Insulators in Finite Electric Fields. Phys. Rev. Lett.; 2002; 89, 117602. [DOI: https://dx.doi.org/10.1103/PhysRevLett.89.117602]

55. Kurakevych, O.O.; Strobel, T.A.; Kim, D.Y.; Muramatsu, T.; Struzhkin, V.V. Na-Si Clathrates Are High-Pressure Phases: A Melt-Based Route to Control Stoichiometry and Properties. Cryst. Growth Des.; 2013; 13, pp. 303-307. [DOI: https://dx.doi.org/10.1021/cg3017084]

56. Jennings, H.M.; Richman, M.H. A Hexagonal (Wurtzite) Form of Silicon. Science; 1976; 193, pp. 1242-1243. [DOI: https://dx.doi.org/10.1126/science.193.4259.1242]

57. Tang, J.; Maurice, J.L.; Fossard, F.; Florea, I.; Chen, W.; Johnson, E.V.; Foldyna, M.; Yu, L.; Roca i Cabarrocas, P. Natural occurrence of the diamond hexagonal structure in silicon nanowires grown by a plasma-assisted vapour–liquid–solid method. Nanoscale; 2017; 9, pp. 8113-8118. [DOI: https://dx.doi.org/10.1039/C7NR01299C]

58. Kasper, J.S.; Richards, S.M. The crystal structures of new forms of silicon and germanium. Acta Crystallogr.; 1964; 17, pp. 752-755. [DOI: https://dx.doi.org/10.1107/S0365110X64001840]

59. Duclos, S.J.; Vohra, Y.K.; Ruoff, A.L. Experimental study of the crystal stability and equation of state of Si to 248 GPa. Phys. Rev. B; 1990; 41, pp. 12021-12028. [DOI: https://dx.doi.org/10.1103/PhysRevB.41.12021]

60. Crain, J.; Ackland, G.J.; Maclean, J.R.; Piltz, R.O.; Hatton, P.D.; Pawley, G.S. Reversible pressure-induced structural transitions between metastable phases of silicon. Phys. Rev. B; 1994; 50, pp. 13043-13046. [DOI: https://dx.doi.org/10.1103/PhysRevB.50.13043]

61. McMahon, M.I.; Nelmes, R.J.; Wright, N.G.; Allan, D.R. Pressure dependence of the Imma Phase Silicon. Phys. Rev. B; 1994; 50, pp. 739-743. [DOI: https://dx.doi.org/10.1103/PhysRevB.50.739]

62. Hanfland, M.; Schwarz, U.; Syassen, K.; Takemura, K. Crystal Structure of the High-Pressure Phase Silicon VI. Phys. Rev. Lett.; 1999; 82, pp. 1197-1200. [DOI: https://dx.doi.org/10.1103/PhysRevLett.82.1197]

63. Okada, Y.; Tokumaru, Y. Precise determination of lattice parameter and thermal expansion coefficient of silicon between 300 and 1500 K. J. Appl. Phys.; 1984; 56, pp. 314-320. [DOI: https://dx.doi.org/10.1063/1.333965]

64. Kasper, J.S.; Hagenmuller, P.; Pouchard, M.; Cros, C. Clathrate Structure of Silicon Na₈Si₄₆ and Na_xSi₁₃₆ (x < 11). Science; 1965; 150, pp. 1713-1714. [DOI: https://dx.doi.org/10.1126/science.150.3704.1713]

65. Momma, K.; Izumi, F. VESTA3 Three-Dimens. Vis. Crystal, Vol. Morphol. Data. J. Appl. Crystallogr.; 2011; 44, pp. 1272-1276. [DOI: https://dx.doi.org/10.1107/S0021889811038970]

66. Parker, J.H.; Feldman, D.W.; Ashkin, M. Raman Scattering by Silicon and Germanium. Phys. Rev.; 1967; 155, pp. 712-714. [DOI: https://dx.doi.org/10.1103/PhysRev.155.712]

67. Beekman, M.; Nolas, G. Synthesis and thermal conductivity of type II silicon clathrates. Phys. B Condens. Matter; 2006; 383, pp. 111-114. [DOI: https://dx.doi.org/10.1016/j.physb.2006.03.070]

68. Härkönen, V.J.; Karttunen, A.J. Ab initio studies on the lattice thermal conductivity of silicon clathrate frameworks II and VIII. Phys. Rev. B; 2016; 93, 024307. [DOI: https://dx.doi.org/10.1103/PhysRevB.93.024307]

69. Norouzzadeh, P.; Krasinski, J.S.; Tadano, T. Thermal conductivity of type-I, type-II, and type-VIII pristine silicon clathrates: A first-principles study. Phys. Rev. B; 2017; 96, 245201. [DOI: https://dx.doi.org/10.1103/PhysRevB.96.245201]

70. Ceperley, D.M.; Alder, B.J. Ground State of the Electron Gas by a Stochastic Method. Phys. Rev. Lett.; 1980; 45, pp. 566-569. [DOI: https://dx.doi.org/10.1103/PhysRevLett.45.566]

71. Perdew, J.P.; Chevary, J.A.; Vosko, S.H.; Jackson, K.A.; Pederson, M.R.; Singh, D.J.; Fiolhais, C. Atoms, molecules, solids, and 117 surfaces: Applications of the generalized gradient approximation for exchange and correlation. Phys. Rev. B; 1992; 46, pp. 6671-6687. [DOI: https://dx.doi.org/10.1103/PhysRevB.46.6671]

72. Perdew, J.P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett.; 1996; 77, pp. 3865-3868. [DOI: https://dx.doi.org/10.1103/PhysRevLett.77.3865]

73. Murnaghan, F.D. The Compressibility of Media under Extreme Pressures. Proc. Natl. Acad. Sci. USA; 1944; 30, pp. 244-247. [DOI: https://dx.doi.org/10.1073/pnas.30.9.244]

74. Birch, F. Finite Elastic Strain of Cubic Crystals. Phys. Rev.; 1947; 71, pp. 809-824. [DOI: https://dx.doi.org/10.1103/PhysRev.71.809]