Coherent control of thermal phonon transport in van der Waals superlattices

Abstract

van der Waals (vdW) heterostructures are a central focus of materials science and condensed matter physics due to the novel physical phenomena and properties obtained by precisely stacking heterogeneous atomically thin layers. vdW heterostructures are expected to allow for the coherent manipulation of THz lattice vibrations and hence heat conduction due to the ability to precisely control chemical composition at the atomic scale, but little work has focused on thermal transport in these materials. Here, we report an ab initio study of thermal transport in vdW superlattices consisting of alternating transition metal dichalcogenide atomic layers. Our calculations show that the lattice vibrational spectrum and scattering rates can be precisely manipulated by the choice of each atomically thin layer, resulting in materials with novel properties such as large thermal anisotropies approaching 200 and ultralow cross-plane thermal conductivities comparable to those of amorphous materials. Our work demonstrates how coherent manipulation of phonons in vdW superlattices can expand the property space beyond that occupied by natural materials and suggests an experimental route to realize these properties.

---

I. Introduction

van der Waals (vdW) heterostructures (HTs) consisting of stacked atomically thin layers in a precise sequence are of intense interest due to the novel physical phenomena and properties that can be realized.^1,2 Owing to the weak vdW interactions between the layers, considerable freedom is allowed to integrate highly diverse layers of metals, semiconductors, superconductors, or insulators. Recent pioneering efforts have enabled many novel devices, including field-effect transistors,^3,4 tunneling transistors,^5,6 photodetectors,^7,8 solar cells,⁹ nonvolatile memory devices^10,11 and light-emitting diodes.¹²

The close proximity of the atomic layers in vdW HTs leads to the coupling of charge carriers, excitons, photons and phonons between adjacent layers that results in fascinating physics. Extensive efforts have been devoted to modifying the electronic spectrum using vdW HTs. For instance, manipulating the coupling of charge carriers by inserting hexagonal boron nitride (hBN) layers between WSe₂/MoS₂ vdW gaps enables electronic band engineering.¹³ vdW HTs consisting of a mono/bilayer graphene coupled to a hBN substrate create long-wavelength moiré superlattices (SLs), leading to the formation of anomalous quantum Hall states and the Hofstadter butterfly at high magnetic fields.^14–16

vdW HTs also offer the possibility to control the atomic vibrational spectrum and hence thermal transport. During the past few decades, classical size effects of phonons in nanostructures have been extensively studied.^17–20 The possibility of coherently manipulating the phonon modes themselves has recently been of intense interest.^21–30 Because thermally occupied phonon modes that transport heat at room temperature possess wavelengths on the scale of the lattice constant, coherent control of phonons requires atomic scale periodicity and precision of composition. Such precision has been reported in oxide superlattices but the observed contrast in thermal properties was on the order of tens of percent due to the small acoustic contrast.²² Prior simulations predict coherent modifications of phonon dispersions in short-period Si/Ge and monolayer in-plane SLs as the periods of the SL approach the lattice constant.^26,28,29 Recently, similar two-dimensional lateral superlattices have been fabricated.^31,32 However, achieving coherent phonon transport using these structures experimentally remains challenging due to structural stability considerations, lattice matching, or smearing of the interface by atomic diffusion.³³

In contrast, vdW HTs are an experimentally realizable platform on which to manipulate the atomic vibrational spectrum at the length scale of single atoms. A few recent studies provide individual perspectives on the in-plane thermal conductivity and thermal conductance across interfaces of vdW HTs.^34–37 However, little work has examined how THz phonon frequencies and lifetimes can be coherently manipulated in vdW SLs to alter macroscopic thermal properties. Considering the rich family of two-dimensional materials and stacking sequences, the degree of freedom in designing artificial thermal properties using vdW SLs is large.

Here, we systematically investigate thermal transport in vdW SLs consisting of alternating transition metal dichalcogenide (TMDC) atomic layers using ab initio calculations. We show that precise manipulation of the lattice vibrational spectrum and scattering rates can be achieved by the choice of each atomically thin layer, yielding novel properties such as large thermal anisotropies approaching 200 and ultralow cross-plane thermal conductivities comparable to those of amorphous materials. Our work demonstrates how vdW SLs can be engineered to possess novel thermal properties that can be realized experimentally.

II. Computational methods

