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

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 κ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 κ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.


Introduction
Among the most important challenges to contemporary science is meeting the everincreasing 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 energyintensive 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].
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 κ 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 κ latt .

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 Γ-centered Monkhorst-Pack k-point meshes [47] listed in Table 1 and a Gaussian smearing with a width σ 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(ν) were computed by interpolating the phonon frequencies onto the regular Γ-centered q-point grids listed in Table 1. The DoS curves were broadened using a Gaussian smearing with a width σ = 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 vib . Phonon dispersion curves ν(q) were computed by interpolating frequencies along strings of 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 q-point grids listed in Table 1.
Finally, the infrared (IR) and Raman activities of the Γ-point phonon modes were computed using the procedure outlined in Refs [51,52]. The Born effective-charge tensors Z * required for the IR activities were calculated using the density-functional perturbation theory (DFPT) routines in VASP [53]. The polarisability derivatives ∂α/∂Q ≡ ∂ε ∞ /∂Q were computed using numerical differentiation with a step size ∆Q chosen such that the largest Cartesian displacement was 10 −2 Å. The ε ∞ 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 . 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 Φ (2) /Φ (3) ; and q-point meshes used to compute the phonon density of states (DoS) and vibrational Helmholtz free energy A vib and to evaluate the lattice thermal conductivity κ latt . Cell definitions are either the primitive ("Prim.") or conventional ("Conv.") cells, for systems where these differ. k-and q-point meshes are given as regular Γ-centered meshes generated using the Monkhorst-Pack scheme [47]. 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. 2 Convergence of the κ latt of C-II with respect to the q-point sampling density is slow, such that a coarser 8 × 8 × 8 mesh yields a κ 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.

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 Si VII K structure reported as the Type-V clathrate in Ref. [34] (mp-1203790), the Si II K structure identified as corresponding to the empty Type-I clathrate (mp-971662), and the Si V K 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]. 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.

Phonon Spectra
The calculated phonon dispersion ν(q) and density of states (DoS) curves g(ν) of the six structures are compared in Figure 2. The primitive cell of bulk diamond Si has a simple unit cell containing n a = 2 atoms, resulting in 3n a = 6 branches at each wavevector q in the dispersion, and the high symmetry of the Fd3m 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 3n a = 102 (C-II), 36 (oC24), 138 (K-II/C-I), 120 (K-V/C-VI) and 204 branches (K-VII/C-V). (c) the "oC24" structure reported in Ref. [36]; (d) the Type-I clathrate reported as Si II K in Ref. [35] (K-II/C-I); (e) the Si V K 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 Si VII K in Ref. [34] (K-VII/C-V). These images were prepared using the VESTA software [65].
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 q = Γ 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] where the sum runs over the 3n a modes j at q = Γ with ordinal frequencies ν j , spectroscopic activities I j and linewidths Γ j , and we have dropped the 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 Γ-point mode eigenvector W j the Cartesian displacement X j,k of the kth atom is given by: where m k is the atomic mass. The infrared activities I IR are given by the change in polarisation along the mode eigenvector and are calculated from the Born effective-charge tensors of the atoms Z * k using the expression: [51] I IR where α and β are the Cartesian directions x, y and z. The Raman activity I Raman depends on the change in polarisability α, which can be computed from the derivative of the high-frequency dielectric constant ε ∞ with respect to the mode amplitude Q: [51] I Raman 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 I Raman that would be measured from a powder sample: [51] where I αβ denote the components of the Raman tensor. In this work we compute the linewidths Γ 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 Γ 3N spanned by the Γ-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 Fm3m and Pm3m spacegroups adopted by d-Si, C-II and K-II/C-I both have the O 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 D 2h 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× 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 P6 3 /mmm spacegroup of the K-VII/C-V structure has the D 6h 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.

