Anisotropic phonon thermal transport in two-dimensional layered materials

Abstract

Two-dimensional layered materials (2DLMs) have attracted growing attention in optoelectronic devices due to their intriguing anisotropic physical properties. Different members of 2DLMs exhibit unique anisotropic electrical, optical, and thermal properties, fundamentally related to their crystal structure. Among them, directional heat transfer plays a vital role in the thermal management of electronic devices. Here, we use density functional theory calculations to investigate the thermal transport properties of representative layered materials: β-InSe, γ-InSe, MoS₂, and h-BN. We found that the lattice thermal conductivities of β-InSe, γ-InSe, MoS₂, and h-BN display diverse anisotropic behaviors with anisotropy ratios of 10.4, 9.4, 64.9, and 107.7, respectively. The analysis of the phonon modes further indicates that the phonon group velocity is responsible for the anisotropy of thermal transport. Furthermore, the low lattice thermal conductivity of the layered InSe mainly comes from low phonon group velocity and atomic masses. Our findings provide a fundamental physical understanding of the anisotropic thermal transport in layered materials. We hope this study could inspire the advancement of 2DLMs thermal management applications in next-generation integrated electronic and optoelectronic devices.

Graphical abstract

Keywords

thermal conductivity / two-dimensional layered materials / first-principles calculation / Boltzmann transport theory

1 Introduction

Since the successful exfoliation of graphene in 2004, two-dimensional layered materials (2DLMs) have been intensively investigated due to their fascinating physical properties [1–7], such as high carrier mobility, tunable band gap, outstanding photoelectric characteristics, and thermal stability [8–14]. It is expected to become one of the next-generation optoelectronic materials of the most prominent families of 2DLMs: hexagonal boron nitride (h-BN) [15], transition metal dichalcogenides (TMDCs) [16], and post-transition metal chalcogenides (PTMCs) [17–20]. 2DLMs are bonded by strong covalent bonds in the in-plane and coupled by the weak van der Waals (vdW) forces in the out-of-plane direction with strong structural anisotropy. Therefore, layered materials have interesting anisotropic physical properties in terms of electrical [21], optical [22, 23], magnetic [24], and thermal transport properties [25–27]. For example, 2DLMs with high lattice thermal conductivity (κ_l) in the in-plane direction can be used as a heat-conducting material in electronic and optoelectronic devices [28]. On the other hand, the low κ_l of 2DLMs along the out-of-plane direction can be exploited in a variety of applications in thermoelectricity [29, 30].

With the rapid development of 2D electronic devices toward high power density and nanoscale [31–34], it is necessary to provide them with efficient thermal management strategies to achieve directional heat transfer [35]. 2DLMs can provide an ideal material platform for creating such applications due to their special heat transfer mechanism and anisotropic thermal transport properties. For bulk h-BN, the intrinsic anisotropy ratio at room temperature varies from 69 to 111 based on time-domain thermoreflectance measurement [26]. In contrast, in PTMCs, such as β-InSe, the measured anisotropy ratio is small, spanning the range of 8 to 14 [25]. Additionally, the phonon thermal transport properties in the out-of-plane direction of 2DLMs can be effectively controlled by their weak interlayer interactions. For example, Kim et al. [36] modulated the thermal transport of MoS₂ through interlayer rotation, which resulted in a high anisotropy ratio of 900 at room temperature. Chen et al. [37] applied a series of cross-plane strains and reported strongly tunable anisotropy ratio in the range of 5 to 250 for MoS₂. However, in the literature, there is no systematic theoretical study of the thermal transport properties of 2DLMs with different structural units, especially anisotropic effects. In addition, the underlying physical mechanism leading to a wide range of anisotropic thermal transport in 2DLMs is also unexplored. Thus, it is highly desirable to screen the anisotropic thermal properties of 2DLMs for next-generation integrated optoelectronic devices.

Here, through first-principles calculations and lattice dynamics, we explore the thermal transport properties of four 2DLMs, including β-InSe, γ-InSe, MoS₂, and h-BN. In these compounds, the harmonic parameters (i.e., phonon group velocity, phonon frequency, and atomic mass) play a dominant role in determining κ_l. In addition, we found that these compounds exhibit various anisotropic thermal transport properties along different directions (mainly originated from different phonon group velocities). Our work provides an in-depth theoretical understanding of the anisotropic thermal transport in 2DLMs and can guide future experimental work for the practical thermal management application of 2DLMs.

2 Computational methodology

The crystal structures of β-InSe, γ-InSe, MoS₂, and h-BN were optimized utilizing the plane-wave projector augmented-wave (PAW) [38] method implemented in the Vienna Ab-initio Simulation Package (VASP) [39]. We employed the Perdew−Burke−Ernzerhoff (PBE) scheme to describe the exchange-correlation functional within generalized gradient approximation (GGA) [40]. The plane-wave kinetic energy cutoff was set to 520 eV, and the uniform Monkhorst–Pack k-point meshes with a grid spacing of 2π × 0.03 Å⁻¹ were employed in the Brillouin zone (BZ) sampling of the considered systems. The geometry was totally relaxed until the atomic force and energy values reached 10⁻⁴ eV/Å and 10⁻⁸ eV, respectively. Due to the van der Waals weak interaction in layered materials, we tested the optB88-vdW and optB86b-vdW functionals to obtain equilibrium lattice constants [41]. The optimized lattice parameters for bulk β-InSe, γ-InSe, MoS₂, and h-BN obtained with the optB86b-vdW functional are in close agreement with the experimental results (Table S1), with a percent error of less than 3%. Therefore, we will employ the optB86b-vdW functional in our further calculations.

