First-Principles Study of the Electronic, Vibrational Properties and Anharmonic Effects of Some Si-Based Type-II Binary Clathrates

Electronic, vibrational, and anharmonic studies on some binary clathrate AxSi136 (A = Na, K, Rb, Cs; 0 < x ≤ 24) are theoretically presented. The Fermi energy lies in the range of 1.1 eV to 1.4 eV for NaxSi136 and increases as stoichiometry (x) is tuned from 8 to 12 to 16. The determined isotropic “Mexican-hat” shape of the guest-host potential describing Na motion in the Si28 cage indicates the “off-center” position when the temperature is elevated beyond zero. Accordingly, the calculated Na “off-center” displacements correlate well with the X-Ray Diffraction (XRD) data (0.4 Å–0.5 Å) for a similar composition range (0 < x < 24). The lack of first-principles analysis on quartic anharmonicity motivates us to initiate a self-consistent model to examine the temperature-dependent rattling frequency Ω(T) of the guest (Na, Rb). The predicted values of Ω(T) for Na24Si136 at 300 K are significantly higher (approximately six times larger) than the value at absolute zero, which contrasts with the case of Rb8Si136. Moreover, underestimation of the isotropic atomic displacement parameter Uiso is caused by the temperature-dependent quartic anharmonicity of Na, and this discrepancy might be offset by the square of the “off-center” displacement.


Introduction
Type-II inorganic clathrates have attracted considerable interest due to their unique crystal structures and promising thermodynamic applications [1][2][3][4][5][6]. These materials have a face-centered cubic (FCC) geometry and an open framework lattice with 136 atoms in the unit cell. Furthermore, the framework atoms in this structure are all in sp 3 bonding configurations, which form 20-and 28-atom polyhedron "cage" voids that are connected periodically in a 2:1 ratio. Each cage can trap an impurity ("guest") atom, which is usually from column I or column II of the periodic table. Several type-II clathrates have been synthesized with Si, Ge, or Sn on the framework lattice. The encapsulated guests are weakly bound to the host lattice, which means that their vibrational modes are expected to occur at very low frequencies. Several different type-II Si-based binary clathrates have been synthesized, including Na x Si 136 (x = 2.9, 5.1, 8.2, 14.7), K x Si 136 (x = 7.5, 17.8), and Na 16 (Rb,Cs) 8 Si 136 , and many of their properties have been measured [7][8][9][10][11]. One of the interesting results that has been obtained is the fact that (Na,K) 16 Rb 8 Si 136 [12] shows intermetallic behavior. This is caused by the charge transfer from alkaline metal atom guests to the framework conduction band states, which pushes the Fermi energy up into the lower portion of the conduction bands.
In this paper, we first report the results of a systematic, first-principles computational and theoretical study of the electronic properties of some of the binary type-II Si-based clathrates A x Si 136 (x = Na, K, Rb, Cs; 0 < x ≤ 24). Among the properties, we study the dependence of the pseudogap and Fermi level of the filled intermetallic clathrate Na x Si 136 (8 ≤ x ≤ 16) on the guest content (x). We also study the electronic density of states, which indicates the occurrence of temperature-dependent Knight shift behavior, and when possible, we compare our calculations with some experimental results for (Na,K) 16 Rb 8 Si 136 , which share similar characteristics [13].
One interesting and useful property of the type-II clathrate materials is that their novel cage-structured framework lattice is capable of housing impurity atoms that can form host-guest complexes. Of particular interest are the vibrational properties of these host-guest complexes, specifically the existence of low-frequency guest-associated "rattling" modes. Such modes may help to enhance the thermoelectric (TE) performance by minimizing the phonon contribution to the thermal conductivity. In this paper, we mainly report the results of the vibrational properties of several Si-based type-II clathrate materials. We especially focus on the vibrational modes and their anharmonic effects of various guests in the lattice cages.
If the guest atoms in large cages (28-atom cages) of type-II clathrates are small in comparison to the cage size and if their mass is also relatively small, it is possible that the guest-associated vibrational modes can exhibit significant anharmonic vibrations. To our knowledge, there has been very little attention paid to understanding this anharmonic behavior and its associated "off-center" guest behavior. In this work, we also report the results of a systematic, first-principles study of such anharmonic vibrational modes due to alkali atoms in a large Si 28 cage in the clathrates A 8 Si 136 (A = Na, K, Rb, Cs) and Na 24 Si 136 . One goal of this study is to explore the composition (x)-dependence of such a dynamically unstable, anharmonic guest-host potential. We find that the effective potential for this situation can be described with a large anharmonicity (large negative quartic coefficient) when the rattling impurity atom is changed from Cs (Rb, K) to Na. We also note that some previous reports of guest anharmonicity did not include a quartic anharmonicity interaction when examining the low-frequency rattling mode of the guests [14,15].

