Lattice Dynamics Study of Phonon Instability and Thermal Properties of Type-I Clathrate K₈Si46 under High Pressure.

For a further understanding of the phase transitions mechanism in type-I silicon clathrates K₈Si46, ab initio self-consistent electronic calculations combined with linear-response method have been performed to investigate the vibrational properties of alkali metal K atoms encapsulated type-I silicon-clathrate under pressure within the framework of density functional perturbation theory. Our lattice dynamics simulation results showed that the pressure induced phase transition of K₈Si46 was believed to be driven by the phonon instability of the calthrate lattice. Analysis of the evolution of the partial phonon density of state with pressure, a legible dynamic picture for both guest K atoms and host lattice, was given. In addition, based on phonon calculations and combined with quasi-harmonic approximation, the specific heat of K₈Si46 was derived, which agreed very well with experimental results. Also, other important thermal properties including the thermal expansion coefficients and Grüneisen parameters of K₈Si46 under different temperature and pressure were also predicted.

1. Introduction

Nanostructured type I clathrates are composed of two 20-atoms (small) and six 24-atoms (large) cages of group-IV elements which can host different kinds of atoms. The framework atoms admit partial substitution by other atomic species. Such good tailorability enables these clathrate compounds to have extensive potential applications in areas of superconductivity [1], wide-band-gap semiconductors [2], optoelectronics [3], magnetism [4], thermoelectric [5] and photovoltaics [6], etc. As illustrated in Figure 1, there are three distinct Wyckoff symmetry sites, i.e., 6c, 16i and 24k for the framework atoms. Encaptured metal atoms located at the center of small and large cages correspond to 2a sites and 6d sites, respectively. Since type-I Ba8Si46 had been realized in a multi-anvil press by the group of Yamanaka in 2000 [7], high pressure technique has opened up a new door to synthesize these kind of clathrates compounds and provided new ideas to explore the interaction mechanism between guest atoms and the host lattice. During the high pressure experiments in studying the stability and compressibility of these silicon clathrates, some clathrates doped with large guest atoms such as K8Si46 [8,9], Ba8Si46 [10,11,12,13,14,15], I8Si44I2 [16,17], Rb6Si46 [18] exhibit an abrupt pressure induced volume collapse transition while the overall crystal symmetry is preserved. Various mechanisms (e.g., the change of hybridization between guest atoms and frame cages [12] and an electronic topological transition [14,15]) for this high pressure phase transition have been proposed to explain this phenomenon for Ba8Si46. In our recent work [19], detailed investigations have been performed to study the effect of forming lattice vacancies on the mechanical and electronic properties of Ba8Si46 under high pressure. The results indicate that the compressbility of Ba8Si46 is governed by the Si framework and the pressure induced escape of host Si atoms would cause the sudden volume collapse during the compress, which is consistent with the explanation from Iitaka et al. according to their atomistic model [20]. A more coherent picture of the phase transition is given by the newest experimental study of the clathrate collapse in mixed Ba8(Si,Ge)46 clathrates which suggests that the volume collapse is a second order transition through a Landau modeling ,with a symmetry-breaking mechanism related to the Ba displacement initiated either by vacancy creation or by local distortion of the framework structure [21].

Figure 1

The sketch map of ball and stick representation of the type I clathrate K8Si46 with an illustration of Wyckoff sites of the structure.

In the case of K8Si46, a noticeable change in the compressibility had been found at around 15 GPa followed by an abrupt change of volume at around 20 GPa. Above 32 GPa, the sample transformed into an amorphous phase. Earlier ab initio phonon band structure calculations using finite displacement method [22] had found that vibrational frequencies of the K atoms in the large cavities became imaginary at a pressure of 16 GPa which suggested a phase transition triggered by positional disordered guest atoms [8]. However, it is known that the supercell should be constructed in a phonon calculation within this method [22], which is obviously unbearable for an originally large primitive cell of K8Si46. Therefore, in their calculations, just primitive cell was supposed to be used. Although the primitive K8Si46 cell under ground state has a large lattice constant of somewhat more than 10 Å, the magnitude of the force constants beyond that distance is generally expected to be negligible. In a covalent semiconductor and complex structured system, the required cutoff radius should be larger than 10 Å as a rule, thus such a treatment may cause unreliable results when the K8Si46 cell is compressed under high pressure (the lattice constant of K8Si46 under 15 GPa is about 9.68 Å). For a confirmation of this point, we have also performed the phonon calculation within finite displacement method using a unit cell of K8Si46 and indeed found that a large number of imaginary frequencies appeared when the pressure was applied up to 15 GPa. So, in the present work, we employ the density functional perturbation theory (DFPT) [23,24] which avoids the use of a supercell and allows calculation of the dynamical matrix exactly at an arbitrary vector to study the stability of K8Si46 under pressure. For a perfect crystal, the dynamical instability of the phonon is found to be originated from the framework of K8Si46 which is quite different to the results reported by J. S. Tse et al. [8]. From the obtained partial phonon density of the state, a clearer dynamic picture of guest K atoms and host lattice under high pressure is given. Besides, based on the phonon calculation, we derived the thermal properties of K8Si46 from quasi-harmonic approximation which were believed to offer a useful reference to design and synthetize a new ternary type-I Si clathrate based on K8Si46 with enhanced thermoelectric performance.