Thermodynamics
Within the harmonic approximation the temperature-dependent constant-volume (Helmholtz) free energy A(T) is given by: where U latt is the lattice internal energy, equal to the DFT total energy obtained after structural optimisation, and A vib = U vib − TS vib is the vibrational free energy from contributions of the phonons to the free energy and entropy. A vib can be computed directly from the vibrational partition function Z vib using the bridge relation: where N is the number of wavevectors q included in the partition function productequivalent to the number of unit cells in the crystal-so the normalisation yields the A vib per unit cell. The phonon partition function Z vib is computed from the (angular) phonon frequencies ω qj according to: where the product runs over phonon modes with wavevector q and band index j. The vibrational internal energy U vib can be computed according to: where n qj is the phonon occupation number given by the Bose-Einstein distribution. The phonon entropy S vib can then be obtained from A vib and U vib . Table 4 compares the calculated lattice energies and the 300 K A vib of the five framework structures to bulk diamond Si, while Figure 5 compares the free energy differences ∆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 ∆U 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 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 ∆A gives similar results. The differences in internal energy ∆U, including both lattice and phonon contributions, are comparable to the ∆U latt , indicating that differences in vibrational internal energy are much smaller than differences in U latt . The differences in the vibrational entropy −T∆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. 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 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 latt + U vib and entropy −TS. 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 ∆A with temperture to the differences in internal energy ∆U and entropy −T∆S ( Figure 6).
The ∆U all show slight increases with temperature, indicative of a larger U 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 ∆U between them remain roughly constant and the vibrational internal energy U 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 ∆U 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 8 Si 136 , 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 8 Si 46 , and NaSi 6 , 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.