Computational Approach
Our calculated results are obtained by applying a computational quantum mechanical modeling method, which employs the local density approximation (LDA) to density functional theory (DFT). For such first-principles calculations, we focus on using the Vienna ab initio Simulation Package (VASP) [16][17][18][19][20][21], which exploits the Ceperley-Alder exchange-correlation potential, as well as pseudopotentials obtained using the projector augmented wave (PAW) method. The energy cutoff parameter when computing dispersion relations is set to 150 eV for silicon-based materials. This method has been extensively and successively tested in reported calculations. Specifically, K. Biswas et al. previously performed VASP determination of electronic structures regarding Na 16 Rb 8 Si 136 [12]. In their work, the calculated lattice parameter agrees well with the experimental result [13], while the calculated electronic density of states possessing a sharply peaked feature in the vicinity of the Fermi energy level can be qualitatively linked to the temperature-dependent Knight shift observed for the NMR-active nuclei in Na 16 Rb 8 Si 136 . Furthermore, K. Biswas et al. identified the low-frequency guest "rattling" modes from VASP-computed phonon dispersion relations when studying vibrational properties of Na 16 Rb 8 Si 136 [22], among which the estimated isotropic mean-square displacement amplitudes (U iso ) relating to Rb and Na are in good agreement with the experiment [11]. Hence, motivated by these reported works, we continued to follow the usage of the LDA to DFT method in the context of VASP while performing ab initio evaluations on electronic, vibrational, and anharmonic features of binary system A x Si 136 (A = Alkaline metal; 0 < x ≤ 24). For our study, a single crystallographic unit cell which contains 34 Si atoms is selected rather than a large clathrate unit cell structure involving 136 framework atoms.
The beginning step of our first-principles calculation starts with structural optimization. This optimization process is achieved by means of a conjugate gradient (CG) method, which relaxes the internal coordinates of the atoms confined in a fixed volume of the FCC unit cell. In the regime of type-II binary clathrate compounds described by cubic space group symmetry (Pm3d), the encapsulated guest atoms are allowed to move freely from their original point positioned by the cage center. It is worth mentioning that such a process for the relaxation and determination of the optimized structure must be repeated many times to achieve a global total minimum energy. Next, we fit limited pairs of LDA-calculated potential energy vs. volume (E, V) curves to the Birch-Murnaghan equation of state (EOS) to determine the equilibrium energy E 0 and the equilibrium volume V 0 . Accordingly, we initiate Brillouin zone (BZ) integration with the help of a 4 × 4 × 4 Monkhorst-Pack k-point grid [23] for the purpose of performing the relaxation and ultimately characterizing the equilibrium geometry. The total energy convergence criterion is adjusted to 10 −7 eV. As soon as the lattice parameter for a material of interest has been determined, the VASP code can readily perform further calculations, pertaining to the Fermi level position, the electronic band structure, and the associated electronic density of states (EDOS).
In the following, we apply our approach to analyze the vibrational features of A 8 Si 136 (A = Na, K, Rb, Cs), as well as Na 24 Si 136 . The typical technique for obtaining the rattling modes of the above binary clathrates consists of two steps. First, each atom in the single crystallographic unit cell is moved by a small finite displacement (U 0 = 0.02 Å), subject to determination by the dynamical matrix D(q). There are two reasons for designating the dynamical matrix requirement. One is based on the fact that D(q) is constructed at the gamma (Γ) point [q = (0,0,0)]. On the other hand, D(q) can still be calculated for nonzero q if it is assumed that the matrix elements of D(q) vanish for atoms separated by a distance greater than the third nearest neighbor distance [24]. In contrast to the previously mentioned BZ integration method, a new 2 × 2 × 2 k-point grid utilizing the k-points along certain high-symmetry directions is adopted to compute D(q). Second, it is clear that the diagonalization of D(q) allows us to find the vibrational eigenvalues (squared frequencies) and eigenvectors.