The κ_l was obtained by iteratively solving the linearized Boltzmann−Peierls transport equation (BTE) of phonons by summing all the phonon modes in the BZ by the following equation [42]:

(1){\kappa }_l^{\alpha \beta }={\displaystyle \frac{1}{NV}\sum _{\lambda }{C}_{\lambda }{v}_{\lambda }^{\alpha }{v}_{\lambda }^{\beta }{\tau }_{\lambda }.}

Here, N: total number of q points, V: volume of the unit cell, \lambda: is the phonon mode, C_{\lambda }: phonon mode-specific heat capacity, v_{\lambda }^{\alpha } /v_{\lambda }^{\beta }: represents the phonon group velocity along the α and β directions, {\tau }_{\lambda }: the anharmonic phonon relaxation time, which can be determined by combining the anharmonic scattering rate ({\mathrm{\Gamma }}_{{\lambda \lambda }^{{}^{\prime }}{\lambda }^{{}^{\prime }{}^{\prime }}}^{\pm }) and the isotopic scattering rate ({\mathrm{\Gamma }}_{{\lambda \lambda }^{{}^{\prime }}}),

(2){\displaystyle \frac{1}{{\tau }_{\lambda }}}={\displaystyle \sum _{{\lambda }^{{}^{\prime }}{\lambda }^{{}^{\prime }{}^{\prime }}}^{+}{\Gamma }_{{\lambda \lambda }^{{}^{\prime }}{\lambda }^{{}^{\prime }{}^{\prime }}}^{+}+\sum _{{\lambda }^{{}^{\prime }}{\lambda }^{{}^{\prime }{}^{\prime }}}^{-}\frac{1}{2}{\Gamma }_{{\lambda \lambda }^{{}^{\prime }}{\lambda }^{{}^{\prime }{}^{\prime }}}^{-}+\sum _{{\lambda }^{{}^{\prime }}}{\Gamma }_{{\lambda \lambda }^{{}^{\prime }}}.}

Here, the anharmonic scattering ({\mathrm{\Gamma }}_{{\lambda \lambda }^{{}^{\prime }}{\lambda }^{{}^{\prime }{}^{\prime }}}^{\pm }) can be expressed as [43–45]

(3)\begin{array}{rl}{\Gamma }_{{\lambda \lambda }^{{}^{\prime }}{\lambda }^{{}^{\prime }{}^{\prime }}}^{\pm }=& {\displaystyle \frac{\hslash \pi }{8N}\left\{\begin{array}{c}2\left(f_{{\lambda }^{{}^{\prime }}}-f_{{\lambda }^{{}^{\prime }{}^{\prime }}}\right)\\ f_{{\lambda }^{{}^{\prime }}}+f_{{\lambda }^{{}^{\prime }{}^{\prime }}}+1\end{array}\right\}}\\ & \cdot \frac{\delta \left({\omega }_{\lambda }\pm {\omega }_{{\lambda }^{{}^{\prime }}}-{\omega }_{{\lambda }^{{}^{\prime }{}^{\prime }}}\right)}{{\omega }_{\lambda }{\omega }_{{\lambda }^{{}^{\prime }}}{\omega }_{{\lambda }^{{}^{\prime }{}^{\prime }}}}\left|{\mathrm{V}}_{{\lambda \lambda }^{{}^{\prime }}{\lambda }^{{}^{\prime }{}^{\prime }}}^{\pm }\right|.\end{array}

The terms in the upper and lower curly brackets define the absorption and emission of phonons, respectively. The terms \left|V_{{\lambda \lambda }^{{}^{\prime }}{\lambda }^{{}^{\prime }{}^{\prime }}}^{\pm }\right| describe the scattering matrix elements, which is given by [45]

(4)V_{{\lambda \lambda }^{{}^{\prime }}{\lambda }^{{}^{\prime }{}^{\prime }}}^{\pm }={\displaystyle \sum _{i\in u.c.}\sum _{j,k}\sum _{\alpha \beta \gamma }{\Phi }_{ijk}^{\alpha \beta \gamma }\frac{e_{\lambda }^{\alpha }\left(i\right)e_{p^{{}^{\prime }},\pm q^{{}^{\prime }}}^{\beta }\left(j\right)e_{p^{{}^{\prime }{}^{\prime }},-q^{{}^{\prime }{}^{\prime }}}^{\gamma }\left(k\right)}{\sqrt{M_i{M}_j{M}_k}}.}

The weighted phase space is defined as the number of channels in which three phonons can simultaneously satisfy the conservation laws for quasi-momentum and energy [46]:

(5)W_{\lambda }^{\pm }={\displaystyle \frac{1}{2N}}{\displaystyle \sum _{{q,q}^{{}^{\prime }}}\left\{\begin{array}{c}2(f_{{\lambda }^{{}^{\prime }}}-f_{{\lambda }^{{}^{\prime }{}^{\prime }}})\\ f_{{\lambda }^{{}^{\prime }}}+f_{{\lambda }^{{}^{\prime }{}^{\prime }}}+1\end{array}\right\}\frac{\delta ({\omega }_{\lambda }\pm {\omega }_{{\lambda }^{{}^{\prime }}}-{\omega }_{{\lambda }^{{}^{\prime }{}^{\prime }}})}{{\omega }_{\lambda }{\omega }_{{\lambda }^{{}^{\prime }}}{\omega }_{{\lambda }^{{}^{\prime }{}^{\prime }}}}.}

The harmonic phonon transport properties and second-order interatomic force constants (IFCs) were obtained by Phonopy code [47] using the finite displacement method with 5 × 5 × 2 supercell for h-BN, and 4 × 4 × 1 supercells for β-InSe, γ-InSe, and MoS₂, respectively. The ionic dielectric tensor and Born effective charge were calculated using density functional perturbation theory (DFPT) [48], which corrects the longitudinal and transverse optical (LO-TO) splitting of phonon dispersion caused by long-range dipole-dipole interactions (see Table S2 for further details). The anharmonic third-order IFCs of h-BN, MoS_2, and InSe were simulated with thirdorder.py code [49] using a 3 × 3 × 1 supercell, and the range of the considered interatomic three-body interactions was set to 5.99, 7.04, and 8.19 Å, as they have proved very successful according to previous reports [37, 50, 51]. To ensure sufficiently convergent κ_l, we employed the ShengBTE code [49] with q-point meshes of 24 × 24 × 5, 29 × 29 × 4, 25 × 25 × 5, and 30 × 30 × 10 in solving the phonon BTE for β-InSe, γ-InSe, MoS₂, and h-BN, respectively, as illustrated in Fig. S1.

3 Results and discussion

3.1 Anisotropic structure and ratio

The relaxed crystal structure and interlayer differential charge densities of bulk β-InSe, MoS₂, and h-BN are shown in Fig.1(a), whereas the crystal structure of γ-InSe is provided in Fig. S2(a). These materials crystallize in the hexagonal honeycomb structure indexed in the space group P6₃/mmc, except for γ-InSe (R3mH). The intralayer consists of strongly bonded 2D atomic sheets, and the interlayer is coupled to each other through relatively weak vdW interaction, both of which result in structural anisotropy. Furthermore, these structures differ greatly in the in-plane arrangement of surface atoms. Bulk InSe and MoS₂ exhibit buckling configurations with four (Se−In−In−Se) and three (S–Mo–S) atomic layers, respectively, while bulk h-BN has a planar geometry consisting of diatomic B−N. In addition, γ-InSe exhibits ABC stacking, whereas β-InSe, MoS₂, and h-BN have AB stacking order. Moreover, the interlayer differential charge densities can directly reflect the character of the interlayer interaction. Clearly, in the β-InSe and MoS₂ systems, a large number of electrons are localized in the interlayer region, indicating the existence of strong interlayer coupling and consistent with previous results reported elsewhere [52–54]. On the contrary, the charge distribution in the interlayer spacings of h-BN is negligible, indicating relatively weak interlayer coupling.

Next, we consider the trend of the lattice thermal conductivities (κ_l) obtained by iteratively solving the linearized Boltzmann−Peierls transport equation of phonons. The anisotropy ratio is defined as κ_in/κ_out, where κ_in and κ_out are the calculated κ_l along the in-plane and out-of-plane directions, respectively. The κ_l ratio obtained at 300 K is presented in Fig.1(b). It is clear that the anisotropy ratios of β-InSe, γ-InSe, MoS₂, and h-BN are significantly different from each other. In particular, the anisotropy ratios of β-InSe and γ-InSe are 10.4 and 9.4, respectively, indicating that different stacking structures have a less direct influence on κ_l; therefore, the thermal transport coefficients of γ-InSe are shown in Figs. S2(b–d). Subsequently, the anisotropy ratio of MoS₂ is 64.9, and h-BN is 107.7, which are the highest values reported in our work. The anisotropy ratio is mainly determined by the κ_in and κ_out. Therefore, accurate prediction of κ_in and κ_out are the basic premises for quantitatively evaluating thermal transport anisotropy.

3.2 Anisotropic lattice thermal conductivity

Fig.2 compares the variation of the κ_in and κ_out with temperature for β-InSe, MoS₂, and h-BN. The theoretical results and available experimental data of the κ_in and κ_out for β-InSe, MoS₂, and h-BN at 300 K are included in Table S3. Our results are in good agreement with the previous experiments, showing the validity of our calculations. We find that the quantitative relationship between κ_in and κ_out at 300 K follows the order κ_l (β-InSe) < κ_l (MoS₂) < κ_l (h-BN). In particular, MoS₂ and h-BN have relatively high κ_in, over 11 and 50 times higher than that of InSe. For κ_out, both MoS₂ and h-BN show 1.5 and 4.9 times larger κ_l than InSe, and the gap is far less than κ_in. In the whole temperature range, both the in-plane and out-of-plane components decrease with increasing temperature and exhibit strong anisotropy. Interestingly, the relative magnitude of the anisotropy ratio is weakly dependent on temperature.