Thermal Transport
Using the single-mode relaxation-time approximation (SM-RTA) model, the macroscopic thermal conductivity κ latt is computed as a sum of contributions from individual phonon modes as: Here N is the number of wavevectors in the summation and V is the unit-cell volume. The C qj are the modal heat capacities given by: The ν qj are the mode group velocities given by the derivative of the frequencies ω qj with respect to the wavevector: The τ qj are the phonon lifetimes, which are related to the phonon linewidths Γ qj by: As some of the clathrates have non-cubic spacegroups and therefore show anisotropic thermal transport, we also compute a scalar average κ ave according to: where κ xx , κ yy and κ zz are the three diagonal elements of the κ 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 κ 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 κ 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 κ 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 (Figures 2 and 3), the high thermal conductivity of d-Si can be explained by the wide dispersion of the acoustic modes and consequent large group velocities ν qj , 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 κ 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 κ latt is written as the product of a harmonic term and a weighted-average lifetime τ CRTA : where we have again omitted the temperature dependence for brevity. Comparing the harmonic term in the summand and the τ CRTA allows the difference in κ 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 τ 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 κ 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 κ 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 κ/τ CRTA at 2.3 W m −1 K −1 ps −1 (54% smaller), and the remaining framework structures all have κ/τ 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 κ/τ CRTA .
The variation in the lifetime term is more intricate. At 300 K, d-Si has a τ 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 κ 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 κ/τ CRTA .
The low thermal conductivity of the C-II structure, on the other hand, appears to be characterised by both a low κ/τ CRTA and a short τ 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 τ CRTA is also the shortest of the six compounds and is around 50% smaller than that of d-Si. 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 κ xx , κ yy and κ zz of the κ latt tensor, corresponding to transport along the Cartesian x, y and z directions, are shown together with the average κ ave defined in Equation (15). Decompositions of the κ ave into a harmonic term κ/τ CRTA and lifetime term τ CRTA as defined in Equation (16) are also shown, together with the scaled averaged three-phonon iteraction strengthsP defined in Equation (18). The phonon lifetimes are the inverse of the phonon linewidths Γ qj (c.f. Equation (14)). Within the SM-RTA model the Γ qj are calculated as a sum of contributions from threephonon scattering processes between triplets of modes qj, q j and q j according to: [50] Γ qj (T) = 18π The contribution from each triplet is determined by the three-phonon interaction strengths Φ −qj,q j ,q j , the occupation numbers n qj , 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 Φ −qj,q j ,q j [50]). The three-phonon interaction strengths are temperature independent and calculated from the harmonic frequencies and eigenvectors ω qj /W qj and the third-order (anharmonic) force constants Φ (3) , while the conservation of energy expression is calculated solely from the ω qj and is temperature dependant through the n qj . Based on the analysis in Ref. [50], we approximate the linewidths as: whereP is a weighted squared average three-phonon interaction strength and N 2 (q, ω, T) is a weighted two-phonon joint density of states (w-JDoS) defined as: with: For a given mode with a wavevector q and frequency ω these functions count respectively the number of possible energy-and momentum-conserving Class 1 and Class 2 events. As per Equation (11) the κ latt is proportional to the τ qj and therefore inversely proportional to the Γ qj , so that when the squared three-phonon interaction strengths are replaced by a constantP the κ latt is proportional to 1/P. A suitableP to reproduce the scalar average κ latt can therefore be obtained from a linear fit of the κ latt as a function of 1/P. As noted in Ref. [39], although the Φ −qj,q j ,q j in Equation (17) are temperature independent, replacing them with a single value cannot in general capture the temperature dependence of the κ latt , so theP are weakly temperature dependent and the fit must be performed at a chosen temperature T. Both theP and N 2 (q, ω, 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 (q, ω, T) it is prudent to average over q to obtain a function of frequency only, i.e.,: We note that both the N 2 (q, ω, T)/N 2 (ω, T) andP 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 (3n a ) 2 and theP scaled up by (3n a ) 2 .
The scaledN 2 (ω) 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, theN 2 (ω) 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.  (18) in the text. Plot (a) compares the scaled weighted two-phonon joint density of states (w-JDoS) functionsN 2 (ω) at T = 300 K (Equation (22)). Plot (b) compares the averaged lifetimes τ CRTA defined in Equation (16) to the inverse averaged squared phonon-phonon interaction strengths 1/P. 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.
A value ofP was fitted for each of the six structures at T = 300 K (see Supplementary Materials), and the scaled values are listed alongside the κ latt and CRTA components in Table 5. Excluding C-II, a plot of 1/P against the averaged lifetime τ CRTA (Figure 9b) shows a strong linear correlation across all the framework structures (R 2 = 0.997), providing clear evidence that differences in the τ CRTA can predominantly be attributed to differences in the phonon-phonon interaction strengths. TheP 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 largerN 2 (ω) (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 τ qj 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 κ 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 1 Si 136 and Na 8 Si 136 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 κ 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 κ 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 κ 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 κ xx , κ yy and κ zz components of the κ 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 κ latt components correspond to transport along the crystallographic a, b and c axes. The oC24 structure shows relatively large anisotropy. The largest κ 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 κ latt differ by a factor of two. The anisotropy in the K-V/C-VI structure is much smaller, with the κ latt components ranging from κ xx = 38.24 W m −1 K −1 to κ 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 κ 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 κ xx and κ yy components are equal and smaller than the κ zz component, corresponding to transport along the longer c axis, by around 2 W m −1 K −1 .  (15). For the hexagonal K-VII/C-V structure shown in (c) the κ xx and κ yy are equivalent, and κ xx = κ yy is shown together with κ zz and κ ave .

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 κ 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 κ latt . In this direction, the low κ 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.
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 Γ-point vibrational modes of the six systems examined in this work including frequencies, irreducible representations, intensities, and roomtemperature linewidths; and analysis of the room-temperature phonon linewidths/lifetimes in each of the six structures including unscaled w-JDoS functionsN (1) 2 (ω, T)/N (2) 2 (ω, T) and linear fits to obtain averaged three-phonon interaction strengthsP. Refs [70][71][72][73][74] are cited in the supplementary materials.

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.