We calculated the vibrational and thermal properties of vdW SLs from first-principles combined with the Boltzmann transport equation (BTE). To calculate the lattice thermal conductivity, phonon frequencies and velocities are calculated by diagonalizing the dynamic matrix determined from harmonic interatomic force constants (IFCs) while lifetimes include the contribution of intrinsic phonon–phonon and isotope scatterings. Instead of using the virtual crystal approximation that is typically applied for SLs,^26,27 we derived harmonic and cubic anharmonic IFCs of the SL directly from density functional theory (DFT) based on the unit cell spanning one period of the SL, which therefore includes the coherent interference of lattice waves. As a well-established method, good agreement with measured results has been reported for many materials, such as Si,^38,39 diamond,⁴⁰ GaN,⁴¹ CoSb₃ ^42,43 and PbTe.⁴⁴ More technique details can be found in previous publications.^38,40,41,43 In this study, we determine the IFCs from density functional theory calculations using the projector augmented wave method, as implemented in the Vienna ab initio simulation package (VASP).⁴⁵ For all materials, both relaxation and static energy calculations use the van der Waals density functional optB88-vdW and the Perdew–Burke–Ernzerhof (PBE) generalized gradient approximation. The unit cells were fully relaxed using a cutoff energy of 400 eV and 12 × 12 × 3 k-point sampling grid. The obtained in-plane lattice constants are 3.17, 3.17 and 3.30 Å for MoS₂, WS₂ and MoSe₂, respectively, which therefore produce perfect lattice matching in the MoS₂/WS₂ SL but 2% tensile and compressive strain in the WS₂ and MoSe₂ layers of the MoSe₂/WS₂ SL, respectively. 3 × 3 × 2 and 3 × 3 × 1 supercells were constructed to calculate the IFCs of structures with bi- and multi-layers in the unit cell, respectively. An interaction cutoff distance of 6.3 Å was selected for the calculation of anharmonic IFCs to ensure the convergence of κ_L. To integrate the κ_L over Brillouin zone using the tetrahedron method, we used a 27 × 27 × 7 q-mesh for structures with bilayers in the unit cell and appropriately adapted ones for those with multilayers.

III. Results and discussion

The SLs we consider consist of MX₂/M′X′₂ (M/M′ = transition metals, X/X′ = chalcogens) with two periodically alternating TMDC single-layers. Two periods of this SL structure are shown in Fig. 1a. Very recently, a similar structure has been experimentally synthesized on a wafer scale.⁴⁶ The primitive unit cell, shown in Fig. 1b possesses a bilayer structure with three atoms in each layer, characterized by strong intralayer bonding and weak interlayer vdW interactions, similar to those of naturally occurring TMDCs. Therefore, they share the hexagonal Brillouin zone shown in Fig. 1c. In our study we constrain the two materials to have the same in-plane lattice constants, leading to negligible strain for the MoS₂/WS₂ SL but 2% tensile and compressive in-plane strains in the WS₂ and MoSe₂ layers for the MoSe₂/WS₂ SL, respectively. Our calculations indicate that strain will not affect our conclusions since the cross-plane thermal conductivities of this SL are little affected while the in-plane ones are only decreased by ∼9% (see ESI Note 1† for details).

Fig. 1 Crystal structures and phonon dispersions. (a) Schematic illustration of the AB-type vdW SL, showing two periods of the SL, (b) primitive unit cell and (c) Brillouin zone of the TMDC SLs MX₂/M′X′₂ (M/M′ = transition metals, X/X′ = chalcogens). Phonon dispersions of bulk (d) MoS₂, (e) WS₂, (f) MoS₂/WS₂ and (g) MoSe₂/WS₂ SLs. Calculations (solid lines) and experiments (red circles).⁵⁰ The lowest six branches in the MoS₂/WS₂ SL exhibit distinct band splitting along the in-plane directions and band gaps (indicated by the red arrow) at the zone boundary of the cross-plane direction (Γ–A), both of which occur to a lesser degree in the MoSe₂/WS₂ SL.

We first validate the accuracy of our calculations for bulk MoS₂, WS₂ and MoSe₂. The fully relaxed lattice constants obtained from the vdW-DF functional deviate from experimental values by less than 1%.^47–49Fig. 1d and e show the phonon dispersions of bulk MoS₂ and WS₂ (that of MoSe₂ is given in ESI Fig. S2†), exhibiting good agreement with the data from inelastic neutron scattering.⁵⁰ The in-plane longitudinal sound velocity is calculated to be 6730 m s⁻¹ for bulk MoS₂, which is close to our experimental value 7000 ± 40 m s⁻¹ measured using ultrafast transient grating spectroscopy⁵¹ and previous ab initio predictions of 6500–6700 m s⁻¹.⁵² All materials feature six dispersive phonon branches at the low-frequency range, which are separated from the relatively flat optical modes by a frequency gap. The gaps differ among the materials, with the largest one equaling 2.9 THz in WS₂, MoSe₂ lacking a gap, and that of MoS₂ in between at 1.0 THz. The value of the gap is primarily due to the mass difference between the heavier transition metals and lighter chalcogens (Atomic mass (amu): Mo-95.94, W-183.84, S-32.065, Se-78.96), which dominate the vibrations of lower and higher branches, respectively (see ESI Fig. S3† for the density of states (DOS)). Another notable feature is the much lower cross-plane group velocities compared to the in-plane values, as indicated by the phonon dispersions along the Γ–A direction. This observation is expected considering the weak vdW interaction between layers.

The phonon dispersions of the MoS₂/WS₂ and MoSe₂/WS₂ SLs are presented in Fig. 1f and g, respectively. While they inherit many characteristics of those in bulk TMDCs, an obvious band splitting arises in both SLs due to the mass difference. For example, the large frequency gap featured in WS₂ is occupied by three branches contributed by the MoSe₂ layer in the MoSe₂/WS₂ SL. In particular, the lowest six branches in the MoS₂/WS₂ SL exhibit distinct band splitting along the in-plane directions and band gaps (indicated by the red arrow) at the zone boundary of the cross-plane direction (Γ–A), both of which occur to a lesser degree in the MoSe₂/WS₂ SL.

