Effects of van der Waals interactions on the phonon transport properties of tetradymite compounds

During the past few decades, tetradymites have attracted much attention in the field of thermoelectric (TE) materials and topological insulators (TIs) due to the existence of heavy elements and small band gaps [1]. The well-known examples are the binary tetradymites Bi₂Te₃, Bi₂Se₃, Sb₂Te₃ and Sb₂Se₃, which are found to be well-established TE materials and their electronic and transport properties have been subject to numerous studies, both experimentally and theoretically [2–9]. Meanwhile, all these four compounds have been confirmed to be three-dimensional (3D) TIs with similar band inversion caused by the spin–orbit coupling, and feature inherent gapless surface states which behave as a single Dirac cone in the Brillouin zone [10–13].

Although extensive experimental and theoretical studies have been performed on the tetradymites, most of them are focused on the electronic and topological properties. The investigations of their phonon transports are less known. By performing classical molecular dynamics (CMD) simulations, Huang and Kaviany [14] firstly established a set of complicated interatomic potentials to calculate the lattice thermal conductivity of Bi₂Te₃. Using the same approach, Qiu and Ruan [15] developed a simpler two-body potentials for Bi₂Te₃ by fitting the results from density functional calculations. It should be mentioned that the precision of CMD strongly depends on appropriate selection of interatomic potentials [16], which may lead to certain discrepancy in the calculated lattice thermal conductivity compared with the experimentally measured value [17]. On the other hand, Sosso et al [18] theoretically predicted the Raman, infrared ray and phonon spectra of crystalline Sb₂Te₃. Although their results are overall consistent with the experimental data, the lack of appropriate exchange and correlation functionals tends to underestimate the frequency of longitudinal acoustic phonons. In principle, due to the presence of van der Waals (vdW) interactions in the tetradymites with layered structure, one should include vdW functionals for high-level density functional theory (DFT) calculations. However, the effects of such weak interactions are generally assumed to be neglectable, and standard DFT calculations without consideration of vdW functional can well reproduce/predict most physical results for many decades. Only in recent years, it was found that the seemingly weak vdW forces may have real effects on certain physical properties of some layered systems. For instance, using the GW method with optB86b-vdW functional, Cheng et al predicted a band gap of 157 meV for the Bi₂Te₃ that is much closer to the experimental value [19]. Surprisingly, the rhombohedral Sb₂Se₃ was proven to be also topologically non-trivial if the vdW interactions are explicitly taken into account [13, 20]. All these findings make it reasonable to expect that the intrinsic vdW may also play an important role in precisely determining other fundamental properties of tetradymites such as the phonon transport, which is so far less known to the best of our knowledge.

In this work, using four binary tetradymites (Bi₂Te₃, Bi₂Se₃, Sb₂Te₃ and Sb₂Se₃) as prototypical examples, we give a first-principles study on their phonon transport properties, and the results with and without consideration of the vdW forces are both discussed for comparison. It is surprising to find that the seemingly weak interlayer interactions actually play a very important role in predicting their lattice thermal conductivities, which are obviously higher than those obtained without vdW. The underlying physical mechanisms are explored in terms of the structural, harmonic and anharmonic properties of such layered systems.

The optimized structures of the four binary tetradymites can be obtained by performing DFT calculations [21–23], where the projector-augmented wave [24, 25] method is adopted as implemented in the Vienna ab initio simulation package [26]. The exchange–correlation energy is in the form of Perdew–Burke–Ernzerhof (PBE) with the generalized gradient approximation [27]. The energy cutoff is set to 500 eV and the Brillouin zone is sampled with an 11 × 11 × 11 Monkhorst–Pack k meshes [28]. To include the effect of vdW, the so-called optB86b-vdW [13] is selected after carefully checking various functional forms.

For the phonon properties, we use the finite displacement method to derive the harmonic and anharmonic interatomic force constants (IFCs) by using 4 × 4 × 4 and 3 × 3 × 3 supercells, respectively. The codes are embedded in the so-called PHONOPY [29] package and the THIRDORDER.PY script. Besides, the fourth nearest neighbors are imposed to achieve converged third-order IFCs. The thermal transport properties can be computed by solving the phonon Boltzmann transport equation, which is implemented in the ShengBTE program [30]. A fine q -mesh of 25 × 25 × 25 is used to obtain converged lattice thermal conductivity.