3.3 Anharmonic phonon lifetime

We further investigated the phonon lifetime, phonon spectrum, scattering phase space, and phonon group velocity to fully understand the anisotropic thermal transport in InSe, MoS₂, and h-BN. All of these are closely related to the phonon thermal transport properties, as simplified by Eq. (1). The relationship between phonon lifetime ({\tau }_{\lambda }) and frequency for β-InSe, MoS₂, and h-BN at 300 K is depicted in Fig.3(a−c). It can be seen that {\tau }_{\lambda } of β-InSe, MoS_2, and h-BN in the in-plane direction is higher than that of the out-of-plane component, showing the same direction-dependent anisotropy similar to κ_l. In addition, we also calculated the average {\tau }_{\lambda } of all phonon modes, and the obtained values are 5.3, 8.3, and 3.6 ps for β-InSe, MoS₂, and h-BN, respectively. However, the effect of {\tau }_{\lambda } is not good enough to explain the significant difference in κ_l.

3.4 Anisotropic phonon dispersion

Subsequently, we turn our attention to the phonon dispersion spectrum and phonon group velocity related to the harmonic properties of phonons. Accurate phonon dispersion is critical to determine the required group velocity and the number of allowed three-phonon scattering channels for κ_l calculations. Hence, we tested the effect of supercell size on the phonon spectrum (Fig. S3) and selected large supercells of size 5 × 5 × 2 for h-BN, and 4 × 4 × 1 for β-InSe, γ-InSe, and MoS₂ to encompass a wider range of second-order interatomic force constants to describe the phonon dispersion spectrum better.

The phonon dispersion curves of β-InSe, MoS₂, and h-BN are shown in Fig.4(a) and (c). The absence of imaginary modes indicates the dynamic stability of the system. The cumulative κ_l as a function of phonon frequency reveals that the κ_l is dominated by the low-frequency acoustic phonon branch (see Fig. S4). Therefore, here we focus on the three acoustic heat-carrying phonon modes, i.e., the in-plane longitudinal acoustic (LA), transverse acoustic (TA), and the out-of-plane flexural acoustic (ZA) branch [Fig.4(a, c)]. The bond strength in the in-plane direction of 2DLMs is much stronger than the weak interlayer interaction (vdW) in the out-of-plane direction. Accordingly, high-frequency phonon modes are generated in the in-plane, and low-frequency phonon modes are generated in the out-of-plane direction.

We further investigated the difference in phonon frequencies of the in-plane and out-of-plane directions of InSe, MoS₂, and h-BN in Fig.4(a) and (c). In the in-plane direction, taking the LA phonon branch as an example, at point M, the corresponding frequencies of β-InSe, MoS₂, and h-BN are 1.7, 6.9, and 34.3 THz, respectively. Such a high difference in frequency may be due to different interatomic bonding. In addition, we noticed that there are also differences in the shapes of dispersion curves; the acoustic phonons of β-InSe have significant cross-overlaps with the optical phonons, followed by MoS₂. Meanwhile, the acousto−optic coupling of h-BN is small. This is due to the large number of atoms in the primitive cell of InSe, resulting in a rapid increase in the proportion of optical phonon branches. Specifically, there are 8 atoms per unit cell of β-InSe, which results in 21 optical branches, while the primitive cells of MoS₂ and h-BN contain 6 and 4 atoms; thus, the respective phonon spectra contains 15 and 9 optical branches. Therefore, suppressed phonon frequencies and strong acoustic−optic coupling will likely enhance the number of scattering channels available for phonon−phonon interactions in β-InSe, which favors low κ_l [56]. Fig.4(b) shows the value of the scattering phase space conforms to β-InSe >MoS₂ > h-BN, which is in contrast to the numerical trend of κ_in.

On the other hand, the shape of the dispersion curves for the out-of-plane direction of β-InSe, MoS₂, and h-BN are very similar, as shown in Fig.4(c). Only the magnitude of the phonon frequency is different. At point A, the LA acoustic branch frequencies corresponding to β-InSe, MoS₂, and h-BN are 0.7, 1.2, and 2.4 THz, respectively. It is worth mentioning that though the phonon frequencies in the out-of-plane direction are different, the numerical difference between the three structures is much less pronounced than that in the in-plane direction. This corresponds to the fact that the numerical difference in the phonon group velocity becomes less prominent than in the in-plane direction. We attribute the frequency difference in the out-of-plane direction primarily to the average atomic mass (Table S4). Specifically, β-InSe has a lower interlayer phonon frequency due to its heavier average atomic mass, followed by MoS₂. In contrast, h-BN has a lighter average atomic mass, leading to high interlayer phonon frequencies. In addition, the scattering phase space [shown in Fig.4(d)] also satisfies the quantitative relationship of β-InSe > MoS₂ > h-BN, further supporting the predicted κ_out in these materials.

3.5 Anisotropic phonon group velocity