We then calculate the thermal conductivity κ_Lversus temperature of bulk MoS₂, MoSe₂ and WS₂ with the naturally occurring isotopic composition, as shown in Fig. 2a. The calculated values at 300 K are 102.6, 40.9 and 153.5 W m⁻¹ K⁻¹ for the in-plane thermal conductivity κ_In and 5.1, 2.6 and 3.9 W m⁻¹ K⁻¹ for the cross-plane κ_Cross, respectively. The experimental data are scattered for MoS₂, falling in the range of 82–112 W m⁻¹ K⁻¹ and 2.0–4.8 W m⁻¹ K⁻¹ along the in- and cross-plane directions, respectively.^53,54 A recent time-domain thermoreflectance measurement reported κ_In = 35 W m⁻¹ K⁻¹ and κ_Cross = 2.6 W m⁻¹ K⁻¹ for MoSe₂.⁵⁴ For WS₂, the measured in-plane data are close to 124 W m⁻¹ K⁻¹ while the cross-plane values vary from 1.7 to 2.9 W m⁻¹ K⁻¹.^54,55 The decreasing orders WS₂ > MoS₂ > MoSe₂ for κ_In and MoS₂ > WS₂ > MoSe₂ for κ_Cross from our calculations are consistent with measurements⁵⁴ and prior ab initio calculations.⁵² We find that the impurity scattering more strongly suppresses the κ_In due to the significant contribution of higher frequency modes along this crystal plane, possibly explaining why our results agree better with the κ_Cross reported by Jiang et al.⁵⁴ (see ESI Note 2† for details).

Fig. 2 κ _L and its spectral contribution of MoS₂, WS₂, MoSe₂, MoS₂/WS₂ and MoSe₂/WS₂ crystals. (a) Cross-plane (bottom) and in-plane (top) κ_Lversus temperature. (b) Thermal anisotropy, κ_In/κ_Cross, versus temperature. (c) Cross-plane and (d) in-plane thermal conductivity contribution κ_fversus frequency. The thermal anisotropy ratio of the MoS₂/WS₂ SL is enhanced by a factor of ∼2 at different temperatures due to the mainly preserved κ_In and substantially reduced κ_Cross while that of the MoSe₂/WS₂ SL is decreased by ∼20%.

We next calculate the thermal conductivity of the MoS₂/WS₂ and MoSe₂/WS₂ SLs. Since MoS₂ has much higher κ_In and κ_Cross than those of MoSe₂, classical effective medium theories would predict that a composite with MoS₂ would have higher κ_L than one with MoSe₂. However, our calculation shows that the MoS₂/WS₂ SL has much lower κ_Cross and larger thermal anisotropy than those of the MoSe₂/WS₂ SL. As shown in Fig. 2a, the MoS₂/WS₂ SL exhibits the lowest κ_Cross over the entire temperature range, which is substantially reduced compared to the corresponding values of either constituent material. However, the κ_Cross of the MoSe₂/WS₂ SL is very close to that of MoSe₂, especially at low temperatures. For example, with respect to the lower value of the corresponding bulk materials, the κ_Cross of the MoS₂/WS₂ and MoSe₂/WS₂ SLs is reduced by 62.4% and 6.7% at 300 K, respectively.

In contrast to the reduction of the κ_Cross, the κ_In of both SLs lies in between those of corresponding bulk materials. We also note that the κ_In of the MoS₂/WS₂ SL is closer to the average value of the corresponding bulk constituents, lower by 7.7% at 300 K, whereas that of the MoSe₂/WS₂ SL has a value 32.1% lower. Due to the different variation of κ_In and κ_Cross from those of the bulk materials, both SLs exhibit altered thermal anisotropy ratios κ_In/κ_Cross, as illustrated in Fig. 2b. In the present work, we evaluate the variation of thermal anisotropy in SLs using the constituent material with larger thermal anisotropy ratio as the reference (WS₂ for both SLs). The thermal anisotropy ratio of the MoS₂/WS₂ SL is enhanced by a factor of ∼2 while that of the MoSe₂/WS₂ SL is decreased by ∼20% throughout our temperature range.

To understand the unusual thermal anisotropy in these SLs, we first plot the thermal conductivity distribution κ_fversus frequency at 300 K in Fig. 2c and d. We observe that the cross-plane κ_f is mainly contributed by phonons within a narrow frequency range (f < ∼2 THz) for each material. For both SLs, the reduction of the cross-plane κ_f contribution with respect to those of the bulk counterparts decreases as f approaches zero. Above this low-frequency range, the κ_f of the MoS₂/WS₂ SL is substantially suppressed and becomes negligible above 2.1 THz while that of the MoSe₂/WS₂ SL is very close to the value of MoSe₂ above 1.2 THz.