As shown in figure 1, the binary tetradymites Bi₂Te₃, Bi₂Se₃, Sb₂Te₃ and Sb₂Se₃ all have the same rhombohedral crystal structure with space group of ${D}_{3d}^{5}$ ($R\overline{3}m$). The primitive cell consists of five atoms, and it is usually convenient to show the layered structure using the hexagonal unit cell with 15 atoms. The system can be viewed as five atomic layers arranged along the c direction in the sequence of Te(Se)1–Bi(Sb)–Te(Se)2–Bi(Sb)–Te(Se)1 held by covalent bonds, known as the quintuple layers (QLs). Unlike the strong coupling between two atomic layers inside one QL, the vdW interactions between adjacent QLs are much weaker. However, it should be noted that the vdW is a dominant factor in correctly determining the interlayer spacing (d_eqm), which was previously found to have pronounced effects on the electronic and topological properties of tetradymites [13, 18, 19]. As summarized in table 1, the lattice parameters (a) predicted by the standard DFT calculations with PBE functional tend to be overestimated compared with the experimental values [31–34]. Such a problem can be effectively addressed by considering the vdW functionals, where the form of optB86b-vdW was confirmed to be the best that can well reproduce the lattice parameters of some tetradymites [13]. It is thus reasonable to choose such a functional exclusively to predict the real crystal structures of all the four binary tetradymites. Indeed, we see from table 1 there is obvious reduction of the interlayer spacing when the vdW is considered, which also leads to slightly decreased volume of the primitive cell.

Zoom In Zoom Out Reset image size

It is interesting to check if the weak vdW interactions could also affect the thermal properties of these tetradymites. Figure 2 shows the phonon dispersion relations of the Bi₂Te₃, Bi₂Se₃, Sb₂Te₃ and Sb₂Se₃. Since there are five atoms in the primitive cell of each system, we observe 15 phonon branches where three are the acoustic modes (including LA, TA, and ZA) and 12 are the optical modes (OP). There is no imaginary frequency which suggests that all the four binary tetradymites are dynamically stable, no matter the vdW functional is included or not. Due to the larger mass of constituent atoms, the maximum phonon frequencies of those four compounds are all smaller than 210 cm⁻¹. Among them, the Bi₂Te₃ exhibits the lowest value (140 cm⁻¹) and the Sb₂Se₃ has the highest (205 cm⁻¹), which is negatively correlated with their atomic mass as generally assumed. However, we find that the Bi₂Se₃ has a relatively higher cutoff frequency (175 cm⁻¹) compared with the Sb₂Te₃ (165 cm⁻¹), despite the former possess a larger total atomic mass. It is thus reasonable to expect that there exists stronger interatomic bonding of Bi–Se than that of Sb–Te [1, 35]. Furthermore, we find some hybridizations between acoustic and low-frequency optical branches of all the four compounds, which suggests that the inherent phonon–phonon scatterings are more likely to occur and the systems could have lower thermal conductivities [36]. On the other hand, we see that the dispersion relations with and without the consideration of vdW functional are almost the same for the optical branches. However, the presence of vdW tends to increase the phonon group velocity of the acoustic branches along the Γ–Z direction so that a larger lattice thermal conductivity can be expected. In another word, such effect should become more pronounced along the z direction, as discussed in the following.

Zoom In Zoom Out Reset image size

Table 2 summarizes the calculated lattice thermal conductivities (κ_l) of the binary tetradymites Bi₂Te₃, Bi₂Se₃, Sb₂Te₃ and Sb₂Se₃. For comparison, the results without the vdW functional are also included. It should be mentioned that in the present work, we focus on the effects of the vdW interactions on the phonon transport and only consider the lattice thermal conductivity at 300 K (the effect of thermal expansion was previously investigated by Hellman et al using temperature-dependent IFCs [37]). We see that the in-plane thermal conductivities are obviously larger than those along the out-of-plane direction, which is generally found for layered structures [38–40]. It should be mentioned that, if the vdW is explicitly considered in the calculation, the κ_l of the Bi₂Te₃, Bi₂Se₃ and Sb₂Te₃ agree well with those found previously [1, 37, 41]. On the other hand, including the vdW leads to increased thermal conductivities for all the four investigated systems, especially for those along the out-of-plane direction, which is consistent with obviously enhanced phonon dispersion relations along the Γ–Z direction discussed above. For example, the thermal conductivities of Sb₂Se₃ increase respectively by 39% and 84% along the in-plane and out-of-plane directions when the vdW functional is explicitly included in the calculations. Such a surprising finding indicates that the generally assumed weak vdW forces are actually of vital importance in predicting the phonon transport properties of layered systems, and a deep understanding of the physical mechanisms involved is quite necessary.

We now focus on the Sb₂Se₃ where the effect of vdW is the strongest among the four investigated binary tetradymites. Table 3 summarizes the contributions of each phonon branch to the lattice thermal conductivity of the Sb₂Se₃ calculated with and without consideration of the vdW, where we see that the heat transport is mainly governed by the three acoustic branches, especially for the out-of-plane direction. Together with figure 2(d), such a finding can be used to qualitatively explain the obvious increase of the lattice thermal conductivity when the vdW is taken into considerations. To go into the details, we know that the lattice thermal conductivity can be expressed as ${\kappa }_{\mathrm{l}}=\frac{1}{{\Omega}}{\sum }_{i}{c}_{i}{v}_{i}^{2}{\tau }_{i}$, where Ω, c, v and τ are the volume of the primitive cell, the heat capacity, the group velocity, and the relaxation time, respectively. The summation is over all the phonon modes. As shown in figure 2(d), the phonon dispersion relations of three acoustic branches along the Γ–Z direction are obviously enhanced when the vdW functional is included, which means that they have a larger group velocity. Indeed, our calculated v increases from 954 m s⁻¹ (747 m s⁻¹) to 2490 m s⁻¹ (1875 m s⁻¹) for the LA (TA) mode at the Γ point. Table 4 summarizes the average values of the relaxation time ($\overline{\tau }$) for the LA, TA, ZA and OP phonon modes. It is found that considering the vdW leads to a significantly larger $\overline{\tau }$, especially for the LA and TA modes. Indeed, the overall average relaxation time is calculated to be 6.06 ps and 3.13 ps for the system with and without vdW forces, respectively. All these observations suggest that the effect of vdW is to increase both the phonon group velocity and the relaxation time so that a larger lattice thermal conductivity of Sb₂Se₃ can be expected.