Electronic Properties
For the purpose of revealing how guest-framework interactions affect the electronic and vibrational properties of clathrate lattices, we initiate our beginning step to provide an EDOS study for the optimized geometries. A previous comprehensive description of the electronic structure associated with Na 16 Rb 8 Si 136 and K 16 Rb 8 Si 136 is presented in [12], where K. Biswas et al. calculate the electronic density of states, which displays two sharp peaks and a dip near the Fermi level. Such neighboring peaks have a separation (∆E) at an energy scale of several k B T, which is indicative of the temperature-dependent Knight shift effect observed experimentally for the NMR-active nuclei in Na 16 Rb 8 Si 136 [25]. Therefore, motivated by this model describing the two-peak form of the EDOS in the vicinity of the Fermi level, we perform a first-principles analysis of the EDOS of Na x Si 136 (x = 8, 16), which shows an analogous quasi-activation energy ∆E when taking the T-dependent Knight shift into consideration. Figure 1 illustrates the electronic density of states in the lower portion of the conduction band for Na 8 Si 136 and Na 16 Si 136 , where the corresponding Fermi energies are denoted by E F,2 and E F,4 . Here, quasi-activation energy ∆E 2,4 is predicted to lie in the range of approximately 0.4 eV to 0.46 eV. Moreover, the predicted ∆E 2 of Na 8 Si 136 and ∆E 4 of Na 16 Si 136 are similar to that of K 16 Rb 8 Si 136 (∆E~0.41 eV) from a previous work [12].
In view of the filled intermetallic clathrate Na 12 Si 136 , we found the existence of a direct "pseudogap," which is defined as the energy difference between the top of the valence band and the bottom of the conduction band at the same L point from the band structure (BS) perspective. Specifically, our DFT computation shows that the pseudogap is approximately equal to 0.48 eV. Figure 2 shows the calculated BS spanning over various Brillouin zone points, illustrating the direct pseudogap behavior. In view of the filled intermetallic clathrate Na12Si136, we found the existence of a direct "pseudogap," which is defined as the energy difference between the top of the valence band and the bottom of the conduction band at the same L point from the band structure (BS) perspective. Specifically, our DFT computation shows that the pseudogap is approximately equal to 0.48 eV. Figure 2 shows the calculated BS spanning over various Brillouin zone points, illustrating the direct pseudogap behavior. One of the electronic characteristics of the NaxSi136 system is understood through the variation of the Fermi energy level, whose position is shifted into a deep portion of the conduction band as x ranges from 8 to 12 and 16 in Figure 3. Furthermore, the shapes of the predicted EDOS profiles for these three filled clathrates are roughly identical and remain nearly independent of the guest  In view of the filled intermetallic clathrate Na12Si136, we found the existence of a direct "pseudogap," which is defined as the energy difference between the top of the valence band and the bottom of the conduction band at the same L point from the band structure (BS) perspective. Specifically, our DFT computation shows that the pseudogap is approximately equal to 0.48 eV. Figure 2 shows the calculated BS spanning over various Brillouin zone points, illustrating the direct pseudogap behavior. One of the electronic characteristics of the NaxSi136 system is understood through the variation of the Fermi energy level, whose position is shifted into a deep portion of the conduction band as x ranges from 8 to 12 and 16 in Figure 3. Furthermore, the shapes of the predicted EDOS profiles for these three filled clathrates are roughly identical and remain nearly independent of the guest One of the electronic characteristics of the Na x Si 136 system is understood through the variation of the Fermi energy level, whose position is shifted into a deep portion of the conduction band as x ranges from 8 to 12 and 16 in Figure 3. Furthermore, the shapes of the predicted EDOS profiles for these three filled clathrates are roughly identical and remain nearly independent of the guest composition x. In particular, the energy states occupying the valence band in three materials seem to be populated in a significantly consistent manner.
Materials 2019, 12, x FOR PEER REVIEW 5 of 12 composition x. In particular, the energy states occupying the valence band in three materials seem to be populated in a significantly consistent manner.

Vibrational Properties and Anharmonic Effects
Calculations of the "rattling" vibrational modes of the alkali metal atom guests in AxSi136 (A = Na, K, Rb, Cs; 0 < x < 24) can reveal very interesting basic physics, such as lattice dynamics and guesthost coupling. The main idea of the "rattler" concept originates from the fact that loosely bound guest atoms encapsulated in the "oversized" (28-atom) cages in the type-II clathrates vibrate and produce "localized" modes that are capable of efficiently scattering heat-carrying acoustic phonons [26,27]. Thus, the rattling behavior of the alkali metal atom guests can potentially participate in reducing the material thermal conductivity to a glass-like level, as suggested by Slack's "Phonon Glass Electron Crystal" (PGEC) criteria [28].
Our study of the vibrational properties of the filled clathrate AxSi136 is divided into two parts, as follows. First, the VASP code is used to calculate the phonon dispersion curves, along with the effective potential energy diagrams, which may contain apparent anharmonicity at zero temperature, depending on the guest (A) type. Second, the significant anharmonicity of the motion of the Na rattling guests in the Si28 cages in NaxSi136 is investigated in detail with the aid of a self-consistent phonon model to reveal its temperature dependence. The results predict that these anharmonic effects can strongly affect the force constants for the guest vibrational modes.
The LDA-calculated result correlating to the guest vibrational mode from the phonon dispersion relations in Na4Si136 agrees fairly well with the data obtained from inelastic neutron scattering (INS) experiments [29]. Specifically, the low-lying vibrational ("rattling") frequencies of the Na guests occur in a narrow band centered at approximately 50 cm −1 , which are approximately 4% lower than that of the experimentally observed mode (52 cm −1 ). In addition, an effective force constant K for the rattling frequency ω in the harmonic approximation (HA) may be obtained by assuming ω = (K/M) 1/2 , where M is the atomic mass of the guest. For Na4Si136, this gives K = 0.44 eV/Å 2 from the first-principles viewpoint for the Na vibrations in the Si28 cages.
To gain insight into the anharmonic effects associated with the Na guest vibrations in the Si28 cages, we carried out a method that can be summarized as follows. Using the LDA to generate the