Along the in-plane direction, the major κ_f contribution comes from a broader frequency range (up to ∼5 THz) for all materials, as shown in Fig. 2d. Similarly, the in-plane κ_f of the two SLs almost overlaps with that of corresponding bulk crystals below ∼1 THz. Above this range, the κ_f of both SLs lies in between that of constituent materials although the κ_f of the MoSe₂/WS₂ SL is significantly lower than the average value of its constituent materials. Furthermore, the lowest six branches contribute more than 94% of the total κ_In and κ_Cross in each material.

These calculations indicate that the origin of the increased thermal anisotropy in the MoS₂/WS₂ SL lies in the mainly preserved and decreased contribution of modes with frequency less than around 7 THz along the in- and cross plane directions, respectively. This difference among directions must be due to changes in either harmonic properties such as group velocities and phonon scattering phase space or the strength of anharmonic coupling determined by the magnitude of the cubic force constants.

We first probe the role of scattering phase space and anharmonicity by examining the lifetimes of the MoS₂/WS₂ SL. We find that the lifetimes of the lowest six branches can be divided into two sets of three, for which each have lifetimes similar to either those of bulk MoS₂ and WS₂. In Fig. 3a and b, we show the lifetimes of the first transverse acoustic branch B1 (from the WS₂ layer) and the highest low-lying branch B6 (from the MoS₂ layer), as typical examples. This splitting of the lifetimes between each constituent material arises from the splitting of previously degenerate phonon modes, as shown in Fig. 3c. Due to the bilayer structure of the unit cell, the vibrations of the MoS₂ and WS₂ layers separately contribute to three branches, which, importantly, almost overlap with those in bulk MoS₂ and WS₂. Therefore, the κ_In of the SLs approaches the average value of bulk constituent materials considering that the individual κ_L contributions of the six branches are comparable.

Fig. 3 Phonon lifetimes and dispersions in MoS₂, WS₂ and MoS₂/WS₂ crystals. Lifetimes of (a) the first transverse acoustic (B1) and (b) highest low-lying (B6) branches. (c) Phonon dispersions of the six low-lying branches B1–B6 along the in-plane directions in the MoS₂/WS₂ SL, along with three separate branches from bulk WS₂ and MoS₂, respectively. The splitting of the lifetimes between each constituent material arises from the splitting of phonon dispersions. (d) Lifetimes of B1–B6 along the Γ–A direction. The relatively suppressed lifetimes of the MoS₂/WS₂ SL near the zone edge have only a minor contribution to the reduction of the total κ_Cross.

The splitting of phonon dispersions changes phonon scattering channels and thus lifetimes, which depend on the frequency range of acoustic (including the six low-lying branches as they all have acoustic-like behavior) and optical branches. As shown in Fig. 3a and b, some lifetimes of B1 are relatively smaller than those in bulk WS₂ while most lifetimes of B6 are very close to the corresponding values in bulk MoS₂. The former behavior is mainly due to the enhanced aoo scattering, or scattering between one acoustic and two optical phonons, which is the dominant process for these modes of B1. For example, 85.3% of the scattering rate for the mode M1 (q = [0, 8/27, 0], 2.9 THz) from B1 is contributed by aoo scattering in the MoS₂/WS₂ SL, which is larger by 14.5% with respect to that in WS₂. The enhancement can be explained by the widening of optical spectrum, from 4.3 THz in WS₂ to 5.4 THz in the MoS₂/WS₂ SL, resulting in more aoo scattering channels. In contrast, the widths of acoustic and optical spectra in this SL are very close to those in MoS₂, resulting in similar scattering rates of B6. For the MoSe₂/WS₂ SL, the width of optical spectrum further increases to 7.6 THz, which thus causes stronger enhancement of aoo scattering for phonon modes contributed by WS₂ layers. Meanwhile, the appearance of three branches in the band gap allows more aao channels. Consequently, the κ_In of the MoSe₂/WS₂ SL is further lowered from the average value of corresponding bulk ones.

The large thermal anisotropy of the SLs is due to mainly preserved κ_In, the origin of which is discussed above, and due to reduced κ_Cross. We now examine if the lifetimes of phonons along the Γ–A direction can explain this reduction in κ_Cross. Fig. 3d shows that, as the q vector approaches the zone edge, the lifetimes of the MoS₂/WS₂ SL become relatively smaller than those in bulk MoS₂ and WS₂. However, recalculating κ_Cross of WS₂ using these lifetimes leads to only a decrease of 0.23 W m⁻¹ K⁻¹, accounting for 9.6% of the overall reduction (from 3.9 in WS₂ to 1.5 W m⁻¹ K⁻¹ in the MoS₂/WS₂ SL). Therefore, the major reduction of the κ_Cross must result from alterations to the phonon dispersion.