2. Computational Methods

Phonon dispersion spectrum was calculated by using the linear-response method [25,26] within the density functional perturbation theory (DFPT) [23,24]. The full phonon dispersion curves were obtained through Fourier interpolation. Norm-conserving pseudopotentials generated using the kinetic energy optimization scheme developed by Lin et al. [27] and Lee [28] were employed to describe the ion-electron interactions with an energy cutoff of 800 eV to expand the valence electronic wave functions. Monkhorst-Pack k-points mesh of 6 × 6 × 6 had been chosen. The electronic exchange-correlation interactions were treated within the local density approximation (LDA) [29]. During the pseudopotentials calculations, pseudo atomic calculations were performed for K (3s, 3p, 4s) and Si (3s, 3p). By employing the Parrinello–Rahman method [30,31], hydrostatic pressure load on crystal was realized within the variable cell approach. Both the cell parameters and the atomic internal coordinates were fully relaxed at each target external pressure by applying Broyden, Fletcher, Goldfarb and Shanno (BFGS) scheme [32]. All these total energy electronic structure calculations were carried out by Cambridge Serial Total Energy Package (CASTEP) code [33,34].

3. Results and Discussion

3.1. Phonon Spectra

Within the framework of DFPT method, phonon frequencies are computed as second-order derivatives of the total energy with respect to a given perturbation in the form of atomic displacements. The force constants matrix can be obtained by differentiating the Hellmann-Feynman forces on atoms with respect to ionic coordinates. Based on the interatomic force constants, we can obtain the phonon spectra by using Fourier interpolation with specific treatment of the long-range dipole-dipole interaction [35]. The phonon band structure calculations were performed up to 40 GPa with an interval of 10 GPa. A more careful calculation at 39 GPa was performed for determining the exact pressure value when the phonon became instable.

The obtained ground state phonon dispersion relations in K8Si46 crystal are shown in Figure 2a. Clear flat vibrational bands can be easily recognized in the low frequency region around 100 cm−1. By analysis of partial phonon density of state (PPDOS) as shown in Figure 3a, we find that these low frequencies are corresponding to a sharp peak located at about 98.8 cm−1 which are originated from the localized vibration of 6d sites K atoms at Si24 cages. Another somewhat weaker but still clear peak at 128.2 cm−1 also from K(6d) can be attributed to the asymmetry of the large Si24 cage which takes ellipsoidal shape. A similar phenomenon had also been observed in another type-I silicon-clathrate Na8Si46 reported by Li et al. [36]. Moreover, a strong and broad peak centered at 172.1 cm−1 is found to be contributed by the mixing vibrations of K atoms at Si20 cages and the framework Si atoms indicates the intense interaction hybridization between them, this is because the size of Si20 cage is much smaller than Si24 cage which yields a shorter interatomic distance of K-Si. Experimentally, the measured Raman spectrum of K8Si46 showed noticeable peaks at 94, 119 and 177 cm−1 related to the vibration of K atoms [9], which is consistent with our results. By inelastic neutron scattering, Mélinon et al. [37] and Reny et al. [38] obtained very similar vibrational spectrum of K8Si46 and found two identified peaks centered at about (100 cm−1, 170 cm−1) and (98 ± 2 cm−1, 161 ± 5 cm−1) respectively. Our calculated results perfectly reproduced the main characteristics of experimental observation. Under a pressure of 30 GPa, due to the shrinkage of the silicon cages under compress, the PPDOS shows strong mixing of K and Si vibrations. The frequency of K atoms’ motion at both 2a and 6d sites became higher, which indicates a further localization of these atoms, as shown in Figure 3b. However, it is noted that the appearance of massive low frequencies that originate from the framework atoms suggests a collective “soften” of Si-Si bonds which will finally make the host lattice unstable. As the pressure is applied at a value of 40 GPa, it can be seen clearly from Figure 2c that the frequencies around the M symmetry point become imaginary, which clearly shows the instability of clathrate framework. As illustrated in in Figure 3c, the PPDOS under 40 GPa also shows a dramatic reduction of frequencies from host lattice as presented in the phonon spectrum. Our calculation results show that a mechanical instability of the silicon framework is believed to be responsible for the pressure-induced volume collapse at about 20 GPa and subsequent amorphization at 32 GPa of K8Si46 observed experimentally.

Figure 2