Vibrational Properties and Anharmonic Effects
Calculations of the "rattling" vibrational modes of the alkali metal atom guests in A x Si 136 (A = Na, K, Rb, Cs; 0 < x < 24) can reveal very interesting basic physics, such as lattice dynamics and guest-host coupling. The main idea of the "rattler" concept originates from the fact that loosely bound guest atoms encapsulated in the "oversized" (28-atom) cages in the type-II clathrates vibrate and produce "localized" modes that are capable of efficiently scattering heat-carrying acoustic phonons [26,27]. Thus, the rattling behavior of the alkali metal atom guests can potentially participate in reducing the material thermal conductivity to a glass-like level, as suggested by Slack's "Phonon Glass Electron Crystal" (PGEC) criteria [28].
Our study of the vibrational properties of the filled clathrate A x Si 136 is divided into two parts, as follows. First, the VASP code is used to calculate the phonon dispersion curves, along with the effective potential energy diagrams, which may contain apparent anharmonicity at zero temperature, depending on the guest (A) type. Second, the significant anharmonicity of the motion of the Na rattling guests in the Si 28 cages in Na x Si 136 is investigated in detail with the aid of a self-consistent phonon model to reveal its temperature dependence. The results predict that these anharmonic effects can strongly affect the force constants for the guest vibrational modes.
The LDA-calculated result correlating to the guest vibrational mode from the phonon dispersion relations in Na 4 Si 136 agrees fairly well with the data obtained from inelastic neutron scattering (INS) experiments [29]. Specifically, the low-lying vibrational ("rattling") frequencies of the Na guests occur in a narrow band centered at approximately 50 cm −1 , which are approximately 4% lower than that of the experimentally observed mode (52 cm −1 ). In addition, an effective force constant K for the rattling frequency ω in the harmonic approximation (HA) may be obtained by assuming ω = (K/M) 1/2 , where M is the atomic mass of the guest. For Na 4 Si 136 , this gives K = 0.44 eV/Å 2 from the first-principles viewpoint for the Na vibrations in the Si 28 cages.
To gain insight into the anharmonic effects associated with the Na guest vibrations in the Si 28 cages, we carried out a method that can be summarized as follows. Using the LDA to generate the effective guest-host potential energy for an Na guest in an Si 28 cage acts as the first step. Then, to approximately treat the lowest order anharmonic effects present in the guest-host interaction, we expand this effective potential in a power series in the guest atom displacements, keeping only the second and fourth order terms. That is, we assume that the guest-host potential energy V eff can be approximately written as V eff = 1 2 K iso r 2 + 1 4 iso r 4 , where r is the displacement from equilibrium and K iso and iso represents empirical constants determined by fitting to the effective potential energy. Only even order terms are kept in this polynomial expansion because of the symmetry of the energy function.
In this paper, the procedure used to determine the guest-host framework effective potential energy curves displaces one Na guest atom inside an Si 28 cage while keeping all remaining atoms in the unit cell fixed. More concretely, each Na is displaced a distance from the center of the Si 28 cage. In all cases, this displacement distance is strictly restricted to the range of 0 to ∆r, where the upper limit ∆r describes the "excess" radius of that cage. Our previous work [15,26] showed that the "excess" radius of the Si 28 cage is ∆r = 1.83 Å for Na vibrations. Therefore, our calculation is performed in the Na x Si 136 system for Na guests moving in the X direction, but limited to the range of 0 to approximately 1.2 Å and repeated for guest compositions x = 4, 8, and 24. In the coordinate system considered here, the X and Y axes are chosen to be aligned in a plane that is parallel to the plane of the hexagon, while the Z axis is perpendicular to the hexagon in the 28-atom cage. The resulting effective potential energy curves for the compositions x = 4, 8, and 24 are shown in Figure 4. The results for motion along the X direction are shown. By symmetry, the effective potential energies for the Y-directed and Z-directed motions are nearly identical to that for X-directed motion, indicating the isotropic behavior of the guest vibration.
Materials 2019, 12, x FOR PEER REVIEW 6 of 12 effective guest-host potential energy for an Na guest in an Si28 cage acts as the first step. Then, to approximately treat the lowest order anharmonic effects present in the guest-host interaction, we expand this effective potential in a power series in the guest atom displacements, keeping only the second and fourth order terms. That is, we assume that the guest-host potential energy Veff can be approximately written as Veff = ½ Kisor 2 + ¼ Ɛisor 4 , where r is the displacement from equilibrium and Kiso and Ɛiso represents empirical constants determined by fitting to the effective potential energy.
Only even order terms are kept in this polynomial expansion because of the symmetry of the energy function.
In this paper, the procedure used to determine the guest-host framework effective potential energy curves displaces one Na guest atom inside an Si28 cage while keeping all remaining atoms in the unit cell fixed. More concretely, each Na is displaced a distance from the center of the Si28 cage. In all cases, this displacement distance is strictly restricted to the range of 0 to Δr, where the upper limit Δr describes the "excess" radius of that cage. Our previous work [15,26] showed that the "excess" radius of the Si28 cage is Δr = 1.83 Å for Na vibrations. Therefore, our calculation is performed in the NaxSi136 system for Na guests moving in the X direction, but limited to the range of 0 to approximately 1.2 Å and repeated for guest compositions x = 4, 8, and 24. In the coordinate system considered here, the X and Y axes are chosen to be aligned in a plane that is parallel to the plane of the hexagon, while the Z axis is perpendicular to the hexagon in the 28-atom cage. The resulting effective potential energy curves for the compositions x = 4, 8, and 24 are shown in Figure 4. The results for motion along the X direction are shown. By symmetry, the effective potential energies for the Y-directed and Z-directed motions are nearly identical to that for X-directed motion, indicating the isotropic behavior of the guest vibration. Here, for Na4Si136 and Na8Si136, the results of this procedure predict that the effective potential energy has an obvious "bump" at the cage center (zero Na displacement in Si28). Additionally, each of these energy curves in NaxSi136 (x = 4, 8) has a "Mexican-hat" (double valley) shape. This form of the potential energy shape shows that Na guests in the Si28 cages are dynamically unstable at the cage center. Thus, these calculations predict that, at elevated temperatures, thermal excitation can induce Na guests to move off the cage center to a site with a lower potential energy. In other words, these Here, for Na 4 Si 136 and Na 8 Si 136 , the results of this procedure predict that the effective potential energy has an obvious "bump" at the cage center (zero Na displacement in Si 28 ). Additionally, each of these energy curves in Na x Si 136 (x = 4, 8) has a "Mexican-hat" (double valley) shape. This form of the potential energy shape shows that Na guests in the Si 28 cages are dynamically unstable at the cage center. Thus, these calculations predict that, at elevated temperatures, thermal excitation can induce Na guests to move off the cage center to a site with a lower potential energy. In other words, these calculations find that the "localized" off-center minimum valley of the effective potential curve is approximately 0.017 eV and 0.01 eV below zero for Na 4 Si 136 and Na 8 Si 136 , respectively. That is, the effective potential energy valley depth is predicted to be approximately 0.66 k B T and 0.39 k B T at room temperature for Na 4 Si 136 and Na 8 Si 136 , respectively. It is likely that the minimum energy valley causes the "off-center" position of Na guest to form due to increased temperature. This kind of anharmonic effect might be correlated with the measured strong temperature dependence of the isotropic atomic displacement parameter U iso , as reported in References [5,30].
This predicted "off-center" position for Na in the Si 28 cages in Na x Si 136 can, in principle, be correlated with the experimentally measured guest equilibrium location at various temperatures. However, the disorder of the guest position makes the understanding of such "off-center" sites complex. Similar effects have been observed in some type-I clathrates, including Sr 8 Ga 16 Ge 30 , Ba 8 Ga 16 Ge 30 , and Eu 8 Ga 16 Ge 30 [5,11,13,15,25,27,31]. The arrow in the inset of Figure 4 indicates that the "localized" off-center minimum valley is shifted towards the cage center when the Na content is changed from 4 to 8 to 24, despite the fact that apparently weakened double-valley behavior exists in the effective potential energy curve of Na 24 Si 136 compared to others.
Considerable experimental work on Na x Si 136 by Beekman and collaborators focused on detecting the isotropic, "off-center" displacement of Na in Si 136 and studying the system in significant detail [4,[32][33][34][35]]. An early report by them on the Na vibrations in Si 28 showed a nearly x-independent "off-center" displacement measured to be approximately 0.4-0.5 Å [32]. Our first-principles calculations correlate well with these experimental results. This is shown in Figure 5, where both the XRD data and our first-principles predictions of the Na displacements in Si 28 are shown as a function of guest stoichiometry amount x. For cases x = 12 and x = 16, these calculations predict an isotropic displacement that ranges from 0.50 Å for Na 12 Si 136 to 0.42 Å for Na 16 Si 136 . Another report [34] states elsewhere that an attractive interaction occurs between the neighboring Na guests (situated in adjacent polyhedron cage Si 28 ) in the presence of the above discussed anharmonic guest-host potential energy, which causes the "off-center" displacement of the guests to appear. calculations find that the "localized" off-center minimum valley of the effective potential curve is approximately 0.017 eV and 0.01 eV below zero for Na4Si136 and Na8Si136, respectively. That is, the effective potential energy valley depth is predicted to be approximately 0.66 kBT and 0.39 kBT at room temperature for Na4Si136 and Na8Si136, respectively. It is likely that the minimum energy valley causes the "off-center" position of Na guest to form due to increased temperature. This kind of anharmonic effect might be correlated with the measured strong temperature dependence of the isotropic atomic displacement parameter Uiso, as reported in References [5,30]. This predicted "off-center" position for Na in the Si28 cages in NaxSi136 can, in principle, be correlated with the experimentally measured guest equilibrium location at various temperatures. However, the disorder of the guest position makes the understanding of such "off-center" sites complex. Similar effects have been observed in some type-I clathrates, including Sr8Ga16Ge30, Ba8Ga16Ge30, and Eu8Ga16Ge30 [5,11,13,15,25,27,31]. The arrow in the inset of Figure 4 indicates that the "localized" off-center minimum valley is shifted towards the cage center when the Na content is changed from 4 to 8 to 24, despite the fact that apparently weakened double-valley behavior exists in the effective potential energy curve of Na24Si136 compared to others.
Considerable experimental work on NaxSi136 by Beekman and collaborators focused on detecting the isotropic, "off-center" displacement of Na in Si136 and studying the system in significant detail [4,[32][33][34][35]]. An early report by them on the Na vibrations in Si28 showed a nearly x-independent "offcenter" displacement measured to be approximately 0.4-0.5 Å [32]. Our first-principles calculations correlate well with these experimental results. This is shown in Figure 5, where both the XRD data and our first-principles predictions of the Na displacements in Si28 are shown as a function of guest stoichiometry amount x. For cases x = 12 and x = 16, these calculations predict an isotropic displacement that ranges from 0.50 Å for Na12Si136 to 0.42 Å for Na16Si136. Another report [34] states elsewhere that an attractive interaction occurs between the neighboring Na guests (situated in adjacent polyhedron cage Si28) in the presence of the above discussed anharmonic guest-host potential energy, which causes the "off-center" displacement of the guests to appear.  Utilizing a similar method to study the anharmonic effects associated with the rattling guest atoms K, Rb, and Cs in Si 136 , we show the isotropic guest-framework effective potentials for K 8 Si 136 , Rb 8 Si 136 , and Cs 8 Si 136 (see Figure 6). The effective potential energy is again symmetric, so only the spring harmonic constant K iso and the quartic anharmonic coefficient iso need to be considered. That is, the LDA effective potential energy can be fitted to the equation for V eff , as discussed above. We find that increasing the atomic mass of the rattling guest results in a higher potential energy for a given displacement. After fitting the LDA effective potential to V eff , the second order and fourth order coefficients for the guest-framework interaction potential (Figures 4 and 6) are obtained; the results are summarized in Table 1. harmonic constant Kiso and the quartic anharmonic coefficient Ɛiso need to be considered. That is, the LDA effective potential energy can be fitted to the equation for Veff, as discussed above. We find that increasing the atomic mass of the rattling guest results in a higher potential energy for a given displacement. After fitting the LDA effective potential to Veff, the second order and fourth order coefficients for the guest-framework interaction potential (Figures 4 and 6) are obtained; the results are summarized in Table 1. Table 1. Fitting parameters of the guest-host interaction energy potentials for A8Si136 (A = Na, K, Rb, Cs) and Na24Si136.
To examine the anharmonic effects in the guest-host interaction, the self-consistent phonon (SCP) model [36] is used to acquire the anharmonic nature of the phonon rattling modes at nonzero temperature. Accordingly, the temperature-dependent anharmonic rattling frequency (Ω(T)) at 300 K appearing in Table 1 (denoted by Ω(T = 300 K)) can be quantitatively understood by correlating the thermal displacement Uiso to the spring constant and quartic anharmonic coefficient arising from potential energy Veff: MΩ 2 (T) = Kiso + Ɛiso<u 2 (k,j)>. Here, the average mean square amplitude of the phonon atomic displacement <u 2 (k,j)> (= Uiso) relating to the wavevector (k) and phonon branch (j) is defined via the Debye model [13], and M denotes the atomic mass of the guest rattler. Figure 6. The potential energy curves for alkali atoms K, Rb, and Cs inside the Si framework. X represents the motion direction with respect to the alkali atom. Table 1. Fitting parameters of the guest-host interaction energy potentials for A 8 Si 136 (A = Na, K, Rb, Cs) and Na 24 Si 136 .