Fig. 4a and b show the z-component group velocities normalized by the magnitude of the group velocity of the second transverse acoustic (B2) and the B6 branches in the WS₂, MoS₂/WS₂ and MoSe₂/WS₂ crystals. The corresponding values in MoS₂ and MoSe₂ are comparable to or higher than those in WS₂. Clearly, the MoS₂/WS₂ SL has substantially smaller normalized cross-plane group velocities compared to WS₂ except for the low-frequency range, where the group velocity vectors of a few modes lie directly along the cross-plane direction. Above ∼2 THz, the z-component of the group velocity is only less than 3% of the magnitude of overall group velocities for most phonons in the MoS₂/WS₂ SL, and those with larger z-components at the high-frequency range also have small absolute group velocities, explaining the corresponding negligible κ_L contribution. Conversely, there is no significant suppression in the normalized z-component group velocities of the MoSe₂/WS₂ SL.

Fig. 4 Normalized cross-plane group velocities and isoenergy surfaces. Cross-plane group velocities normalized by the magnitude of group velocities of (a) the second transverse acoustic (B2) and (b) the B6 branches for the WS₂, MoS₂/WS₂ and MoSe₂/WS₂ crystals. Isoenergy surfaces of B2 in the A–Γ–M plane for (c) MoS₂/WS₂ and (d) MoSe₂/WS₂ SLs. The substantially smaller normalized cross-plane group velocities of the MoS₂/WS₂ SL arise from phonon focusing effects, as indicated by the smaller curvatures of the isoenergy surfaces.

We further reveal the phonon focusing effect by visualizing the isoenergy surfaces of B2 in the A–Γ–M plane for both SLs, as plotted in Fig. 4c and d. The comparison clearly indicates smaller curvatures of the isoenergy surfaces in the MoS₂/WS₂ SL. In the middle range (2–3 THz), the isoenergy surfaces of the MoS₂/WS₂ SL are almost completely flat, resulting in phonon focusing along the in-plane direction. Above 3 THz, the MoSe₂/WS₂ SL exhibits dispersive isoenergy surfaces, in contrast to the relatively flat ones in the MoS₂/WS₂ SL. The flattening of the isoenergy surfaces thus leads to a decrease in the slope of the surface and hence the cross-plane group velocities, increasing the thermal anisotropy. A similar effect has also been reported in graphite by tuning the in-plane/cross-plane bonding strength, resulting in a negative correlation between the bonding strength and κ_L along the orthogonal direction.⁵⁶

We now seek to clarify the physical mechanisms for the flattening of the isoenergy surfaces. Early works suggested that phonon group velocity vectors are preferentially focused along the crystalline directions with larger elastic constants.⁵⁷ However, elastic anisotropy cannot explain the enhanced thermal anisotropy in the MoS₂/WS₂ SL, which has a ratio of elastic constants C₁₁/C₃₃ = 4.9, compared to 4.5 for MoS₂ and 5.1 for WS₂ from ab initio calculations.

Zone folding effects have been used to explain the decrease of cross-plane group velocities in SLs.^58–61 The present SLs have approximately the same Brillouin zone size as bulk TMDCs, suggesting that zone-folding does not occur. However, if we observe that approximating the three atoms in each layer of the unit cell as an equivalent atom simplifies bulk TMDCs to parallel monatomic chains, while those of the SLs are diatomic chains, the Brillouin zone of the SLs must be folded once along the cross-plane direction. As a result, flattened bands and band gaps are expected to appear near the zone edges where the Bragg reflection condition is satisfied.⁵⁹

These features are indeed observed in the dispersion of the MoS₂/WS₂ SL along the Γ–A direction, as shown in Fig. 1f. Due to the Bragg reflection, the group velocity vanishes near the zone edge. In addition, anticrossing effects occur inside the zone due to Bragg reflections⁶⁰ (ESI Fig. S4†), further suppressing the cross-plane group velocities. The band gaps at zone edges and anticrossings yield sharp peaks in the DOS projected along the z axis for the SL, a unique feature distinct from bulk TMDCs (ESI Fig. S5†). Furthermore, the optical phonons of the MoS₂ and WS₂ layers have no common frequencies and thus become localized, resulting in flat dispersions. Due to the negligible κ_L contribution of these branches, this effect has little influence on the thermal anisotropy of TMDC SLs. However, it may cause significant enhancement of thermal anisotropy of vdW materials for which these modes contribute, such as for SnSe and SnS.⁶²

In contrast, the MoSe₂/WS₂ SL shows no clear band gaps along the Γ–A direction, as shown in Fig. 1g. The major difference between the two SLs is the mass ratio of neighboring layers, being 1.55 for the MoS₂/WS₂ SL compared to 1.02 for the MoSe₂/WS₂ SL. As in the simple diatomic chain model, the frequency gap vanishes as the masses become equal. We validate this explanation by computing how phonon dispersions change as the mass of one of the bilayers is varied while keeping the same force constants as in MoS₂. Modifying the mass of Mo or S significantly alters the magnitude of the band gaps at the zone edge (ESI Fig. S6†).

A key requirement for coherent manipulation of phonons is that the coherence length must be much longer than the imposed periodicity. Our calculation shows that for the MoS₂/WS₂ SL, phonons with MFPs longer than 200 nm contribute 88% to the κ_Cross, indicating that most heat is indeed carried by phonons with coherence lengths far exceeding SL periodicity of 12.42 Å (ESI Fig. S7†). As the period increases, scattering events will cause decoherence, preventing interference and thereby leading to the effective medium limit. For the MoSe₂/WS₂ SL, although phonons with MFPs longer than 200 nm contribute 82% to the κ_Cross, coherent interference has little effect on the thermal transport due to the negligible mass difference between layers.