The slope of the phonon dispersion curve gives information about the group velocity \left(v_{\lambda }\right), and the square of its value ({v_{\lambda }}^2) is essential in determining the κ_l [see Eq. (1)]. We first present the v_{\lambda } as a function of frequency for β-InSe, MoS₂, and h-BN in Fig. S5. As mentioned earlier, the v_{\lambda } along the in-plane direction is significantly larger than the out-of-plane direction over the entire frequency range, finding agreement with the fact that the phonon frequency along the Γ−M path is more dispersed than that along the Γ−A path. Moreover, in Fig.5(a) we compare the axial average v_{\lambda } of β-InSe, MoS₂, and h-BN to understand the difference in anisotropic κ_l directly. For the in-plane direction, the calculated v_{\lambda } of β-InSe is 309 m/s, which is lower than that of MoS₂ (753 m/s), and much lower than that of h-BN (3007 m/s). Furthermore, the large difference in group velocity in the in-plane direction is well understood from the elastic constant. We calculated the elastic constants based on the stress-strain relationship [57], and the ultimate results are summarized in Table S4. The elastic constant of h-BN is 3.8 and 13.9 times larger than that of MoS₂ and β-InSe, whereas the elastic constant of MoS₂ is 3.6 times larger than that of β-InSe. The calculated results confirm our preliminary predictions that the low in phonon frequency in the in-plane direction is mainly attributed to the weak bonding strength, leading to the smaller phonon group velocity in β-InSe. For the out-of-plane direction, the average v_{\lambda } of β-InSe is 106 m/s, relatively smaller than MoS₂ (120 m/s) and h-BN (326 m/s).

We also present the anisotropy ratio relationship between κ_l and the square of average v_{\lambda } in Fig.5(b). We discovered that the anisotropy ratio of the square of average v_{\lambda } for β-InSe, MoS₂ and h-BN are 8.5, 51.6 and 85.1, respectively, contributing around 80 % of the total anisotropy ratio of the κ_l. This further justifies that the phonon group velocity plays an essential role in the direction dependence of κ_l. In addition, it can also be deduced from Fig.5(b) that the slope of the line fitted by the anisotropy ratio between κ_l and v_{\lambda } for all the considered structures is less than unity, which is due to the anisotropic {\tau }_{\lambda } explained above; together contributing towards the anisotropic κ_l. We find that v_{\lambda } is the main reason for the difference in the κ_in and κ_out of β-InSe, MoS₂, and h-BN. In addition, v_{\lambda } is also the dominant factor causing diverse anisotropic thermal transport in β-InSe, MoS_2, and h-BN. Finally, we also studied the anisotropic lattice thermal conductivity of β-InSe, γ-InSe, MoS₂, and h-BN as a function of sample size (Fig. S6). This is particularly useful for thermal conductivity measurements of experimentally limited thin-film samples ranging from nanometer to micrometer scale.

4 Conclusions and outlook

In summary, we employed ab-initio method and lattice dynamics to investigate the origin of the anisotropic κ_l of 2DLMs: InSe, MoS₂, and h-BN. We found that the κ_l of InSe, MoS_2, and h-BN have significant differences in the in-plane and out-of-plane directions, resulting in a wide range of anisotropy ratios (9.4−107.7). Then, the physical mechanism affecting the κ_l is elucidated from the perspective of phonon transport modes. The κ_l of these materials is mainly dominated by phonon frequency, phonon group velocity, and atomic mass. On the other hand, the anisotropic thermal transport mainly originates from the difference of phonon group velocity in the in-plane and out-of-plane direction, which is caused by the anisotropic phonon dispersion.

Our results reveal the physical mechanism of anisotropic thermal transport in 2DLMs and provide a routine method for the rational design of the anisotropic heat transport mechanism in 2DLMs. Three key features are important, including (i) crystal structure; (ii) harmonic parameters; and (iii) interlayer coupling. For example, Yuan et al. [58] reported the anisotropic thermal transport behaviour of layered InSe under hydrostatic pressure. With increasing pressure, the strong interlayer coupling may help to reduce the anisotropy ratio of InSe [25, 58], and the ratio of κ_in to κ_out decreases from 6.95 at 0 GPa to 1.26 at 8 GPa, exhibiting a transition from anisotropic to isotropic. Furthermore, it is relatively common to manipulate the anisotropy of thermal transport by tuning the interfacial phonon properties of vertically stacked heterojunctions in 2DLMs. For example, the results of molecular dynamics simulations by Ren et al. [59] showed that the thermal transport at the graphene/h-BN heterojunction interface could be effectively controlled by tuning the interlayer rotation angle between graphene and h-BN layers. Nobakht et al. [60] reported that by applying a pressure of 2.6 GPa to the graphene interlayer, the thermal resistance between the graphene/silicon interface could be reduced by 50 % without significantly changing the κ_in. In practical device applications, the substrate also significantly impacts the κ_l at the interface of 2DLMs. For example, Ni et al. [61] reported the average interlayer thermal resistance of suspended and supported few-layer graphene (FLG) using equilibrium molecular dynamics simulations. The results show that the FLG with SiO₂ substrate can be significantly improve the κ_out compared to suspended FLG. Studying the phonon transport of the substrate at the interface is also an important direction and is worthy of further exploration in future theoretical studies. To sum up, our work may drive a wide tuning range of anisotropic thermal transport in 2DLMs by tailoring their lattice thermal transport properties for high-tech thermal management applications, such as integrated electronics and optoelectronics.