Material
K iso (eV/Å 2 ) ε iso (eV/Å 4 ) ω ph (cm −1 ) Ω (T = 300 K) (cm −1 ) (K iso /M) 1/2 (cm −1 ) To examine the anharmonic effects in the guest-host interaction, the self-consistent phonon (SCP) model [36] is used to acquire the anharmonic nature of the phonon rattling modes at nonzero temperature. Accordingly, the temperature-dependent anharmonic rattling frequency (Ω(T)) at 300 K appearing in Table 1 (denoted by Ω (T = 300 K) ) can be quantitatively understood by correlating the thermal displacement U iso to the spring constant and quartic anharmonic coefficient arising from potential energy V eff : MΩ 2 (T) = K iso + iso <u 2 (k,j)>. Here, the average mean square amplitude of the phonon atomic displacement <u 2 (k,j)> (= U iso ) relating to the wavevector (k) and phonon branch (j) is defined via the Debye model [13], and M denotes the atomic mass of the guest rattler.
In contrast to the vibration of the Na rattling mode in the Si 28 polyhedron cage, a very slight variation in temperature-dependent Ω(T) is observed for the localized Rb impurity in Figure 7. The guest "off-center" displacement has a larger impact on the derived anharmonic rattling frequency and correlated thermal average of <u 2 (k,j)>. Our numerical results describing Ω(T) in Figure 7 give approximately 42.55 cm −1 (≈5.27 meV) for the Rb vibration frequency at T = 0, which is approximately 4.2% smaller than that obtained at T = 300 K. However, the predicted values of Ω(T) at 300 K are significantly higher than that at absolute zero for sodium filled clathrate. This strong temperature dependence of Ω(T) in Na 24 Si 136 indicates that static disorder effects might impact the isotropic atomic displacement parameters (ADPs). In other words, the average mean square isotropic atomic displacement for Na that resides in the "off-center" system has a boosted value, compared to <u 2 (k,j)> for the "on-center" vibration at the same temperature. Our work predicts that this difference in thermal average <u 2 (k,j)> is nearly close to the square of the "off-center" displacement investigated above (see Figure 5). Intensive discussion on the discrepancy happening to U iso is given below. In contrast to the vibration of the Na rattling mode in the Si28 polyhedron cage, a very slight variation in temperature-dependent Ω(T) is observed for the localized Rb impurity in Figure 7. The guest "off-center" displacement has a larger impact on the derived anharmonic rattling frequency and correlated thermal average of <u 2 (k,j)>. Our numerical results describing Ω(T) in Figure 7 give approximately 42.55 cm −1 (≈ 5.27 meV) for the Rb vibration frequency at T = 0, which is approximately 4.2% smaller than that obtained at T = 300 K. However, the predicted values of Ω(T) at 300 K are significantly higher than that at absolute zero for sodium filled clathrate. This strong temperature dependence of Ω(T) in Na24Si136 indicates that static disorder effects might impact the isotropic atomic displacement parameters (ADPs). In other words, the average mean square isotropic atomic displacement for Na that resides in the "off-center" system has a boosted value, compared to <u 2 (k,j)> for the "on-center" vibration at the same temperature. Our work predicts that this difference in thermal average <u 2 (k,j)> is nearly close to the square of the "off-center" displacement investigated above (see Figure 5). Intensive discussion on the discrepancy happening to Uiso is given below. A previous XRD study by M. Beekman et al. [5] stated that the strong T-dependence of <u 2 (k,j)> for Na@Si28 in NaxSi136 might be attributable to thermal excitation, which causes the guest to explore the minimum of the "Mexican-hat" shape potential (see Figure 4). That is, the anharmonic and dynamically unstable guest-host interaction makes it possible for the Na rattler to seek a larger vibration volume at elevated temperatures. In other words, the experimentally measured atomic displacement parameters are closely related to the temperature-dependent Ω(T) at room temperature instead of the DFT-calculated rattling mode ωph when considering the harmonic oscillator model [14], even if there is a large discrepancy between the Ω(T = 300 K) and ωph values.
For the specific case of Na24Si136, our VASP-calculated result predicts a harmonic rattling frequency that is approximately 15% larger than the anharmonic counterpart Ω(T) obtained at 300 K. This difference in frequency can lead to differences with respect to the corresponding atomic displacement parameters Uiso. To account for this discrepancy, this effect is interpreted in the picture of the thermal displacement parameter while taking "off-center" displacement due to static disorder into account. Mathematically, we postulate that the relation Uiso(Ω(T = 300 K)) ≈ Uiso(ωph) + d 2 can be used to account for the existence of the strong ADPs observed by means of XRD, INS, or temperaturedependent heat capacity (Cp) measurements [5,31,37]. Numerical values of Ω(T = 300 K) computed by the A previous XRD study by M. Beekman et al. [5] stated that the strong T-dependence of <u 2 (k,j)> for Na@Si 28 in Na x Si 136 might be attributable to thermal excitation, which causes the guest to explore the minimum of the "Mexican-hat" shape potential (see Figure 4). That is, the anharmonic and dynamically unstable guest-host interaction makes it possible for the Na rattler to seek a larger vibration volume at elevated temperatures. In other words, the experimentally measured atomic displacement parameters are closely related to the temperature-dependent Ω(T) at room temperature instead of the DFT-calculated rattling mode ω ph when considering the harmonic oscillator model [14], even if there is a large discrepancy between the Ω (T = 300 K) and ω ph values.
For the specific case of Na 24 Si 136 , our VASP-calculated result predicts a harmonic rattling frequency that is approximately 15% larger than the anharmonic counterpart Ω(T) obtained at 300 K. This difference in frequency can lead to differences with respect to the corresponding atomic displacement parameters U iso . To account for this discrepancy, this effect is interpreted in the picture of the thermal displacement parameter while taking "off-center" displacement due to static disorder into account. Mathematically, we postulate that the relation U iso (Ω (T = 300 K) ) ≈ U iso (ω ph ) + d 2 can be used to account for the existence of the strong ADPs observed by means of XRD, INS, or temperature-dependent heat capacity (C p ) measurements [5,31,37]. Numerical values of Ω (T = 300 K) computed by the self-consistent phonon model and ω ph given by the VASP code can be plugged into the classical equation U iso (ω) = k B T/Mω 2 within a high-temperature approximation [38]. For the case of Na 4 Si 136 , the calculation leads to the corresponding ADPs, which give approximately 0.497 Å 2 and 0.122 Å 2 for U iso (Ω (T = 300K) ) and U iso (ω ph ) in Na 4 Si 136 , respectively. Therefore, the calculated Na "off-center" displacement d in Na 4 Si 136 is approximately equal to 0.61 Å, which is reasonably comparable to the experimental counterpart (~0.52 Å) found in Figure 4. Similarly, our calculated d value by virtue of the suggested model for Na 24 Si 136 is 0.28 Å, which is in good agreement with the results estimated by Rietveld analysis in Ref. [32]. Moreover, the DFT-calculated "off-center" displacement of Na in Na 24 Si 136 (see Figure 5) is roughly equal to 0.36 Å, which is slightly higher than our numerically determined d parameter.