The thermal anisotropy in TMDC SLs can be further enhanced by additional zone folding. For example, we construct two representative SLs with four atomic layers in each period: WS₂/2MoS₂/WSe₂ (mass ratio = 1 : 0.65 : 0.65 : 1.38) and WS₂/2MoSe₂/WSe₂ (mass ratio = 1 : 1.02 : 1.02 : 1.38), producing a thermal anisotropy of 176 and 50 at 300 K, respectively, as shown in Fig. 5. The large thermal anisotropy of the former SL is comparable to that (100–300) of high-quality graphite.⁶³ Meanwhile, its κ_Cross decreases to 0.4 W m⁻¹ K⁻¹, which is close to that of amorphous polymers (∼0.35 W m⁻¹ K⁻¹ for bulk polyethylene).⁶³ These results agree well with our expectation considering the dominant role of band splitting and phonon focusing (ESI Fig. S8†).

Fig. 5 Coherent manipulation of thermal transport for SLs with four atomic layers in each period (WS₂/2MoS₂/WSe₂ with mass ratio = 1 : 0.65 : 0.65 : 1.38 and WS₂/2MoSe₂/WSe₂ with mass ratio = 1 : 1.02 : 1.02 : 1.38). (a) Cross-plane (bottom) and in-plane (top) κ_L with respect to temperature. (b) Thermal anisotropy, κ_In/κ_Cross, versus temperature. The thermal anisotropies of the WS₂/2MoS₂/WSe₂ SL are substantially further enhanced by additional zone folding.

In this work, the very similar lattice structures of different TMDCs result in band splitting and hence mainly preserved κ_In, leading to enhanced thermal anisotropy. Constructing vdW HTs as a versatile strategy to enhance thermal anisotropy will benefit designing efficient thermal management of electronics and exploring extreme values of thermal conductivity such as low values for thermoelectrics. We expect that for vdW SLs consisting of materials with dissimilar crystal structures, a more dramatic modification of phonon scattering phase space may have stronger influence on thermal conductivity along different crystallographic directions. For instance, an unusual enhancement of thermal conductivity of the short-period Si/Ge SL compared to those of both its constituents has been observed due to the suppressed scattering of acoustic phonons by optical phonons.²⁶ For a given vdW SL, twisting plane rotation angles can further manipulate the vibrational spectrum and scattering strength.⁶⁴

It should be noted that the modification of phonon dispersions in vdW SLs also depends on weak vdW interactions, which allow for the simplification of diatomic chain models. Such a treatment may fail for traditional short-period SLs coupled via strong covalent bonds, in which atomic interactions beyond nearest neighbors are not negligible. Therefore, the strength of vdW interactions, which can be tuned by for example applying pressure,⁶⁵ offers another degree of freedom for the design of thermal properties using vdW SLs.

IV. Conclusions

In summary, we report an ab initio study of coherent phonons and thermal transport in vdW SLs as the periodicity approaches the atomically thin limit. Our calculations show that these materials achieve novel thermal properties such as large thermal anisotropies approaching 200 and ultralow cross-plane thermal conductivities comparable to those of amorphous materials. Further coherent control of the thermal transport may be achieved using vdW SLs composed of atomic layers of diverse lattice structures, leading to strongly modified phonon scattering phase space and anharmonicity. Our study indicates a promising route to experimentally realize solids with thermal properties that do not occur naturally in atomically precise vdW superlattices.

Author contributions

R. G. performed the calculations. All authors contributed to the analysis and discussion of the results. R. G. and A. J. M. wrote the manuscript.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by the “GIST-Caltech Research Collaboration” Project through a grant provided by Gwangju Institute of Science and Technology in 2017 and the DARPA MATRIX program under Grant No. HR0011-15-2-0039.

References