Table 4. Average values of the phonon relaxation time ($\overline{\tau }$, in unit of ps) and Grüneisen parameter ($\overline{\gamma }$) for each phonon branch of the Sb₂Se₃, calculated with and without the consideration of vdW interactions.

Variable	vdW	LA	TA	ZA	OP
$\overline{\tau }$	With	8.43	8.21	5.50	2.11
	W/o	3.94	3.03	4.57	0.97
$\overline{\gamma }$	With	0.97	0.66	1.58	1.92
	W/o	4.77	2.49	2.67	2.10

In principle, the phonon relaxation time of a given system is closely related to its anharmonicity, which can be characterized by the so-called Grüneisen parameter (γ). As known [42–44], a smaller γ means a weaker phonon scattering strength and therefore a longer relaxation time which consequently leads to a higher lattice thermal conductivity. Indeed, we see from figures 3(a) and (b) that the Sb₂Se₃ with vdW considered exhibits obviously smaller γ, especially for the LA, TA, and ZA branches which are previously confirmed to be the dominant contribution to the lattice thermal conductivity. As can be also found from table 4, the average Grüneisen parameter ($\overline{\gamma }$) of these three acoustic branches shows a significant reduction when the vdW functional is explicitly included in the calculations. In particular, the calculated total γ_tot is also decreased from 2.29 to 1.69, suggesting the reduced phonon anharmonicity and thus enhanced lattice thermal conductivity. As a further evidence, we plot in figures 3(c) and (d) the energy change (ΔE) with respect to the relative lattice constant (a/a₀) of the Sb₂Se₃. For comparison, the results with and without vdW are both shown. It is found that ΔE can be well fitted by a cubic polynomial, where the coefficient of the third-order term is respectively 72 and 92 for the calculations with and without vdW. That is, the system with vdW explicitly included exhibits relatively smaller anharmonicity as discussed above. It is noted that similar results can be found in other binary tetradymites, which are summarized in supplementary table 1 and figure 1 (https://stacks.iop.org/NJP/23/083002/mmedia). It should be noted that the inclusion of the vdW functional also enhances the in-plane thermal conductivity, which is more pronounced for the Sb₂Se₃. The reason is that the vdW force is not completely along the out-of-plane direction, and its in-plane component could also affect the corresponding phonon transport. To have a better understanding, we plot in supplementary figure 2 the atomic displacement parameters along the in-plane direction for the Se atom at the QL edge of Sb₂Se₃ (i.e. the Se1 in figure 1). We see that the in-plane vibration of the Se atom is obviously smaller when the vdW force is included, which suggests a lower anharmonicity and thus a higher thermal conductivity. Of course, the major effect of the vdW force on the phonon transport is along the out-of-plane direction.

Zoom In Zoom Out Reset image size

Before concluding our discussions, it should be mentioned that the effect of vdW on the lattice thermal conductivity can be also understood by a simple physical picture. For the layered structures such as the Sb₂Se₃, the coexistence of strong covalent bonding inside the QL and much weaker vdW force between neighboring QLs leads to inherently strong phonon anharmonicity, which will be diminished at smaller interlayer distance when the seemingly weak vdW force is presented. As a consequence, lower lattice thermal conductivity can be expected.

In conclusion, we demonstrate by DFT calculations and Boltzmann transport theory that the weak vdW interactions actually have pronounced effects on the phonon transport properties of the layered tetradymites. The physical origin is that considering the vdW leads to obviously reduced interlayer distance, which can effectively weaken the strong phonon anharmonicity of the system, as confirmed by the decrease of the Grüneisen parameter, as well as the increase of the group velocity and relaxation time. As a consequence, all the four investigated tetradymites exhibit larger lattice thermal conductivities in the presence of vdW, especially for those along the out-of-plane direction. Our work emphasizes the very important role of the usually ignored vdW interactions in accurately predicting the lattice thermal conductivity of the tetradymite compounds, and should be also applied in other layered systems such as transition metal dichalcogenides and various heterostructures of two-dimensional materials.

We thank financial support from the National Natural Science Foundation of China (Grant Nos. 51772220 and 62074114). The numerical calculations in this work have been done on the platform in the Supercomputing Center of Wuhan University.

The data that support the findings of this study are available upon reasonable request from the authors.