The phonon dispersion relation of K8Si46 at (a) 0 GPa; (b) 30 GPa and (c) 40 GPa.

Figure 3

The partial phonon density of state of K8Si46 at (a) 0 GPa; (b) 30 GPa and (c) 40 GPa.

However, it is noted that the unstable pressure given in the present work is obviously larger than the experimental observation. This is because our calculations are performed using a perfect K8Si46 crystal from the view of lattice dynamics without consideration of other possible transition mechanisms (e.g., vacancies formation [20], local symmetry-breaking [21] or an electronic topological transition [14,15], etc.) associated with the pressure collapse of the clathrate structure. If multi-mechanisms are involved in phase transition, the value of transition pressure would be affected a lot. For instance, the lattice vacancies are actually very likely to be produced in these type-I clathrates, especially in 6c sites. The experimental observed Cs8Sn44 was formed from the missing two Sn atoms in the 6c sites [39]. Besides, theoretical calculations by Iitaka et al. also showed that 6c sites lattice vacancies formation under high pressure was indeed energetically preferable. By the model of partially occupied Si sites they explained the transition pressure and change of Raman spectra of both K8Si46 and Ba8Si46 [20]. Moreover, our recent work showed that 6c sites lattice vacancies increased the compressibility of clathrate greatly while guest atoms vacancies hardly had any influence on this property [19]. Experimentally, by performing high quality in situ high-pressure angle-dispersive X-ray powder diffraction measurements, Li et al. found a highly disordered Si framework from analysis of the obtained anomalously large Si thermal parameters [13]. Also, present lattice dynamics calculation for K8Si46 shows that the Si-Si bond “softens” under high pressure which again provides theoretical possibility for the formation of lattice vacancies. If one considers this mechanism, the volume collapse pressure of K8Si46 is believed to be reduced. Thus, in view of these results, guest K atoms displacement induced phonon instability from earlier ab initio phonon band structure calculations [8] is indeed more likely to be caused by the disadvantages of the finite displacement method in treating these large cell calthrate compound, especially under high pressure. Our calculated results based on DFTP method obviously gives a more convincing and clear physical picture for the instability of K8Si46 under high pressure.

3.2. Thermal Properties from Quasi-Harmonic Approximation

The results of calculated phonon spectra and phonon density of state can be used to compute the thermodynamic properties using the quasi-harmonic approximation (QHA) [40]. In the QHA, the phonon Helmholtz free energy is given by: where kB is the Boltzmann constant. ħ is Planck’s constant and f(ω) is the phonon density of states (PDOS). Through a series calculation of PDOS of K8Si46 with different volumes, the volume dependence of Helmholtz free energy Fvib(V,T) can be obtained. Then the vibrational contribution to the entropy, the specific heat at constant volume and isothermal bulk modulus can be derived by:

The Grüneisen parameters can be computed by the volume derivative of (−TS):

Then, the volume coefficient of thermal expansion and constant pressure heat capacity (Cp) follows:

From QHA calculation, zero-point energy Fvib (T = 0) of K8Si46 compound is determined to be 2.688 eV. Moreover, the calculated variation of volume thermal expansion coefficient α with temperature under different pressures are illustrated in Figure 4 from which we can find that α increases rapidly with temperature below about 200 K and pressure imposes a strong restraint on the lattice expansion. At room temperature and zero pressure, the α is predicted as 6.26 × 10−5 K−1, corresponding to a linear thermal expansion coefficient α as 2.09 × 10−5 K−1 which is very close to the experimental value of Na8Si46 (about 2.0 × 10−5 K−1) given by Qiu et al. [41]. In addition, another important thermodynamic quantity of Grüneisen parameters which is difficult to determine experimentally can also be predicted by QHA method. The obtained temperature dependencies of Grüneisen parameters under pressure of 0, 10, 20 and 30 GPa are presented in Figure 5. It can be found that the Grüneisen parameters become almost constant in relation to the varied temperature under high pressure. At ambient conditions, the Grüneisen parameter is found to be 2.47. For congener compound Na8Si46 under same condition, this value was reported as 2.68 [41].

Figure 4

The calculated variation of volume coefficient of thermal expansion of K8Si46 with temperature under 0, 10, 20 and 30 GPa.

Figure 5

The calculated variation of Grüneisen parameters of K8Si46 with temperature under 0, 10, 20 and 30 GPa.