1. A. K. Geim and I. V. Grigorieva, Nature, 2013, 499, 419–425.
2. Y. Liu, N. O. Weiss, X. Duan, H.-C. Cheng, Y. Huang and X. Duan, Nat. Rev. Mater., 2016, 1, 16042.
3. W. J. Yu, Z. Li, H. Zhou, Y. Chen, Y. Wang, Y. Huang and X. Duan, Nat. Mater., 2013, 12, 246–252.
4. W. Yang, G. Chen, Z. Shi, C.-C. Liu, L. Zhang, G. Xie, M. Cheng, D. Wang, R. Yang and D. Shi, Nat. Mater., 2013, 12, 792–797.
5. T. Georgiou, R. Jalil, B. D. Belle, L. Britnell, R. V. Gorbachev, S. V. Morozov, Y.-J. Kim, A. Gholinia, S. J. Haigh and O. Makarovsky, Nat. Nanotechnol., 2013, 8, 100–103.
6. L. Britnell, R. V. Gorbachev, R. Jalil, B. D. Belle, F. Schedin, A. Mishchenko, T. Georgiou, M. I. Katsnelson, L. Eaves, S. V. Morozov, N. M. Peres, J. Leist, A. K. Geim, K. S. Novoselov and L. A. Ponomarenko, Science, 2012, 335, 947–950.
7. W. Zhang, C. P. Chuu, J. K. Huang, C. H. Chen, M. L. Tsai, Y. H. Chang, C. T. Liang, Y. Z. Chen, Y. L. Chueh, J. H. He, M. Y. Chou and L. J. Li, Sci. Rep., 2014, 4, 3826.
8. R.-J. Shiue, Y. Gao, Y. Wang, C. Peng, A. D. Robertson, D. K. Efetov, S. Assefa, F. H. Koppens, J. Hone and D. Englund, Nano Lett., 2015, 15, 7288–7293.
9. W. J. Yu, Y. Liu, H. Zhou, A. Yin, Z. Li, Y. Huang and X. Duan, Nat. Nanotechnol., 2013, 8, 952–958.
10. S. Bertolazzi, D. Krasnozhon and A. Kis, ACS Nano, 2013, 7, 3246–3252.
11. M. S. Choi, G.-H. Lee, Y.-J. Yu, D.-Y. Lee, S. H. Lee, P. Kim, J. Hone and W. J. Yoo, Nat. Commun., 2013, 4, 1624.
12. F. Withers, O. D. Pozo-Zamudio, A. Mishchenko, A. Rooney, A. Gholinia, K. Watanabe, T. Taniguchi, S. Haigh, A. Geim and A. Tartakovskii, Nat. Mater., 2015, 14, 301–306.
13. H. Fang, C. Battaglia, C. Carraro, S. Nemsak, B. Ozdol, J. S. Kang, H. A. Bechtel, S. B. Desai, F. Kronast, A. A. Unal, G. Conti, C. Conlon, G. K. Palsson, M. C. Martin, A. M. Minor, C. S. Fadley, E. Yablonovitch, R. Maboudian and A. Javey, Proc. Natl. Acad. Sci. U. S. A., 2014, 111, 6198–6202.
14. B. Hunt, J. D. Sanchez-Yamagishi, A. F. Young, M. Yankowitz, B. J. LeRoy, K. Watanabe, T. Taniguchi, P. Moon, M. Koshino, P. Jarillo-Herrero and R. C. Ashoori, Science, 2013, 340, 1427–1430.
15. C. Dean, L. Wang, P. Maher, C. Forsythe, F. Ghahari, Y. Gao, J. Katoch, M. Ishigami, P. Moon and M. Koshino, Nature, 2013, 497, 598–602.
16. L. Ponomarenko, R. Gorbachev, G. Yu, D. Elias, R. Jalil, A. Patel, A. Mishchenko, A. Mayorov, C. Woods and J. Wallbank, Nature, 2013, 497, 594–597.
17. G. Chen, Nanoscale energy transport and conversion, Oxford University Press, New York, 2005.
18. C. J. Vineis, A. Shakouri, A. Majumdar and M. G. Kanatzidis, Adv. Mater., 2010, 22, 3970–3980.
19. T. Luo and G. Chen, Phys. Chem. Chem. Phys., 2013, 15, 3389–3412.
20. D. G. Cahill, P. V. Braun, G. Chen, D. R. Clarke, S. Fan, K. E. Goodson, P. Keblinski, W. P. King, G. D. Mahan and A. Majumdar, Appl. Phys. Rev., 2014, 1, 011305.
21. M. N. Luckyanova, J. Garg, K. Esfarjani, A. Jandl, M. T. Bulsara, A. J. Schmidt, A. J. Minnich, S. Chen, M. S. Dresselhaus, Z. Ren, E. A. Fitzgerald and G. Chen, Science, 2012, 338, 936–939.
22. J. Ravichandran, A. K. Yadav, R. Cheaito, P. B. Rossen, A. Soukiassian, S. Suresha, J. C. Duda, B. M. Foley, C.-H. Lee and Y. Zhu, Nat. Mater., 2014, 13, 168–172.
23. J. Maire, R. Anufriev, R. Yanagisawa, A. Ramiere, S. Volz and M. Nomura, Sci. Adv., 2017, 3, e1700027.
24. M. Maldovan, Nature, 2013, 503, 209–217.
25. M. Maldovan, Nat. Mater., 2015, 14, 667–674.
26. J. Garg, N. Bonini and N. Marzari, Nano Lett., 2011, 11, 5135–5141.
27. J. Garg and G. Chen, Phys. Rev. B: Condens. Matter, 2013, 87, 140302.
28. X. Mu, T. Zhang, D. B. Go and T. Luo, Carbon, 2015, 83, 208–216.
29. C. da Silva, F. Saiz, D. A. Romero and C. H. Amon, Phys. Rev. B: Condens. Matter, 2016, 93, 125427.
30. B. Latour and Y. Chalopin, Phys. Rev. B: Condens. Matter, 2017, 95, 214310.
31. Z. Zhang, P. Chen, X. Duan, K. Zang, J. Luo and X. Duan, Science, 2017, 357, 788–792.
32. K. Chiu, K. Huang, C. Chen, Y. Lai, X. Zhang, E. Lin, M. Chuang, J. Wu and Y. Lee, Adv. Mater., 2018, 30, 1704796.
33. P. Chen, N. Katcho, J. Feser, W. Li, M. Glaser, O. Schmidt, D. G. Cahill, N. Mingo and A. Rastelli, Phys. Rev. Lett., 2013, 111, 115901.
34. B. Liu, J. A. Baimova, C. D. Reddy, A. W.-K. Law, S. V. Dmitriev, H. Wu and K. Zhou, ACS Appl. Mater. Interfaces, 2014, 6, 18180–18188.
35. B. Liu, F. Meng, C. D. Reddy, J. A. Baimova, N. Srikanth, S. V. Dmitriev and K. Zhou, RSC Adv., 2015, 5, 29193–29200.
36. Z. Ding, Q.-X. Pei, J.-W. Jiang, W. Huang and Y.-W. Zhang, Carbon, 2016, 96, 888–896.
37. Z. Yan, L. Chen, M. Yoon and S. Kumar, Nanoscale, 2016, 8, 4037–4046.
38. D. Broido, M. Malorny, G. Birner, N. Mingo and D. Stewart, Appl. Phys. Lett., 2007, 91, 231922.
39. X. Wang and B. Huang, Sci. Rep., 2014, 4, 6399.
40. A. Ward, D. Broido, D. A. Stewart and G. Deinzer, Phys. Rev. B: Condens. Matter, 2009, 80, 125203.
41. L. Lindsay, D. Broido and T. Reinecke, Phys. Rev. Lett., 2012, 109, 095901.
42. W. Li and N. Mingo, Phys. Rev. B: Condens. Matter, 2014, 90, 094302.
43. R. Guo, X. Wang and B. Huang, Sci. Rep., 2015, 5, 7806.
44. N. Shulumba, O. Hellman and A. J. Minnich, Phys. Rev. B: Condens. Matter, 2017, 95, 014302.
45. G. Kresse and J. Furthmuller, Phys. Rev. B: Condens. Matter, 1996, 54, 11169–11186.
46. K. Kang, K.-H. Lee, Y. Han, H. Gao, S. Xie, D. A. Muller and J. Park, Nature, 2017, 550, 23905.
47. V. Kalikhman, Izv. Akad. Nauk SSSR, Neorg. Mater., 1983, 19, 1060–1065.
48. K. Bronsema, J. D. Boer and F. Jellinek, Z. Anorg. Allg. Chem., 1986, 540, 15–17.
49. W. J. Schutte, J. L. D. Boer and F. Jellinek, J. Solid State Chem., 1987, 70, 207–209.
50. N. Wakabayashi, H. Smith and R. Nicklow, Phys. Rev. B: Solid State, 1975, 12, 659.
51. T. Kim, D. Ding, J.-H. Yim, Y.-D. Jho and A. J. Minnich, APL Mater., 2017, 5, 086105.
52. D. O. Lindroth and P. Erhart, Phys. Rev. B: Condens. Matter, 2016, 94, 115205.
53. J. Liu, G.-M. Choi and D. G. Cahill, J. Appl. Phys., 2014, 116, 233107.
54. P. Jiang, X. Qian, X. Gu and R. Yang, Adv. Mater., 2017, 29, 1701068.
55. A. Pisoni, J. Jacimovic, R. Gal, B. Nfrdi, H. Berger, Z. Rvay and L. Forr, Scr. Mater., 2016, 114, 48–50.
56. Z. Wei, Y. Chen and C. Dames, Appl. Phys. Lett., 2013, 102, 011901.
57. B. Taylor, H. Maris and C. Elbaum, Phys. Rev. Lett., 1969, 23, 416.
58. C. Colvard, T. Gant, M. Klein, R. Merlin, R. Fischer, H. Morkoc and A. Gossard, Phys. Rev. B: Condens. Matter, 1985, 31, 2080.
59. D. Hurley, S. Tamura, J. Wolfe and H. Morkoc, Phys. Rev. Lett., 1987, 58, 2446.
60. S. Tamura and J. Wolfe, Phys. Rev. B: Condens. Matter, 1987, 35, 2528.
61. M. Simkin and G. Mahan, Phys. Rev. Lett., 2000, 84, 927–930.
62. R. Guo, X. Wang, Y. Kuang and B. Huang, Phys. Rev. B: Condens. Matter, 2015, 92, 115202.
63. Y. S. Touloukian, Thermophysical properties of matter, IFI/Plenum, New York, 1970.
64. H. Li, H. Ying, X. Chen, D. L. Nika, A. I. Cocemasov, W. Cai, A. A. Balandin and S. Chen, Nanoscale, 2014, 6, 13402–13408.
65. X. Meng, T. Pandey, S. Fu, J. Yang, J. Jeong, K. Chen, A. Singh, F. He, X. Xu and A. K. Singh, 2017, arXiv preprint arXiv:1708.03849.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/C8NR02150C

---