Conclusions
We have employed ab initio DFT to systematically investigate the electronic, vibrational, and anharmonic properties of the binary clathrate system A x Si 136 . One of the main electronic characteristics is that the electronic density of states remains nearly unchanged as the Na composition x changes from 8 to 12 to 16 in Na x Si 136 . Furthermore, the determined Fermi energy level E F is shifted to a deep portion of the conduction band in the same range (8 ≤ x ≤ 16), while a similar effect characterized by a "double-peak" along with a dip near the Fermi energy in the EDOS calculation is observed in Na 8 Si 136 and Na 16 S i136 , in analogy to the (Na,K) 16 Rb 8 Si 136 case reported previously. Using an effective potential-energy curve to explore the Na vibration in the Si 28 cage, we point out the existence of a dynamically unstable state in the "double-well" structure. This is also indicative of an "off-center" displacement for guests subjected to elevated temperatures. In addition, the guest vibrations are nearly isotropic in the same stoichiometry Na x Si 136 (0 < x ≤ 24), independent of the composition x. Furthermore, the DFT-determined "off-center" displacement for the Na vibration in the Si 28 cages is nearly independent of guest content, in fair agreement with the XRD data. It is noted that the quartic anharmonicity effect induced by the guest-framework interaction becomes comparably weakened when the ratio of the lowest even anharmonic coefficient to the spring constant decreases in the A 8 Si 136 system. Making use of the SCP model, the determined temperature-dependent rattling frequency Ω(T) of Na in Na 24 Si 136 implies strong static disorder effects associated with the average mean square isotropic atomic displacement U iso . This is manifested by the fact that the values of Ω(T) for Na at 300 K are approximately six times larger than the value at absolute zero. To account for the experimental overestimation of U iso at room temperature, we suggest that the existence of strong ADPs observed by means of XRD, INS, or temperature-dependent heat capacity (C p ) measurements in apparently anharmonic clathrate Na 24 Si 136 is offset by the square of the Na "off-center" displacement in the Si 28 cage. Considering the ability of tuning the alkaline guest (A) species and its composition x in such a "cage-structured" configuration of A x Si 136 , theoretical investigation on the corresponding rattling behavior, as well as guest "off-center" displacement, may contribute to boosting the search for novel thermoelectric materials. A good performance of TE material is manifested and parameterized by high electric conductivity and low thermal conductivity, according to "Phonon Glass Electron Crystal" criteria [28]. Therefore, fabrication of this class of clathrates may be used in TE technology involving the conversion of heat to electricity, and in optoelectronics technology, including the conversion of light to electricity. Accordingly, recent work on alloyed Si-Ge type II clathrate [39] presents many promising optoelectronic applications, which are pertinent to single junction PV cells with optimized band gaps and Si-Ge multijunction cells.
Author Contributions: D.X., the first author, undertook all of the (sometimes computationally intense) calculations that are reported in the paper. He also wrote the majority of the first draft and participated in the editing of that draft which produced the final, submitted version. C.W.M., the second author, is a Professor of Physics and Astronomy at Texas Tech University. He is the PhD research advisor to D.X.
Funding: This research received no external funding.