These similarities in thermdynamic properties of type I clathrates doped with same group elements have also been reported by many experimental works, for example, the thermal expansion coefficients of Ba(Sr)8Ga16Ge30 and Sr8Ga16 Ge30 [42] were just found to be almost identical to each other, and so were Rb8Sn44□2 and Cs8Sn44□2 (□ means lattice vacancy) [43]. Consequently, in general, the thermal expansion coefficient of intermetallic clathrates is assumed to mainly depend on the bonding of the framework atoms. The nature of the guest atoms just make a small contribution due to their weaker ionic bonding to the host structure. However, in the case of Ba8Si46, which has been attracting extensive attention due to the discovery of its superconductivity, the measured α at room temperature was found to be only 1.2 × 10−5 K−1. Besides, a drop of α of Ba8Si46 occurred at the superconducting transition temperature and the frequency-dependent Grüneisen parameter of Ba8Si46 indicated strong anharmonicity of the lattice vibrations for low energy mode with a value of up to 8.6 while higher energy modes are much less anharmonic with a value somewhat below 2 [44]. These novel results of Ba8Si46 are quite different from those of Na(K)8Si46 which reveals the non-negligible role of different hybridization between the guest and host atoms in studying the vibrational properties for doped clathrate compound. Moreover, the vacancies which are very likely to be formed in clathrate compounds can also significantly influence the system vibrational properties. Their presence can decrease the average bonding strength among host atoms and induce the displacement of guest atoms from their ideal positions. Thus, the anharmonicity of system is enhanced which would yeild an increment of the thermal expansion coefficient and Grüneisen parameters at zero pressure that had been identified by the experimental work for type-I Ge-based clathrates [45]. However, this effect became more unintelligible under high pressure because the formation of vacancies can also increase the compressbility of the clathrate which would lead the volume collapse under compress. So, further theoretical simulation with considering the formation of vacancies in hosts framework is expected to explore the effects of these vacancies on the vibrational properties of clathrate system under high pressure.

Because the experimental specific heat was conducted at constant pressure namely Cp, in Figure 6, the resulting dependence of Cp on temperature calculated from QHA method at zero pressure is illustrated, the experimental data up to room temperature from Stefanoski have also been plotted together for a comparison [46]. It can be found that they agree with each other excellently (for instance, the specific heat at room temperature determined by present work is 1151.5 J·mol−1·K−1, experimental measured value is 1159.6 J·mol−1·K−1) which indicates the validity of present lattice dynamic simulation within DFTP. At low temperatures, the temperature dependence of specific heat presents a similar behavior to thermal expansion coefficient. When the temperature is higher than about 400 K, specific heat gradually approaches the Dulong-Petit limit, i.e., 3nNAkB (about 1397.2 J·mol−1·K−1) which is followed by all solid at high temperatures. In our former study [47], we found the specific heat of Na8Si46 predicted by quasi-harmonic Debye model was underestimated obviously at low temperature which revealed the limitation of the Debye model in dealing with these doped clathrate compounds. The experimental heat capacities of Na8Si46 were finally reproduced by treating the special “rattle” modes of captured Na atoms in cages as Einstein oscillators. In the present work, a full phonon calculation within DFPT method can give a rather exact description of vibrational properties of these clathrate compounds which avoids a discriminatory treatment of the host lattice and encapsulated atoms. The specific heat of K8Si46 under pressure of 30 GPa is also presented in Figure 6, from which it can be found that pressure can decrease the specific heat considerably due to the suppression of lattice vibration.

Figure 6

The value of CP of K8Si46 as a function of temperature at 0 and 30 GPa, in comparison with the experimental data [41].

The comparison of the calculated specific heat to that predicted by the Debye model leads to the concept of the temperature dependent Debye temperature. The obtained Debye temperature at the high temperature limit is 550 K, which is consistent with the reported experimental result of 577 K [48]. In addition, the thermal conductivity ĸL contributed by the lattice can be estimated by the Debye equation [49]: where C is the volumetric heat capacity, λ is the mean free path of phonons, assumed as the average distance between the guest atoms, Vm is the velocity of sound which can be derived from Debye temperature [50]. In this way, the lattice thermal conductivity is given as 1.64 W m−1·K−1 which is comparable to other reported room temperature thermal conductivities of type-I silicon based clathrate compound [48].

4. Conclusions

In summary, we have given a detailed study of the vibrational and thermal properties of type-I calthrate K8Si46. By employing the linear-response method within the density functional perturbation theory (DFPT), we obtained the phonon dispersion relation in K8Si46 crystal as well as phonon density of states which contained detail information about the vibrational properties of both guest atoms and framework lattice. Under a pressure of 40 GPa, the frequencies from the motions of framework silicon atoms become imaginary around the M symmetry point suggests the instability of K8Si46 crystal. The determined phonon instable pressure by lattice dynamic simulation is larger than the transition pressure observed experimentally. However, actual phase transition may be driven by multi-mechanisms which are expected to reduce the phonon instable pressure. By using quasi-harmonic approximation, the thermal properties of K8Si46 including heat specific, thermal expansion coefficient and Grüneisen parameters are predicted. From a comparison of Na8Si46, we find that the thermal properties of these type-I calthrates are very insensitive to the different nature of guest atoms in the same group of the periodic table.