References

[1] A. K. Geim, K. S. Novoselov. The rise of graphene. Nat. Mater., 2007, 6(3): 183

[2] J. C. Meyer, A. K. Geim, M. I. Katsnelson, K. S. Novoselov, T. J. Booth, S. Roth. The structure of suspended graphene sheets. Nature, 2007, 446(7131): 60

[3] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, A. A. Firsov. Electric field effect in atomically thin carbon films. Science, 2004, 306(5696): 666

[4] E. Scalise, M. Houssa, G. Pourtois, V. Afanas’ev, A. Stesmans. Strain-induced semiconductor to metal transition in the two-dimensional honeycomb structure of MoS₂. Nano Res., 2012, 5(1): 43

[5] D. Chiappe, E. Scalise, E. Cinquanta, C. Grazianetti, B. van den Broek, M. Fanciulli, M. Houssa, A. Molle. Two-dimensional Si nanosheets with local hexagonal structure on a MoS₂ surface. Adv. Mater., 2014, 26(13): 2096

[6] Z. Yang, W. Jie, C. H. Mak, S. Lin, H. Lin, X. Yang, F. Yan, S. P. Lau, J. Hao. Wafer-scale synthesis of high-quality semiconducting two-dimensional layered InSe with broadband photoresponse. ACS Nano, 2017, 11(4): 4225

[7] K. Zhang, Y. Feng, F. Wang, Z. Yang, J. Wang. Two-dimensional hexagonal boron nitride (2D-hBN): Synthesis, properties and applications. J. Mater. Chem. C, 2017, 5(46): 11992

[8] L. Wu, J. Shi, Z. Zhou, J. Yan, A. Wang, C. Bian, J. Ma, R. Ma, H. Liu, J. Chen, Y. Huang, W. Zhou, L. Bao, M. Ouyang, S. T. Pantelides, H. J. Gao. InSe/hBN/graphite heterostructure for high-performance 2D electronics and flexible electronics. Nano Res., 2020, 13(4): 1127

[9] G. W. Mudd, S. A. Svatek, L. Hague, O. Makarovsky, Z. R. Kudrynskyi, C. J. Mellor, P. H. Beton, L. Eaves, K. S. Novoselov, Z. D. Kovalyuk, E. E. Vdovin, A. J. Marsden, N. R. Wilson, A. Patanè. High broad-band photoresponsivity of mechanically formed InSe-graphene van der Waals heterostructures. Adv. Mater., 2015, 27(25): 3760

[10] D. Buckley, Z. R. Kudrynskyi, N. Balakrishnan, T. Vincent, D. Mazumder, E. Castanon, Z. D. Kovalyuk, O. Kolosov, O. Kazakova, A. Tzalenchuk, A. Patanè. Anomalous low thermal conductivity of atomically thin InSe probed by scanning thermal microscopy. Adv. Funct. Mater., 2021, 31(11): 2008967

[11] H. Arora, Y. Jung, T. Venanzi, K. Watanabe, T. Taniguchi, R. Hübner, H. Schneider, M. Helm, J. C. Hone, A. Erbe. Effective hexagonal boron nitride passivation of few-layered InSe and GaSe to enhance their electronic and optical properties. ACS Appl. Mater. Interfaces, 2019, 11(46): 43480

[12] B. Radisavljevic, A. Radenovic, J. Brivio, V. Giacometti, A. Kis. Single-layer MoS₂ transistors. Nat. Nanotechnol., 2011, 6(3): 147

[13] M. Brotons-Gisbert, D. Andres-Penares, J. Suh, F. Hidalgo, R. Abargues, P. J. Rodriguez-Canto, A. Segura, A. Cros, G. Tobias, E. Canadell, P. Ordejón, J. Wu, J. P. Martínez-Pastor, J. F. Sánchez-Royo. Nanotexturing to enhance photoluminescent response of atomically thin indium selenide with highly tunable band gap. Nano Lett., 2016, 16(5): 3221

[14] Y. Y. Hui, X. Liu, W. Jie, N. Y. Chan, J. Hao, Y. T. Hsu, L. J. Li, W. Guo, S. P. Lau. Exceptional tunability of band energy in a compressively strained trilayer MoS₂ sheet. ACS Nano, 2013, 7(8): 7126

[15] S.MoonJ. KimJ.ParkS.ImJ.Kim I.HwangJ. K. Kim, Hexagonal boron nitride for next‐generation photonics and electronics, Adv. Mater. 2204161 (2022)

[16] M. Wu, Y. Xiao, Y. Zeng, Y. Zhou, X. Zeng, L. Zhang, W. Liao. Synthesis of two‐dimensional transition metal dichalcogenides for electronics and optoelectronics. InfoMat, 2021, 3(4): 362

[17] S. Lei, L. Ge, S. Najmaei, A. George, R. Kappera, J. Lou, M. Chhowalla, H. Yamaguchi, G. Gupta, R. Vajtai, A. D. Mohite, P. M. Ajayan. Evolution of the electronic band structure and efficient photo-detection in atomic layers of InSe. ACS Nano, 2014, 8(2): 1263

[18] D. J. Late, B. Liu, J. Luo, A. Yan, H. S. S. R. Matte, M. Grayson, C. N. R. Rao, V. P. Dravid. GaS and GaSe ultrathin layer transistors. Adv. Mater., 2012, 24(26): 3549

[19] P. Hu, Z. Wen, L. Wang, P. Tan, K. Xiao. Synthesis of few-layer GaSe nanosheets for high performance photodetectors. ACS Nano, 2012, 6(7): 5988

[20] J. Xie, L. Zhang. Optical emission enhancement of bent InSe thin films. Sci. China Inf. Sci., 2021, 64(4): 140405

[21] T. Hu, H. Zhang, J. Wang, Z. Li, M. Hu, J. Tan, P. Hou, F. Li, X. Wang. Anisotropic electronic conduction in stacked two-dimensional titanium carbide. Sci. Rep., 2015, 5(1): 16329

[22] M. J. Hamer, J. Zultak, A. V. Tyurnina, V. Zólyomi, D. Terry, A. Barinov, A. Garner, J. Donoghue, A. P. Rooney, V. Kandyba, A. Giampietri, A. Graham, N. Teutsch, X. Xia, M. Koperski, S. J. Haigh, V. I. Fal’ko, R. V. Gorbachev, N. R. Wilson. Indirect to direct gap crossover in two-dimensional InSe revealed by angle-resolved photoemission spectroscopy. ACS Nano, 2019, 13(2): 2136

[23] T. V. Shubina, W. Desrat, M. Moret, A. Tiberj, O. Briot, V. Yu. Davydov, A. V. Platonov, M. A. Semina, B. Gil. InSe as a case between 3D and 2D layered crystals for excitons. Nat. Commun., 2019, 10(1): 3479

[24] T. J. Williams, A. A. Aczel, M. D. Lumsden, S. E. Nagler, M. B. Stone, J. Q. Yan, D. Mandrus. Magnetic correlations in the quasi-two-dimensional semiconducting ferromagnet CrSiTe₃. Phys. Rev. B, 2015, 92(14): 144404

[25] A. Rai, V. K. Sangwan, J. T. Gish, M. C. Hersam, D. G. Cahill. Anisotropic thermal conductivity of layered indium selenide. Appl. Phys. Lett., 2021, 118(7): 073101

[26] P. Jiang, X. Qian, R. Yang, L. Lindsay. Anisotropic thermal transport in bulk hexagonal boron nitride. Phys. Rev. Mater., 2018, 2(6): 064005

[27] P. Jiang, X. Qian, X. Gu, R. Yang, Probing anisotropic thermal conductivity of transition metal dichalcogenides MX₂ (M = Mo, W = S. Se) using time‐domain thermoreflectance. Adv. Mater., 2017, 29(36): 1701068

[28] N. Yang, F. Pei, J. Dou, Y. Zhao, Z. Huang, Y. Ma, S. Ma, C. Wang, X. Zhang, H. Wang, C. Zhu, Y. Bai, H. Zhou, T. Song, Y. Chen, Q. Chen. Improving heat transfer enables durable perovskite solar cells. Adv. Energy Mater., 2022, 12(24): 2200869

[29] S. Chandra, P. Dutta, K. Biswas. High-performance thermoelectrics based on solution-grown SnSe nanostructures. ACS Nano, 2022, 16(1): 7

[30] J. D. G. Greener, E. de Lima Savi, A. V. Akimov, S. Raetz, Z. Kudrynskyi, Z. D. Kovalyuk, N. Chigarev, A. Kent, A. Patané, V. Gusev. High-frequency elastic coupling at the interface of van der Waals nanolayers imaged by picosecond ultrasonics. ACS Nano, 2019, 13(10): 11530

[31] J. Peng, X. He, C. Shi, J. Leng, F. Lin, F. Liu, H. Zhang, W. Shi. Investigation of graphene supported terahertz elliptical metamaterials. Physica E, 2020, 124: 114309

[32] J. Leng, J. Peng, A. Jin, D. Cao, D. Liu, X. He, F. Lin, F. Liu. Investigation of terahertz high Q-factor of all-dielectric metamaterials. Opt. Laser Technol., 2022, 146: 107570

[33] X.HeF.Lin F.LiuW. Shi, 3D Dirac semimetals supported tunable terahertz BIC metamaterials, Nanophotonics 11(21), 4705 (2022)

[34] X. He, W. Cao. Tunable terahertz hybrid metamaterials supported by 3D Dirac semimetals. Opt. Mater. Express, 2023, 13(2): 413

[35] J. Yang, C. Liu, H. Xie, W. Yu. Anisotropic heat transfer properties of two-dimensional materials. Nanotechnology, 2021, 32(16): 162001

[36] S. E. Kim, F. Mujid, A. Rai, F. Eriksson, J. Suh, P. Poddar, A. Ray, C. Park, E. Fransson, Y. Zhong, D. A. Muller, P. Erhart, D. G. Cahill, J. Park. Extremely anisotropic van der Waals thermal conductors. Nature, 2021, 597(7878): 660

[37] S. Chen, A. Sood, E. Pop, K. E. Goodson, D. Donadio. Strongly tunable anisotropic thermal transport in MoS₂ by strain and lithium intercalation: First-principles calculations. 2D Mater., 2019, 6(2): 025033

[38] P. E. Blöchl. Projector augmented-wave method. Phys. Rev. B, 1994, 50(24): 17953

[39] G. Kresse, J. Furthmüller. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. Comput. Mater. Sci., 1996, 6(1): 15

[40] J. P. Perdew, K. Burke, M. Ernzerhof. Generalized gradient approximation made simple. Phys. Rev. Lett., 1996, 77(18): 3865

[41] J. Klimeš, D. R. Bowler, A. Michaelides. Chemical accuracy for the van der Waals density functional. J. Phys. :Condens. Matter, 2010, 22(2): 022201

[42] M. Omini, A. Sparavigna. Beyond the isotropic-model approximation in the theory of thermal conductivity. Phys. Rev. B, 1996, 53(14): 9064

[43] D. A. Broido, M. Malorny, G. Birner, N. Mingo, D. A. Stewart. Intrinsic lattice thermal conductivity of semiconductors from first principles. Appl. Phys. Lett., 2007, 91(23): 231922

[44] A.WardD. A. BroidoD.A. StewartG.Deinzer, Ab-initio theory of the lattice thermal conductivity in diamond, Phys. Rev. B 80(12), 125203 (2009)

[45] W. Li, L. Lindsay, D. A. Broido, D. A. Stewart, N. Mingo. Thermal conductivity of bulk and nanowire Mg₂Si_xSn_1-x alloys from first principles. Phys. Rev. B, 2012, 86(17): 174307

[46] W. Li, N. Mingo. Ultralow lattice thermal conductivity of the fully filled skutterudite YbFe₄Sb₁₂ due to the flat avoided-crossing filler modes. Phys. Rev. B, 2015, 91(14): 144304

[47] A. Togo, I. Tanaka. First principles phonon calculations in materials science. Scr. Mater., 2015, 108: 1

[48] P.GiannozziS. de GironcoliP.PavoneS.Baroni, Ab-initio calculation of phonon dispersions in semiconductors, Phys. Rev. B 43(9), 7231 (1991)

[49] W. Li, J. Carrete, N. A. Katcho, N. Mingo. ShengBTE: A solver of the Boltzmann transport equation for phonons. Comput. Phys. Commun., 2014, 185(6): 1747

[50] B. Niu, L. Zhong, W. Hao, Z. Yang, X. Duan, D. Cai, P. He, D. Jia, S. Li, Y. Zhou. First-principles study of the anisotropic thermal expansion and thermal transport properties in h-BN. Sci. China Mater., 2021, 64(4): 953

[51] W. Wan, S. Zhao, Y. Ge, Y. Liu. Phonon and electron transport in Janus monolayers based on InSe. J. Phys. :Condens. Matter, 2019, 31(43): 435501

[52] Y. Sun, S. Luo, X. G. Zhao, K. Biswas, S. L. Li, L. Zhang. InSe: A two-dimensional material with strong interlayer coupling. Nanoscale, 2018, 10(17): 7991

[53] Y. Sun, Y. Li, T. Li, K. Biswas, A. Patanè, L. Zhang. New polymorphs of 2D indium selenide with enhanced electronic properties. Adv. Funct. Mater., 2020, 30(31): 2001920

[54] J. Q. Hu, X. H. Shi, S. Q. Wu, K. M. Ho, Z. Z. Zhu. Dependence of electronic and optical properties of MoS₂ multilayers on the interlayer coupling and van Hove singularity. Nanoscale Res. Lett., 2019, 14(1): 288

[55] J. Liu, G. M. Choi, D. G. Cahill. Measurement of the anisotropic thermal conductivity of molybdenum disulfide by the time-resolved magneto-optic Kerr effect. J. Appl. Phys., 2014, 116(23): 233107

[56] W. Li, N. Mingo. Thermal conductivity of fully filled skutterudites: Role of the filler. Phys. Rev. B, 2014, 89(18): 184304

[57] S. K. S. S. Karthikeyan, P. Santhoshkumar, Y. C. Joe, S. H. Kang, Y. N. Jo, H. S. Kang, C. W. Lee. Understanding of the elastic constants, energetics, and bonding in dicalcium silicate using first-principles calculations. J. Phys. Chem. C, 2018, 122(42): 24235

[58] K. Yuan, X. Zhang, Z. Chang, Z. Yang, D. Tang. Pressure-induced anisotropic to isotropic thermal transport and promising thermoelectric performance in layered InSe. ACS Appl. Energy Mater., 2022, 5(9): 10690

[59] W. Ren, Y. Ouyang, P. Jiang, C. Yu, J. He, J. Chen. The impact of interlayer rotation on thermal transport across graphene/hexagonal boron nitride van der Waals heterostructure. Nano Lett., 2021, 21(6): 2634

[60] A. Y. Nobakht, S. Shin. Anisotropic control of thermal transport in graphene/Si heterostructures. J. Appl. Phys., 2016, 120(22): 225111

[61] Y. Ni, Y. A. Kosevich, S. Xiong, Y. Chalopin, S. Volz. Substrate-induced cross-plane thermal propagative modes in few-layer graphene. Phys. Rev. B, 2014, 89(20): 205413

RIGHTS & PERMISSIONS

Higher Education Press