Skip to main content

Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.

Extremely anisotropic van der Waals thermal conductors

Abstract

The densification of integrated circuits requires thermal management strategies and high thermal conductivity materials1,2,3. Recent innovations include the development of materials with thermal conduction anisotropy, which can remove hotspots along the fast-axis direction and provide thermal insulation along the slow axis4,5. However, most artificially engineered thermal conductors have anisotropy ratios much smaller than those seen in naturally anisotropic materials. Here we report extremely anisotropic thermal conductors based on large-area van der Waals thin films with random interlayer rotations, which produce a room-temperature thermal anisotropy ratio close to 900 in MoS2, one of the highest ever reported. This is enabled by the interlayer rotations that impede the through-plane thermal transport, while the long-range intralayer crystallinity maintains high in-plane thermal conductivity. We measure ultralow thermal conductivities in the through-plane direction for MoS2 (57 ± 3 mW m−1 K−1) and WS2 (41 ± 3 mW m−1 K−1) films, and we quantitatively explain these values using molecular dynamics simulations that reveal one-dimensional glass-like thermal transport. Conversely, the in-plane thermal conductivity in these MoS2 films is close to the single-crystal value. Covering nanofabricated gold electrodes with our anisotropic films prevents overheating of the electrodes and blocks heat from reaching the device surface. Our work establishes interlayer rotation in crystalline layered materials as a new degree of freedom for engineering-directed heat transport in solid-state systems.

Main

Anisotropic thermal conductors, in which heat flows faster in one direction than in another, can be characterized by the thermal conductivity anisotropy ratio ρ (= κf/κs) between the thermal conductivities along the fast axis (κf) and the slow axis (κs). One common way to engineer ρ in fully dense solids is via nanostructuring6, such as fabricating inorganic superlattices7,8,9,10,11 or designing symmetry-breaking crystal architectures in a single material12. However, such engineered materials have relatively small ρ values of less than 20 at room temperature. Conversely, some natural crystalline materials have an intrinsically large ρ (for example, graphite1 and hexagonal boron nitride (hBN)13, with ρ ≈ 340 and 90 respectively), but they are often difficult to process scalably for thin film integration. Some of these films may also lack the electrical or optical properties necessary for functional device applications.

To design materials with higher ρ that are also suitable for real-world applications, an approach needs to be developed to include three key features: (1) a candidate material with intrinsically high κf, usually one with efficient phonon-mediated thermal transport; (2) a method to substantially reduce κs without affecting κf; and (3) facile, scalable production and integration of such a material with precise control of the material dimensions (for example, film thickness). Layered van der Waals (vdW) materials such as graphite and transition metal dichalcogenides (TMDs) provide an ideal material platform for designing such high-ρ materials. They generally have excellent intrinsic in-plane thermal conductivities (κ||) in single-crystalline form. Previous studies have also measured record-low thermal conductivities in nanocrystalline vdW films (for example, WSe2)14,15,16,17 and heterostructures18. One currently missing capability, however, is an approach for significantly decreasing the out-of-plane thermal conductivity (κ) while maintaining high κ||.

TMD films with interlayer rotations

Here we show that such capability is provided by interlayer rotations, as illustrated in Fig. 1a. Interlayer rotation breaks the through-plane translational symmetry at the atomic scale while retaining in-plane long-range crystallinity in each monolayer, thereby providing an effective means for suppressing only κ. For this, we produce large-area TMD films without interlayer registry (referred to here as r-TMD), which possess long-range crystallinity in-plane and relative lattice rotations at every interlayer interface (Fig. 1b). The films are produced in large-scale using two steps: wafer-scale growth of continuous TMD monolayers (polycrystalline; domain size D) and layer-by-layer stacking in vacuum using previously reported methods19,20 (see Methods).

Fig. 1: Structure of r-TMD films.
figure1

a, Conceptual strategy for engineering thermal anisotropy in a single material system, using random interlayer rotation in polycrystalline vdW layered materials. b, Schematic of an r-MoS2 film with random crystalline orientation. c, Greyscale-inverted TEM electron diffraction patterns probed from a 500 nm × 500 nm area of a monolayer and an N = 10 r-MoS2 film. Inset: darkfield TEM image of a monolayer; the scale bar denotes 400 nm and the colours denote different domain orientations from different crystal domains. d, HAADF-STEM image of a cross-section of an N = 10 r-MoS2 film on AlOx coated with Al, with an interlayer spacing of 6.4 Å. e, Large-area MoS2 films transferred onto 1-inch diameter fused silica substrates.

The transmission electron microscopy (TEM) diffraction (Fig. 1c, left) and darkfield (inset) images from a representative MoS2 monolayer show that it comprises large (D ≈ 1 μm), randomly oriented crystalline domains, which connect laterally to form a continuous polycrystalline film. The vacuum stacking generates r-TMD films with a precise layer number (N) and high-quality interfaces20 with interlayer rotation at every stacked interface. The TEM diffraction pattern of N = 10 r-MoS2 (Fig. 1c, right) shows a ring-like pattern due to the significant increase in the number of diffraction spots, emphasizing the random crystalline orientation in the through-plane direction. Clean and well-defined interfaces can be seen from the cross-sectional high-angle annular darkfield scanning TEM (HAADF-STEM) images of r-MoS2 (Fig. 1d and Extended Data Fig. 1; see Methods). The monolayers have a uniform interlayer spacing d ≈ 6.4 Å (see Methods and Extended Data Fig. 2), which is close to the expected value (6.5 Å) for twisted MoS2 multilayers21. Both the growth and stacking steps are scalable, as shown by the optical images of N = 1 and N = 10 r-MoS2 films (~1 cm2) in Fig. 1e and as demonstrated later in Fig. 4. The large-scale uniformity of these films also enables precise and reproducible measurements with minimal spatial variation (Extended Data Fig. 3a–c). In our experiments, r-MoS2 or r-WS2 films with different N (up to 22) are transferred onto a sapphire wafer for the measurements of κ or suspended over a holey TEM grid (Fig. 3a) for the measurements of κ||.

Ultralow out-of-plane conductivity

In Fig. 2, we illustrate κ of r-TMD films, which is measured using time domain thermoreflectance (TDTR; Fig. 2a, inset; see Methods). A stream of laser pulses (pump) heats up the surface of an Al pad deposited on an r-TMD film on sapphire and produces a temperature-sensitive thermoreflectance signal (−Vin/Vout in Fig. 2a), which is measured with a probe pulse after a varying time delay (for cooling). Figure 2a shows three representative curves measured from r-MoS2 with N = 1, 2 and 10. The curves flatten with increasing N, suggesting that heat dissipation slows down significantly. Fitting these curves using a heat diffusion model (solid lines, Fig. 2a) enables us to obtain RTDTR, the total thermal resistance between the Al transducer layer and sapphire across the r-TMD film for different N.

Fig. 2: Through-plane thermal properties of r-MoS2.
figure2

a, TDTR heat dissipation curves of N-layer r-MoS2 films. Inset: TDTR sample geometry. b, Measured thermal resistances across r-TMD films, where the error bars are the TDTR measurement uncertainties. The thermal conductivities for r-MoS2 and r-WS2 are calculated from the slope using the formula RTDTR = R0 + Nd/κ, whereby R0 is the total interfacial thermal resistance. c, Experiment and MD simulation results of κ(T) of MoS2 and r-MoS2 films. The error bars to the MD simulations originate from the simulation uncertainties. The dotted lines connecting the individual data points are guides to the eye. d, LA (top) and TA (bottom) phonon dispersion curves of r-MoS2 along the Γ–A direction. The dotted lines denote the acoustic curves corresponding to bulk MoS2. e, Lifetime of LA and TA phonons parallel to the Γ–A direction in bulk and r-MoS2. The dashed line is the LA mode vibration period derived from the dispersion curve in d.

Source data.

Figure 2b shows RTDTR versus N for r-MoS2 (N ≤ 22; solid circles) and r-WS2 (N ≤ 10; open circles) measured under ambient conditions. We make two observations. First, RTDTR monotonically increases with N. Second, RTDTR varies linearly with N for N ≥ 2. These observations confirm that the through-plane thermal transport in r-TMD films is diffusive in nature, in contrast to the ballistic transport reported in few-layer single-crystalline MoS2 (as thick as 240 nm)22,23. A single parameter κ characterizes the thermal resistance across r-MoS2 (or r-WS2) using the equation RTDTR = R0 + Nd/κ\(,\) where Nd is the total film thickness, and R0 is a constant corresponding to the total interface resistance (r-TMD/Al and r-TMD/sapphire; see Extended Data Table 1). Therefore, we apply linear fitting to the data (N ≥ 2) in Fig. 2b (solid lines) to determine κ of r-MoS2 alone, regardless of the quality and chemical nature of the top and bottom interfaces (see Methods and Extended Data Fig. 3d), which can potentially be altered by metal deposition24,25. We measure κ = 57 ± 3 mW m−1 K−1 for r-MoS2 and κ = 41 ± 3 mW m−1 K−1 for r-WS2, which are similar to the lowest value ever observed in a fully dense solid15 and comparable to the thermal conductivity of ambient air (~26 mW m−1 K−1). These values are approximately two orders of magnitude smaller than those of single-crystalline MoS2 (2–5 W m−1 K−1)26,27 or WS2 (~3 W m−1 K−1)27, despite the r-TMD films having the same chemical composition as their bulk counterparts as well as clean interfaces (Fig. 1d). This strongly suggests that the main difference, the interlayer rotation, is the principal cause for the ultralow κ in these r-TMD films. Furthermore, repeating similar TDTR experiments on r-MoS2 at different temperatures (T) produces a relatively flat κ(T) curve (green stars, Fig. 2c), a behaviour different from the decreasing κ with T seen in bulk MoS2 (blue squares, lower Fig. 2c).

To understand the microscopic mechanisms that give rise to the dramatic reduction in κ, we carry out homogeneous non-equilibrium molecular dynamics (HNEMD) simulations for the model structures of r-MoS2 and bulk MoS2 (see Methods and Extended Data Table 2)28,29,30. Figure 2c shows κ|| and κ of r-MoS2 (solid circles) and bulk MoS2 (empty circles) calculated from our molecular dynamics (MD) simulations at different temperatures. The calculated κ drops by a factor of more than 20, from 3.7 ± 0.5 W m−1 K−1 in bulk MoS2 to 0.16 ± 0.04 W m−1 K−1 in r-MoS2 at 300 K, and also does not decrease with T, suggesting a transition away from the phonon-limited thermal transport mechanism observed in bulk MoS2.

Further analysis of the vibrational spectrum of r-MoS2 enables us to break down the reduction in κ in terms of the changes in the group velocities (vg) and lifetimes (τ), which are the two factors that determine the thermal conductivity according to Boltzmann transport theory. Figure 2d shows that the vg of the through-plane longitudinal acoustic (LA) mode in r-MoS2 remains similar to that of bulk MoS2 (dashed lines), but the transverse acoustic (TA) modes in r-MoS2 undergo extreme softening with their vgs practically vanishing31. This implies a loss of resistance with respect to lateral shear, consistent with the low-frequency Raman spectra of r-MoS2 films (see Methods and Extended Data Fig. 4) and previous calculations32,33. In addition, the τ of both the LA and the TA modes (Fig. 2e) in r-MoS2 are more than one order of magnitude smaller than in bulk MoS2, with the LA lifetimes being close to the period of the LA mode vibration (dashed line). From these results, the median mean free path \(\widetilde{l}\) = vgτ for the LA modes is estimated to be 2 nm, suggesting that the heat-carrying LA modes are strongly scattered and that a larger D is unlikely to significantly affect κ since D >> \(\widetilde{l}\). Overall, the strongly suppressed TA modes, indicating a loss of resistance to lateral shear, and the overdamping of the LA modes as the main heat carriers, lead to extremely inefficient thermal transport along the through-plane direction in r-MoS2. Along with the nearly temperature-independent κ, this result suggests a glass-like conduction mechanism.

In-plane conductivity and anisotropy

In contrast to κ, κ|| remains high in our simulations with only a modest reduction compared to the ideal bulk crystal (less than a factor of two at 300 K; Fig. 2c). This is indeed what we observe in our Raman thermometry experiments as discussed in Fig. 3 (see Methods). We direct a focused laser spot (λ = 532 nm) at the centre of a suspended r-MoS2 film (Fig. 3a; hole diameter of 5 μm, at 15 torr), which increases its temperature (ΔT) locally upon absorbing laser power Pabs. ΔT is then measured using the temperature-sensitive Raman peak shift (Δω) using a sensitivity factor (|dω/dT|) independently measured for each sample. Examples of Raman spectra measured for N = 2 are shown in Fig. 3b.

Fig. 3: In-plane thermal properties and thermal anisotropy of r-MoS2 films.
figure3

a, 45° SEM micrograph of an N = 4 r-MoS2 film suspended on a TEM grid for Raman thermometry. b, Raman spectra of an N = 2 r-MoS2 film with different absorbed laser powers. Inset: Raman thermometry sample geometry. c, A1g Raman peak shifts versus power absorbed by r-MoS2 films of various N. Inset: layer-dependent thermal conductance values (absorbed power divided by temperature increase) in domain size D = 1 μm and D = 400 nm r-MoS2 films. The error bars are the propagated uncertainties from the calculation of the conductance value for each N. d, Comparison of ρ (y axis), κs (x axis), and κf (diagonal dashed lines) measured for different anisotropic thermal conductors. r-MoS2 has an ultrahigh ρ close to 900, which is larger than bulk MoS2, PG, and disordered layered WSe2. The error bar for ρ of r-MoS2 comes from the propagated uncertainties of the calculated κ and κ|| values.

Source data.

Figure 3c plots Δω versus Pabs for r-MoS2 with different N (2 to 5). The slope of the linear fit (|d(Δω)/dPabs|), which is inversely proportional to the in-plane thermal conductance of the film, is plotted in the inset (solid dots; D ≈ 1 μm). We again observe a linear relation, which indicates that κ|| is well defined for r-MoS2 independent of N, similar to the case of κ. Using a simple diffusion model with radial symmetry (see Methods and Extended Data Fig. 5c for calculation details and other input measurements), we calculate a high κ|| value of 50 ± 6 W m−1 K−1. This value is similar to the predictions of our MD simulations (Fig. 2c) and consistent with previous reports of Raman thermometry on single-crystalline monolayer MoS2 (35–84 W m−1 K−1) at room temperature34,35,36,37. The κ|| of these r-MoS2 films is close to the intrinsic phonon-limited value despite the films being made of polycrystalline monolayers. This result is further supported by our additional measurements on continuous r-MoS2 films with a smaller D ≈ 400 nm (open dots, dashed lines, Fig. 3c inset; Extended Data Fig. 5e). The measured value of κ|| ≈ 44 ± 6 W m−1 K−1 is within the margin of error of that of the D ≈ 1 μm films. This suggests that the phonon mean free path is smaller than 400 nm, which is consistent with previous reports23,38,39,40,41,42. Furthermore, the measured in-plane conductance decreases with T (Extended Data Fig. 6a). This further confirms the phonon-mediated thermal transport mechanism in-plane, in contrast to the glass-like thermal conduction along the through-plane direction.

Our experiments and calculations confirm that interlayer rotation in r-TMD films results in highly directional thermal conductivity and a direction-dependent thermal conduction mechanism. The rotation significantly reduces κ while maintaining high κ||, leading to an ultrahigh value of ρ. We estimate ρ ≈ 880 ± 110 at room temperature for the r-MoS2 films, higher than that of pyrolytic graphite (PG), which is considered to be one of the most anisotropic thermal conductors (ρ ≈ 340)43. In Fig. 3d, we compare our result with other previously reported values of ρ in phonon-based solids15,27,43,44 (for a full comparison, see Extended Data Fig. 6b). Compared to a bulk MoS2 crystal (ρ ≈ 20)27 or disordered layered WSe2 (ρ ≈ 30)15,44, r-MoS2 has a significantly larger ρ because interlayer rotation reduces only κ, as denoted by the grey arrow parallel to the equi-κf lines. This also suggests that ρ can be made even larger by starting with the monolayers of a layered vdW material with a higher κ|| value such as graphene.

Anisotropic vdW heat diffuser

In Fig. 4, we show that the extreme anisotropy of our r-MoS2 films can lead to excellent heat dissipation in-plane from a heat source and drastic thermal insulation in the through-plane direction. Using the COMSOL software, we perform thermal finite-element simulations of a 10-nm-thick r-MoS2 film draped over a nanoscale Au electrode (15 nm tall, 100 nm wide) on a 50 nm SiO2/Si substrate (Fig. 4a). Our simulation results show that for a fixed power of 8 mW supplied to the Au electrode (near thermal breakdown), the temperature rise ∆T of the Au electrode covered by r-MoS2 is 50 K lower than that of the bare electrode, thereby demonstrating our film’s effectiveness at spreading heat due to its excellent κ|| (Fig. 4b, c). Interestingly, the extreme thermal anisotropy of our r-MoS2 films provides thermal insulation in the through-plane direction, with much lower MoS2 surface ∆T values that are only one-third of the value of the bare Au electrode. While single-crystal MoS2 displays similar properties, the insulation effect is stronger in r-MoS2 (Extended Data Fig. 7a). This implies that heat is efficiently directed away from the hot Au electrode laterally through r-MoS2 but not to the surface of r-MoS2, making the surface of the entire device significantly cooler.

Fig. 4: Temperature profiles and heat spreader efficiencies of r-MoS2 films on Au electrodes.
figure4

a, Schematic of the sample configuration of r-MoS2 draped across a current-carrying Au electrode that is 100 nm wide, 15 nm thick and 10 μm long. b, Thermal finite element modelling results of Au electrodes (bare, covered with 10-nm-thick r-MoS2) at constant heating power of 8 mW supplied through Joule heating. c, Lateral profiles of temperature increases across the Au/SiO2 surface (solid dots) and on the r-MoS2 top surface (open circles). Insets: cross-sectional temperature distribution of Au electrodes with and without r-MoS2, using the same colour scale as in b. d, I–V curve of an Au electrode, with and without N = 16 r-MoS2. Inset: optical micrograph of six fabricated Au electrodes. e, Histogram of Ic of Au electrodes with and without an N = 16 r-MoS2 heat spreader and their median values.

Source data.

Our experiments corroborate these simulation results. For this, we fabricate nanoscale Au electrodes with the same geometry and substrate as in our simulation (image shown in Fig. 4d, inset) and transfer N = 16 r-MoS2 (~10 nm thick) using the vacuum stacking process. Both bare and coated Au electrodes show similar resistance at low currents. At higher currents, current-induced Joule heating leads to the thermally activated electromigration process, which causes the electrodes to fail45. Figure 4d compares representative current–voltage (I–V) curves measured from a bare and coated Au electrode, which shows that the Au electrode with r-MoS2 can carry a larger current without breaking. The histogram of critical current Ic (maximum current a Au electrode sustains for at least 20 s) measured from 20 electrodes (10 bare and 10 with r-MoS2) reveals a ~50% increase in the median Ic values (Fig. 4e). These results demonstrate our r-MoS2 film’s ability to efficiently dissipate Joule heat and keep the electrodes cool, as our simulation predicts. As the electromigration process is dominated by the temperature, the observed increase of Ic and maximum power before breaking is in good agreement with our simulation in Fig. 4c. Furthermore, we note that the r-MoS2 film can be integrated with the Au electrodes using mild conditions that do not affect their electrical properties (Extended Data Fig. 7b).

Outlook

We expect interlayer rotation to be an effective and generalizable way to reduce κ and potentially engineer anisotropic thermal properties in a variety of layered materials. Our results call for a systematic study of the exact relation between κ and rotation angle, which could reveal unexpected relationships analogous to the studies of electrical transport in twisted bilayer graphene46. Interlayer rotations can be combined with other parameters (such as pressure or interlayer spacing47,48) and advanced structures (superlattices and heterostructures18) to realize highly tunable ρ, allowing for the customization of thermal transport properties with an unprecedented level of directional and spatial control.

Methods

Sample preparation

Large-area, polycrystalline transition metal dichalcogenide (TMD) (MoS2 and WS2) monolayers were grown on SiO2/Si substrates in a hot-walled tube furnace via metal-organic chemical vapour deposition adapted from a previously reported protocol19. The growth conditions were optimized to produce high-quality monolayer materials with structural characteristics necessary for thermal measurements. These characteristics include large domain size (D ≈ 1 μm and 0.4 μm), full monolayer coverage, and laterally stitched grain boundaries.

Briefly, Mo(CO)6 and W(CO)6 (diluted in N2 to 15 torr) were used as the metal precursors for the MoS2 and WS2 growths, respectively. (C2H5)2S was used as the chalcogen source. All precursors were kept at room temperature. N2 and H2 were used as carrier gases. Typical growth times were 15–20 h for MoS2 at a growth temperature of 525 °C. Typical growth times for WS2 were 2 h at a temperature of 650 °C.

To make the r-TMD films, a TMD monolayer was spin-coated with PMMA A8 (poly-methyl methacrylate, 495 K, 4% diluted in anisole) at 2,800 rpm for 60 s, then baked at 180 °C for 3 min. The PMMA-coated monolayer was stacked onto TMD monolayers layer by layer to a target layer number (N) and transferred to the desired substrates using a previously reported programmed vacuum stacking method20.

TDTR samples

The stacked r-TMD films were transferred onto sapphire substrates (Valley Design, C-plane), which were cleaned with Nanostrip solution for 20 min at 60 °C and then rinsed with deionized water. The PMMA layer on the film was removed by immersing the entire substrate in acetone at 60 °C for 1 h. The film was annealed under a 400/100 SCCM Ar/H2 environment at 350 °C for 4 h. After cleaning, ~90-nm-thick, 90 μm × 90 μm Al pads were deposited onto the TMD films through a holey TEM grid shadow mask using electron-beam evaporation.

Raman thermometry samples

Raman experiments were performed on a different set of films from the TDTR-measured films. First, holey SiNx transmission electron microscopy (TEM) grids were cleaned in a N2/H2 plasma at 100 °C and 180 mtorr for 3 min, followed by the transfer of stacked MoS2 films onto the TEM grids. During the transfer process, the PMMA-coated r-MoS2 was suspended on holey thermal release tape before contacting the TEM grid. The extra PMMA-MoS2 not on the TEM grid was cut away at 180 °C so the PMMA layer was softened. PMMA was removed from r-MoS2 on the TEM grid via annealing the film in 400/100 SCCM Ar/H2 at 350 °C for 4 h.

Cross-sectional STEM

The N = 10 films were coated with Al that was electron beam evaporated onto the surface, whereas the top surface of N = 20 films was bare. The r-MoS2 cross-section was prepared using a Thermo Scientific Strata 400 focused ion beam. Protective layers of carbon (~200 nm) and platinum (~1 µm) were deposited on the sample. A cross-section was milled at a 90° angle from the sample using a Ga ion beam at 30 kV. The cross-section was then polished to ~150 nm thickness with the ion beam at 5 kV.

The cross-section was imaged in a Thermo Scientific Titan Themis scanning transmission electron microscope at 120 kV with a probe convergence angle of 21.4 mrad. The N = 10 film was imaged at a beam voltage of 120 kV, whereas the N = 20 film was imaged at 300 kV. All images were analysed using the open-source software Cornell Spectrum Imager49. The high-angle annular dark field (HAADF) image of the sample (see 'TMD films with interlayer rotation' in the main text; Fig. 1d) shows, from top to bottom, the Al crystal lattice along the [110] zone axis, ten layers of MoS2 (bright bands), followed by an AlOx layer.

TDTR

We used TDTR to measure the thermal conductivity of our r-TMD films. We used a mode-locked Ti:sapphire laser, which produced a train of pulses at a repetition rate of 74.86 MHz, with wavelength centred at 785 nm and a total power of 18 mW. The steady-state temperature rise at the surface of the samples was <4 K for all temperatures. For the low temperature TDTR measurements, an INSTEC stage was used with liquid nitrogen cooling; the other beam conditions were the same. The laser beam was split into pump and probe beams. A mechanical delay stage was used to delay the arrival of the probe with respect to the pump on the sample surface by changing their optical path difference, before they were focused onto the sample surface through an objective lens. The 1/e2 radius of the focused laser beams was 10.7 μm. For our measurements, we modulated the pump beam at a frequency of 9.3 MHz so that the thermoreflectance change at the sample surface could be detected by the probe beam through lock-in detection. The ratio of the in-phase and out-of-phase signals from the lock-in was fitted to a thermal diffusion model. The full details of the TDTR measurement can be found elsewhere50,51.

Calculation of κ

The modelling required material parameters such as heat capacity (C), thickness (h), interface conductance (G) and thermal conductivity (κ) for each layer. Our TDTR samples have three chemically distinct layers with the following structure (from the top): Al/r-TMD/sapphire. In our fitting process, the heat capacities of all materials were adopted from literature52. The thickness of Al layer was obtained from picosecond acoustics using a longitudinal speed of sound of 6.42 nm ps−1 (Extended Data Fig. 3e). The thickness of the r-TMD film was calculated from the product of N and the interlayer spacing (d). The latter was measured by performing grazing-incidence wide-angle X-ray spectroscopy (GIWAXS; see GIWAXS section below in the Methods and Extended Data Fig. 2) on the r-TMD films, which gave d ≈ 0.64 nm. The total thicknesses of the r-MoS2 films were <15 nm; thus, this layer was treated as part of the Al–sapphire interface as a single thermal layer characterized by a single thermal conductance value G1. We used the bulk value of the volumetric heat capacity of 1.89 J K−1 cm−3 for the r-MoS2 layer. The thermal conductivity of the Al layer was calculated from the Wiedemann–Franz law using the electrical resistance of a transducer layer deposited on a bare sapphire substrate as a reference sample. The thermal conductivity of the sapphire substrate, 35 W m−1 K−1, was measured using the same reference sample. Thus, the only remaining free parameter to fit for was G1. To obtain κTMD from G1, we perform TDTR on various N-layer TMD films, then perform a linear fit on the effective thermal resistance (RTDTR, equal to 1/G1) versus N data points; the slope of the linear fit is inversely proportional to the thermal conductivity, whereas the y intercept yields the total interfacial thermal resistances (R0) of the top and bottom interfaces. In Extended Data Table 1, our R0 values match the values reported in literature22,27,53. We note that, although R0 changes depending on the chemical nature of the metal–TMD interface, the slope of the RTDTRN plot (which is used to extract κ) remains constant, despite the use of different transducer metals, as illustrated in Extended Data Fig. 3d.

For highly anisotropic materials, the anisotropy ratio of an in-plane thermal conductivity to a through-plane conductivity should be included in the thermal model. Despite the ultrahigh thermal anisotropy expected of our r-TMD films, our through-plane thermal conductivity measurements were probably not sensitive to the thermal conductivity anisotropy given the thinness of our r-TMD films. Hence, we assumed a one-dimensional thermal transport model and neglected the in-plane thermal transport in our calculations. We found that the effect of the anisotropy was significant only at a smaller modulation frequency (f = 1.12 MHz) and 1/e2 beam radius of ~3.2 μm, and so we deliberately chose a larger f and a 1/e2 beam radius to reduce the sensitivity of our TDTR signal to the in-plane thermal transport.

Raman thermometry

We followed a similar procedure from previous reports34,54,55 with the modification of lower pressures during measurement. All the Raman measurements were performed using a Horiba Raman spectrometer with a laser excitation wavelength of 532 nm and a long-working distance, 50× objective lens (numerical aperture (NA) = 0.5). The r-MoS2 A1g peak shift (ω) versus temperature (T) relation was calibrated using a temperature-controlled, low-vacuum-compatible Linkam stage. For all our Raman measurements, we used the A1g peak since this out-of-plane vibrational mode is less sensitive to in-plane strain56. The ωT calibration measurements were performed at atmospheric pressure and with low laser powers. The stage was purged with dry N2 gas throughout the calibration step to prevent oxidative damage to the film at high temperatures. Extended Data Fig. 5d shows representative ωT calibration curves for N = 2 and N = 4 r-MoS2 films, where a linear fit was performed to obtain the temperature-dependent Raman coefficients. This process was repeated for r-MoS2 films with different domain sizes D (400 nm and 1 μm) for N = 1–3 (Extended Data Fig. 5e).

To measure the in-plane thermal conductivity (κ||) of our films, the laser power (P) was varied and the corresponding Δω values were recorded. The in-plane thermal conductance was obtained from the reciprocal of the slope of the ΔωP linear fit, which is illustrated for r-MoS2 films with N = 2–5 and D = 1 μm in Fig. 3c and for r-MoS2 films with N = 1–3 and D = 400 nm in Extended Data Fig. 5b. As thermal conductivity changes with temperature, laser powers were kept below 250 μW to induce a relatively small ΔT in the film and ensure that the value of κ|| remained relatively constant. This was verified from the observation of a linear ΔωP regime for P < 250 μW. Any higher laser powers caused the Δω–P curve to deviate from the linear regime with \(\frac{{{\rm{d}}}^{2}\omega }{{\rm{d}}{P}^{2}} < 0\). This indicates that the local film temperature increased faster at higher P > 250 μW, which signified that the thermal conductivity could no longer be assumed to be constant. Instead, the thermal conductivity decreased with increasing temperature, consistent with the T-dependent Raman measurements.

The ΔωP measurements were conducted at a pressure of 15 torr to eliminate any heat loss to air. We verified that a lower pressure down to 4 mtorr gave rise to similar Δω values as the measurements at 15 torr (Extended Data Fig. 5a), weighted by the beam spot size.

The other relevant input quantities for our thermal calculations were obtained as follows: the beam spot radius (r0) was estimated using the knife-edge method, whereby a one-dimensional Raman map was taken across a gold step edge on an Au-patterned silicon chip, and the spatial distribution of the integrated peak intensities was fitted to an error function. We measured r0 = 0.71 ± 0.09 μm. The laser powers were measured using a Thorlabs standard silicon photodiode power sensor. The r-MoS2 absorbance \(A=\frac{{\rm{Absorbed\; light\; intensity}}}{{\rm{Incident\; light\; intensity}}}\) was measured at room temperature on a white-light microscope with a 532 nm band-pass filter and a low-NA condenser aperture. We measured the light intensity transmitted through and reflected from a r-MoS2 film suspended on a TEM grid, then compared it against a blank TEM grid. The data were collected using a 12-bit SensiCam QE CCD camera. The pixel intensities were analysed using ImageJ. The values for A were calculated using the formula A = 1 − T − R. We measured A(N) for N = 1–5, then fitted A to a power law. A(N) was found to follow the relation \(A=1-{0.92}^{N}\) (Extended Data Fig. 5c), which matched previous reports20. We use the value A measured at room temperature for our Raman analysis.

Calculation of κ ||

To obtain the value of κ||, we used the two-dimensional thermal diffusion equation with a radial symmetry, following previous reports of Raman thermometry of two-dimensional films54,55. We assumed a Gaussian laser profile \(q(r)=\frac{{PA}}{({\rm{\pi }}{{r}_{0}}^{2})t}\,{\rm{\exp }}\left(-\frac{{r}^{2}}{{{r}_{0}}^{2}}\right)\). We solved for κ|| numerically using the following equations:

$${\kappa }_{{\rm{||}}}\frac{1}{r}\frac{{\rm{d}}}{{\rm{d}}r}\left(r\frac{{\rm{d}}{T}_{{\rm{susp}}}(r)}{{\rm{d}}r}\right)+q(r)=0{\rm{;}}r < R$$
$${\kappa }_{{\rm{||}}}\frac{1}{r}\frac{{\rm{d}}}{{\rm{d}}r}\left(r\frac{{\rm{d}}{T}_{{\rm{supp}}}(r)}{{\rm{d}}r}\right)-\frac{G}{t}[{T}_{{\rm{supp}}}(r)-{T}_{{\rm{a}}}]=0{\rm{;}}r > R$$

applying the boundary conditions

$${T}_{{\rm{susp}}}(R)={T}_{{\rm{supp}}}(R)$$
$$\frac{{\rm{d}}{T}_{{\rm{susp}}}}{{\rm{d}}r}(R)=\frac{{\rm{d}}{T}_{{\rm{supp}}}}{{\rm{d}}r}(R)$$

where \(r\) is the distance from TEM hole centre, \(P\) is the laser power, \(t\) is the film thickness, \(R\) is the TEM hole radius, \(T\) is the film temperature where \({T}_{{\rm{susp}}},{r}\le R\) and \({T}_{{\rm{supp}}},{r}\ge R\), \({T}_{{\rm{a}}}\) is the ambient temperature, \(A\) is the fraction of laser power absorbed, and G = 10 MW m−2 K−1 is the interfacial thermal conductance between r-MoS2 and SiNx.

We solved for the expression of T(r) and obtained an expression for the average temperature measured by the Raman shift

$${T}_{{\rm{m}}}=\frac{{\int }_{0}^{R}T(r)r\,\times \,\exp (-\frac{{r}^{2}}{{{r}_{0}}^{2}}){\rm{d}}r}{{\int }_{0}^{R}r\,\times \,\exp (-\frac{{r}^{2}}{{{r}_{0}}^{2}}){\rm{d}}r}$$

κ|| was obtained by substituting the experimentally measured value for \({T}_{{\rm{m}}}\) and solving the above equation numerically for κ||. We calculated κ|| for each N, and we reported the average value in the main text.

The total measurement uncertainty reported in the main text was calculated based on the error assessment for individual parameters. We used an approximate analytical solution

$${\triangle T=T}_{{\rm{m}}}-{T}_{{\rm{a}}}\approx -\left(\frac{{P}_{{\rm{abs}}}}{{\kappa }_{{||}}}\right)\left(\frac{1}{2{\rm{\pi }}{rt}}\right){\rm{ln}}\left(\frac{R}{{r}_{0}}\right)$$
$${\kappa }_{{\rm{||}}}\approx \left({\left(\frac{\triangle \omega }{\triangle P}\right)}^{-1}\left(\frac{\triangle \omega }{\triangle T}\right)\right)\left(\frac{{A}_{0}}{2{\rm{\pi }}d}\right){\rm{ln}}\left(\frac{R}{{r}_{0}}\right)$$

where \(\omega \) is the Raman frequency of A1g peak, \({A}_{0}\) is the absorption of the monolayer, and \(d\) is the thickness of a monolayer. The difference between the full numerical solution and this analytical form is below 3%. We identified the following independent quantities that carry uncertainty for consideration in our overall uncertainty estimation of κ||.

  1. 1.

    \({\frac{\triangle P}{\triangle T}=(\frac{\triangle \omega }{\triangle P})}^{-1}\left(\frac{\triangle \omega }{\triangle T}\right)\): the associated uncertainty was derived from the error in the linear fit of \((\frac{\triangle \omega }{\triangle P})\) and \(\left(\frac{\triangle \omega }{\triangle T}\right)\) for every sample measured. The total uncertainty in the average \(\frac{\triangle P}{\triangle T}\) value was 9% for D = 1 μm and 8% for D = 400 nm.

  2. 2.

    \(\left(\frac{{A}_{0}}{2{\rm{\pi }}{d}_{0}}\right)\): the uncertainty in A0 from the A(N) fit was 4%.

  3. 3.

    \({\rm{ln}}\left(\frac{R}{{r}_{0}}\right)\): the uncertainty originated from the uncertainty in r0. From 14 repeated measurements of the beam spot size using the knife-edge method, we calculated the standard deviation of r0 to be 12%, which translated to an uncertainty in the expression \({\rm{ln}}\left(\frac{R}{{r}_{0}}\right)\) to be 9%.

Total uncertainty in κ||: 13% (for both D = 1 μm and D = 400 nm).

Variable pressure Raman thermometry measurements

Previous Raman thermometry measurements on graphene films57 and carbon nanotubes58 had shown an appreciable difference between measurements performed in air and at lower pressures, as well as in different gaseous environments. We extended the same precaution and repeated our Raman thermometry measurements at low pressures to reduce heat dissipation to air, an extra heat loss channel that would lead to an overestimation of the thermal conductivity of the r-MoS2 films.

Our ΔωP measurements in Fig. 3 were conducted at a pressure of 15 torr. In Extended Data Fig. 5a, we compared the Δω values of N = 2 r-MoS2 at three different P (1 atm, 15 torr and 4 mtorr), after correcting for the different laser spot sizes.

Temperature-dependent Raman thermometry for κ ||

We performed Raman thermometry while varying the ambient temperature Ta using a Linkam stage. No oxidation or sample damage was detected for any of the temperatures used. We performed the same Δω–P measurements and calculated the κ|| value of for each Ta. We plot a κ|| versus T curve, where the x axis is T = Ta (Extended Data Fig. 6a).

We note that the measured values of κ|| here were lower than the room temperature values reported in the main text. We ascribe this to the sub-optimal growth conditions for the constituent monolayers used for this sample.

r-MoS2 heat spreader experiments (electromigration of Au nanoelectrodes)

All Au electrodes were fabricated on Si substrates with 50 nm dry SiO2 in three nanopatterning and deposition layers: (A) the nanoelectrodes (10 μm long, 100 nm wide, 15 nm thick); (B) the contact pads that would interface with the external electronics (200 μm long, 300 μm wide, 100 nm thick); and (C) the leads connecting the nanoelectrodes and the contact pads (~1,000 μm long, 50 μm wide, 15 nm thick).

We first defined the leads (B) and then the contact pads (C) using standard photolithography, and electron-beam evaporation of Ti (1 nm)/Au and lift-off. The final step was defining the nanoelectrodes (A) using electron-beam lithography, deposition of 15 nm Au, and lift-off.

Electron-beam lithography

We used a bilayer of resists: copolymer P(MMA-MAA 11%) in ethyl lactate and 950 K PMMA A4. The writing was executed with a Raith EBPG 5000 Plus E-beam writer with the beam conditions of 25 nA current, dose of 1,200 μC cm−2, 300 μm aperture size, 100 kV accelerating voltage.

Film transfer

After the nanoelectrodes, leads and pads were fabricated, the device was cleaned with an O2 plasma for 30 s to remove any resist residue and to promote adhesion of the r-MoS2 film to the Au electrodes and the SiO2 surface. A PMMA coated N = 16 r-MoS2 film was transferred onto the electrodes using the same process as the stacking method as outlined above. The PMMA on the r-MoS2 film was removed by immersing the entire chip in toluene at 60 °C for 1 h.

Electrical measurements

All measurements were performed in ambient conditions with a home-built probe station in a two-probe geometry. To measure Ic in Fig. 4d, we swept the voltage bias in only one direction at a rate such that the rate in current increase is 0.05 mA per 20 s.

For comparison, we deposited SiNx onto Au electrodes (10 μm long, 10 nm thick, 100 nm wide) using plasma-enhanced chemical vapour deposition with the following conditions: 10 s deposition at 90 °C and 10 torr and 1,000 W plasma power, with 25 SCCM/35 SCCM SiH4 and N2 as the precursors. The film thickness was measured via ellipsometry to be 16 nm.

Computational methodology

Structural models

Structural models were created according to an algorithm previously described in literature59, which was implemented in Python using the atomic simulation environment package60. The structure models were subsequently relaxed using an analytic bond-order potential61 and implemented in the LAMMPS package62.

The main r-MoS2 model used in the simulations described here comprised 10 randomly stacked layers with a total of 10,152 atoms. The 10 layers came in pairs; each pair was related by a 60° rotation. The four primitive angles present in the stack are 16.1, 25.28, 34.72 and 43.9°. Due to strain, each layer contained a different number of atoms in accordance with strains of around 10%.

The bulk structure used in the MD simulations comprised 40 layers (20 conventional unit-cells) with a total of 26,880 atoms and cell vectors of 44.44, 43.98 and 243.57 Å.

HNEMD simulations

The interatomic potential in our simulations produces the expected slight increase in the interlayer spacing in r-MoS2 and yields thermal conductivities and phonon dispersions of bulk MoS2 that agree with previous experimental observations and Boltzmann transport calculations based on density functional theory38, confirming our MD model’s suitability for this study. The structures described above were driven by an optimized driving force (Extended Data Fig. 8a) and subsequently relaxed, after which the thermal conductivity was computed using HNEMD simulations63 and implemented in the graphics processing units molecular dynamics (GPUMD) package29. We also included the effects of thermal expansion in the simulations (Extended Data Fig. 8b). The calculated κ values of r-MoS2 are higher than the experimental values. We attribute the discrepancy to our neglecting any quantum effects and all boundary scattering in our simulations. Including such effects could further improve results from simulations. Statistics and averages were gathered from ten independent simulation runs for each system and temperature. The other parameters used in these simulations are compiled in Extended Data Table 2.

Phonon dispersion and lifetimes

We first generated the bulk MoS2 phonon dispersion in the harmonic limit as a reference to phonon dispersion calculations using MD simulations. We computed the harmonic (0 K) phonon dispersion using the PHONOPY package64. Forces were computed for 6 × 6 × 2 supercells using the LAMMPS code. Lifetimes were calculated using the lowest applicable order of perturbation theory using the PHONO3PY package65, which also provided us with the thermal conductivity as obtained from a direct solution of the Boltzmann transport equation66. In these calculations, the Brillouin zone was sampled using a 10 × 10 × 10 Γ-centred q-point mesh, which was chosen for consistency with the supercell used in the HNEMD simulations.

Next, we compared the dispersion of bulk MoS2 to that calculated using MD simulations to verify the accuracy of our MD simulations for calculating the phonon dispersion of r-MoS2. For both bulk and r-MoS2, we extracted the phonon dispersions and lifetimes at 300 K by analysing the longitudinal and transverse current correlation functions generated by MD simulations in the microcanonical (NVE) ensemble using dynasor67. The MD simulations details were otherwise identical to the HNEMD simulations. The obtained correlation functions were Fourier transformed and fitted to peak shape functions corresponding to (over)damped harmonic oscillators using the full expressions given in the dynasor paper67 to obtain phonon frequencies and lifetimes.

Finite-element analysis

We used the COMSOL software to simulate the steady-state temperature distribution in a Au electrode on a SiO2/Si substrate. Our geometry contains an Au electrode that is 100 nm wide, 15 nm thick, and 10 μm long on a 50-nm-thick SiO2 layer on a Si substrate. We layer a 10-nm-thick MoS2 film onto the Au electrode and the SiO2 layer. For the thermal anisotropy consideration, we define the thermal conductivity slow axis direction to always be perpendicular to the film’s bottom surface in contact with the substrate or the Au electrode, including the Au electrode side walls.

We supply the Au electrode with 8 mW uniformly over the entire volume as the heat source, matching the power conditions at which the Au electrode fails in our experiments. As the boundary condition, we set the bottom surface on the Si substrate to be at 293.15 K. We also account for all the interfacial thermal resistances between heterogeneous surfaces in our calculations, which include r-MoS2/Au, r-MoS2/SiO2, Au/SiO2, SiO2/Si (refs. 68,69,70,71). All effects of radiation are neglected as they do not affect the temperature values in our simulations.

Low-frequency Raman measurements

The low-frequency Raman spectra of N = 2, 3 and 4 r-MoS2 films, along with the spectrum for MoS2, are shown in Extended Data Fig. 4a. From the layer-dependence of the peak positions, we assigned these to be the breathing modes of MoS2 (ref. 72). We did not observe any peaks corresponding to shear modes in r-MoS2. Our findings agree with theoretical studies of low-frequency Raman modes of twisted MoS2 bilayers, which showed that the shear mode peaks redshift to below the detection capabilities (2 cm−1)73. The positions of the breathing mode peaks of the r-MoS2 films were close to those in exfoliated few layer MoS2 (Extended Data Fig. 4b) as reported in literature72,74. This observation agreed with our MD simulations that suggested that the transverse vibrational mode was suppressed by interlayer rotation, while the longitudinal vibrational mode was retained.

GIWAXS

The GIWAXS measurement was performed using SAXSLAB (XENOCS)’s GANESHA (lab-source Cu Kα, photon flux ~108 photons s−1) to characterize the interlayer spacing of r-MoS2 films. An N = 10 r-MoS2 film was prepared on a SiO2/Si substrate. The incidence angle of the X-ray beam was 0.2° and the integration time was ~60 s. Radially integrating the two-dimensional diffraction images along the out-of-plane direction produced the diffraction spectrum along the c axis shown in Extended Data Fig. 2.

The peak position of 14° corresponded to an interlayer spacing of 6.4 Å in the [0 0 1] direction, which matched previous reports of r-TMD films20.

Data availability

The data that support the findings of this study are available from the corresponding authors on reasonable request. Source data are provided with this paper.

Code availability

All code used in this work is available from the corresponding authors on reasonable request.

References

  1. 1.

    Minnich, A. J. Exploring the extremes of heat conduction in anisotropic materials. Nanoscale Microscale Thermophys. Eng. 20, 1–21 (2016).

    ADS  CAS  Google Scholar 

  2. 2.

    Cui, Y., Li, M. & Hu, Y. Emerging interface materials for electronics thermal management: experiments, modeling, and new opportunities. J. Mater. Chem. C 8, 10568–10586 (2020).

    CAS  Google Scholar 

  3. 3.

    Moore, A. L. & Shi, L. Emerging challenges and materials for thermal management of electronics. Mater. Today 17, 163–174 (2014).

    CAS  Google Scholar 

  4. 4.

    Norley, J., Tzeng, J. J.-W., Getz, G., Klug, J. & Fedor, B. The development of a natural graphite heat-spreader. In Seventeenth Annual IEEE Semiconductor Thermal Measurement and Management Symposium 107–110 (IEEE, 2001).

  5. 5.

    Chung, D. D. L. & Takizawa, Y. Performance of isotropic and anisotropic heat spreaders. J. Electron. Mater. 41, 2580–2587 (2012).

    ADS  CAS  Google Scholar 

  6. 6.

    Termentzidis, K. Thermal conductivity anisotropy in nanostructures and nanostructured materials. J. Phys. Appl. Phys. 51, 094003 (2018).

    ADS  Google Scholar 

  7. 7.

    Lee, S.-M., Cahill, D. G. & Venkatasubramanian, R. Thermal conductivity of Si–Ge superlattices. Appl. Phys. Lett. 70, 2957–2959 (1997).

    ADS  CAS  Google Scholar 

  8. 8.

    Luckyanova, M. N. et al. Anisotropy of the thermal conductivity in GaAs/AlAs superlattices. Nano Lett. 13, 3973–3977 (2013).

    ADS  CAS  Google Scholar 

  9. 9.

    Borca-Tasciuc, T. et al. Thermal conductivity of symmetrically strained Si/Ge superlattices. Superlattices Microstruct. 28, 199–206 (2000).

    ADS  CAS  Google Scholar 

  10. 10.

    Lee, W.-Y. et al. Anisotropic temperature-dependent thermal conductivity by an Al2O3 interlayer in Al2O3/ZnO superlattice films. Nanotechnology 28, 105401 (2017).

    ADS  Google Scholar 

  11. 11.

    Rawat, V., Koh, Y. K., Cahill, D. G. & Sands, T. D. Thermal conductivity of (Zr,W)N/ScN metal/semiconductor multilayers and superlattices. J. Appl. Phys. 105, 024909 (2009).

    ADS  Google Scholar 

  12. 12.

    Sun, B. et al. Dislocation-induced thermal transport anisotropy in single-crystal group-III nitride films. Nat. Mater. 18, 136–140 (2019).

    ADS  CAS  Google Scholar 

  13. 13.

    Jiang, P., Qian, X., Yang, R. & Lindsay, L. Anisotropic thermal transport in bulk hexagonal boron nitride. Phys. Rev. Mater. 2, 64005 (2018).

    CAS  Google Scholar 

  14. 14.

    Chiritescu, C. et al. Low thermal conductivity in nanoscale layered materials synthesized by the method of modulated elemental reactants. J. Appl. Phys. 104, 033533 (2008).

    ADS  Google Scholar 

  15. 15.

    Chiritescu, C. et al. Ultralow thermal conductivity in disordered, layered WSe2 crystals. Science 315, 351–353 (2007).

    ADS  CAS  Google Scholar 

  16. 16.

    Hadland, E. C. et al. Ultralow thermal conductivity of turbostratically disordered MoSe2 ultra-thin films and implications for heterostructures. Nanotechnology 30, 285401 (2019).

    CAS  Google Scholar 

  17. 17.

    Hadland, E. et al. Synthesis, characterization, and ultralow thermal conductivity of a lattice-mismatched SnSe2 (MoSe2)1.32 heterostructure. Chem. Mater. 31, 5699–5705 (2019).

    CAS  Google Scholar 

  18. 18.

    Vaziri, S. et al. Ultrahigh thermal isolation across heterogeneously layered two-dimensional materials. Sci. Adv. 5, eaax1325 (2019).

    ADS  CAS  PubMed  PubMed Central  Google Scholar 

  19. 19.

    Kang, K. et al. High-mobility three-atom-thick semiconducting films with wafer-scale homogeneity. Nature 520, 656–660 (2015).

    ADS  CAS  Google Scholar 

  20. 20.

    Kang, K. et al. Layer-by-layer assembly of two-dimensional materials into wafer-scale heterostructures. Nature 550, 229–233 (2017).

    ADS  PubMed  PubMed Central  Google Scholar 

  21. 21.

    Liu, K. et al. Evolution of interlayer coupling in twisted molybdenum disulfide bilayers. Nat. Commun. 5, 4966 (2014).

    ADS  CAS  PubMed  PubMed Central  Google Scholar 

  22. 22.

    Sood, A. et al. Quasi-ballistic thermal transport across MoS2 thin films. Nano Lett. 19, 2434–2442 (2019).

    ADS  CAS  PubMed  PubMed Central  Google Scholar 

  23. 23.

    Li, Z. et al. Size dependence and ballistic limits of thermal transport in anisotropic layered two-dimensional materials. Preprint at arXiv https://arxiv.org/abs/1711.02772 (2017).

  24. 24.

    Schauble, K. et al. Uncovering the effects of metal contacts on monolayer MoS2. ACS Nano 14, 14798–14808 (2020).

    PubMed  PubMed Central  Google Scholar 

  25. 25.

    Liu, Y. et al. Approaching the Schottky–Mott limit in van der Waals metal–semiconductor junctions. Nature 557, 696–700 (2018).

    ADS  CAS  PubMed  PubMed Central  Google Scholar 

  26. 26.

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

    ADS  Google Scholar 

  27. 27.

    Jiang, P., Qian, X., Gu, X. & Yang, R. Probing anisotropic thermal conductivity of transition metal dichalcogenides MX2 (M = Mo, W and X = S, Se) using time-domain thermoreflectance. Adv. Mater. 29, 1701068 (2017).

    Google Scholar 

  28. 28.

    Evans, D. J. Homogeneous NEMD algorithm for thermal conductivity—application of non-canonical linear response theory. Phys. Lett. A 91, 457–460 (1982).

    ADS  Google Scholar 

  29. 29.

    Fan, Z., Chen, W., Vierimaa, V. & Harju, A. Efficient molecular dynamics simulations with many-body potentials on graphics processing units. Comput. Phys. Commun. 218, 10–16 (2017).

    ADS  CAS  Google Scholar 

  30. 30.

    Xu, K. et al. Thermal transport in MoS2 from molecular dynamics using different empirical potentials. Phys. Rev. B 99, 054303 (2019).

    ADS  CAS  Google Scholar 

  31. 31.

    Li, D., Schleife, A., Cahill, D. G., Mitchson, G. & Johnson, D. C. Ultralow shear modulus of incommensurate [SnSe]n[MoSe2]n layers synthesized by the method of modulated elemental reactants. Phys. Rev. Mater. 3, 043607 (2019).

    CAS  Google Scholar 

  32. 32.

    Huang, S. et al. Probing the interlayer coupling of twisted bilayer MoS2 using photoluminescence spectroscopy. Nano Lett. 14, 5500–5508 (2014).

    ADS  CAS  Google Scholar 

  33. 33.

    Erhart, P., Hyldgaard, P. & Lindroth, D. O. Microscopic origin of thermal conductivity reduction in disordered van der Waals solids. Chem. Mater. 27, 5511–5518 (2015).

    CAS  Google Scholar 

  34. 34.

    Yan, R. et al. Thermal conductivity of monolayer molybdenum disulfide obtained from temperature-dependent Raman spectroscopy. ACS Nano 8, 986–993 (2014).

    CAS  Google Scholar 

  35. 35.

    Sahoo, S., Gaur, A. P. S., Ahmadi, M., Guinel, M. J.-F. & Katiyar, R. S. Temperature-dependent Raman studies and thermal conductivity of few-layer MoS2. J. Phys. Chem. C 117, 9042–9047 (2013).

    CAS  Google Scholar 

  36. 36.

    Zhang, X. et al. Measurement of lateral and interfacial thermal conductivity of single- and bilayer MoS2 and MoSe2 using refined optothermal Raman technique. ACS Appl. Mater. Interfaces 7, 25923–25929 (2015).

    CAS  Google Scholar 

  37. 37.

    Taube, A., Judek, J., Łapińska, A. & Zdrojek, M. Temperature-dependent thermal properties of supported MoS2 monolayers. ACS Appl. Mater. Interfaces 7, 5061–5065 (2015).

    CAS  Google Scholar 

  38. 38.

    Lindroth, D. O. & Erhart, P. Thermal transport in van der Waals solids from first-principles calculations. Phys. Rev. B 94, 115205 (2016).

    ADS  Google Scholar 

  39. 39.

    Jo, I., Pettes, M. T., Ou, E., Wu, W. & Shi, L. Basal-plane thermal conductivity of few-layer molybdenum disulfide. Appl. Phys. Lett. 104, 201902 (2014).

    ADS  Google Scholar 

  40. 40.

    Wei, X. et al. Phonon thermal conductivity of monolayer MoS2: a comparison with single layer graphene. Appl. Phys. Lett. 105, 103902 (2014).

    ADS  Google Scholar 

  41. 41.

    Cai, Y., Lan, J., Zhang, G. & Zhang, Y.-W. Lattice vibrational modes and phonon thermal conductivity of monolayer MoS2. Phys. Rev. B 89, 035438 (2014).

    ADS  Google Scholar 

  42. 42.

    Liu, X., Zhang, G., Pei, Q.-X. & Zhang, Y.-W. Phonon thermal conductivity of monolayer MoS2 sheet and nanoribbons. Appl. Phys. Lett. 103, 133113 (2013).

    ADS  Google Scholar 

  43. 43.

    Ho, C. Y., Powell, R. W. & Liley, P. E. Thermal conductivity of the elements. J. Phys. Chem. Ref. Data 1, 279–421 (1972).

    ADS  CAS  Google Scholar 

  44. 44.

    Mavrokefalos, A., Nguyen, N. T., Pettes, M. T., Johnson, D. C. & Shi, L. In-plane thermal conductivity of disordered layered WSe2 and (W)x(WSe2)y superlattice films. Appl. Phys. Lett. 91, 171912 (2007).

    ADS  Google Scholar 

  45. 45.

    Ho, P. & Kwok, T. Electromigration in metals. Rep. Prog. Phys. 52, 301–348 (1989).

    ADS  CAS  Google Scholar 

  46. 46.

    Cao, Y. et al. Unconventional superconductivity in magic-angle graphene superlattices. Nature 556, 43–50 (2018).

    ADS  CAS  Google Scholar 

  47. 47.

    Sood, A. et al. An electrochemical thermal transistor. Nat. Commun. 9, 4510 (2018).

    ADS  PubMed  PubMed Central  Google Scholar 

  48. 48.

    Chen, S., Sood, A., Pop, E., Goodson, K. E. & Donadio, D. Strongly tunable anisotropic thermal transport in MoS2 by strain and lithium intercalation: first-principles calculations. 2D Mater. 6, 025033 (2019).

    CAS  Google Scholar 

  49. 49.

    Cueva, P., Hovden, R., Mundy, J. A., Xin, H. L. & Muller, D. A. Data processing for atomic resolution electron energy loss spectroscopy. Microsc. Microanal. 18, 667–675 (2012).

    ADS  CAS  Google Scholar 

  50. 50.

    Cahill, D. G. Analysis of heat flow in layered structures for time-domain thermoreflectance. Rev. Sci. Instrum. 75, 5119–5122 (2004).

    ADS  CAS  Google Scholar 

  51. 51.

    Liu, J. et al. Simultaneous measurement of thermal conductivity and heat capacity of bulk and thin film materials using frequency-dependent transient thermoreflectance method. Rev. Sci. Instrum. 84, 034902 (2013).

    ADS  Google Scholar 

  52. 52.

    Ditmars, D. A., Plint, C. A. & Shukla, R. C. Aluminum. I. Measurement of the relative enthalpy from 273 to 929 K and derivation of thermodynamic functions for Al(s) from 0 K to its melting point. Int. J. Thermophys. 6, 499–515 (1985).

    ADS  CAS  Google Scholar 

  53. 53.

    Yue, X. F. et al. Measurement of interfacial thermal conductance of few-layer MoS2 supported on different substrates using Raman spectroscopy. J. Appl. Phys. 127, 104301 (2020).

    ADS  CAS  Google Scholar 

  54. 54.

    Lee, J.-U., Yoon, D., Kim, H., Lee, S. W. & Cheong, H. Thermal conductivity of suspended pristine graphene measured by Raman spectroscopy. Phys. Rev. B 83, 081419 (2011).

    ADS  Google Scholar 

  55. 55.

    Cai, W. et al. Thermal transport in suspended and supported monolayer graphene grown by chemical vapor deposition. Nano Lett. 10, 1645–1651 (2010).

    ADS  CAS  Google Scholar 

  56. 56.

    Wang, Y., Cong, C., Qiu, C. & Yu, T. Raman spectroscopy study of lattice vibration and crystallographic orientation of monolayer MoS2 under uniaxial strain. Small 9, 2857–2861 (2013).

    CAS  Google Scholar 

  57. 57.

    Chen, S. et al. Raman measurements of thermal transport in suspended monolayer graphene of variable sizes in vacuum and gaseous environments. ACS Nano 5, 321–328 (2011).

    Google Scholar 

  58. 58.

    Hsu, I.-K., Pettes, M. T., Aykol, M., Shi, L. & Cronin, S. B. The effect of gas environment on electrical heating in suspended carbon nanotubes. J. Appl. Phys. 108, 084307 (2010).

    ADS  Google Scholar 

  59. 59.

    Hermann, K. Periodic overlayers and moiré patterns: theoretical studies of geometric properties. J. Phys. Condens. Matter 24, 314210 (2012).

    ADS  Google Scholar 

  60. 60.

    Hjorth Larsen, A. et al. The atomic simulation environment—a Python library for working with atoms. J. Phys. Condens. Matter 29, 273002 (2017).

    PubMed  PubMed Central  Google Scholar 

  61. 61.

    Liang, T., Phillpot, S. R. & Sinnott, S. B. Parametrization of a reactive many-body potential for Mo–S systems. Phys. Rev. B 79, 245110 (2009).

    ADS  Google Scholar 

  62. 62.

    Plimpton, S. Fast parallel algorithms for short-range molecular dynamics. J. Comput. Phys. 117, 1–19 (1995).

    ADS  CAS  MATH  Google Scholar 

  63. 63.

    Evans, D. J. Homogeneous NEMD algorithm for thermal conductivity—application of non-canonical linear response theory. Phys. Lett. A 91, 457–460 (1982).

    ADS  Google Scholar 

  64. 64.

    Togo, A. & Tanaka, I. First principles phonon calculations in materials science. Scr. Mater. 108, 1–5 (2015).

    ADS  CAS  Google Scholar 

  65. 65.

    Togo, A., Chaput, L. & Tanaka, I. Distributions of phonon lifetimes in Brillouin zones. Phys. Rev. B 91, 094306 (2015).

    ADS  Google Scholar 

  66. 66.

    Chaput, L. Direct solution to the linearized phonon Boltzmann equation. Phys. Rev. Lett. 110, 265506 (2013).

    ADS  PubMed  PubMed Central  Google Scholar 

  67. 67.

    Fransson, E., Slabanja, M., Erhart, P. & Wahnström, G. dynasor—a tool for extracting dynamical structure factors and current correlation functions from molecular dynamics simulations. Adv. Theory Simul. 4, 2000240 (2021).

    CAS  Google Scholar 

  68. 68.

    Yalon, E. et al. Energy dissipation in monolayer MoS2 electronics. Nano Lett. 17, 3429–3433 (2017).

    ADS  CAS  PubMed  PubMed Central  Google Scholar 

  69. 69.

    Freedy, K. M., Olson, D. H., Hopkins, P. E. & McDonnell, S. J. Titanium contacts to MoS2 with interfacial oxide: interface chemistry and thermal transport. Phys. Rev. Mater. 3, 104001 (2019).

    CAS  Google Scholar 

  70. 70.

    Lombard, J., Detcheverry, F. & Merabia, S. Influence of the electron–phonon interfacial conductance on the thermal transport at metal/dielectric interfaces. J. Phys. Condens. Matter 27, 015007 (2015).

    ADS  CAS  Google Scholar 

  71. 71.

    Kimling, J., Philippi-Kobs, A., Jacobsohn, J., Oepen, H. P. & Cahill, D. G. Thermal conductance of interfaces with amorphous SiO2 measured by time-resolved magneto-optic Kerr-effect thermometry. Phys. Rev. B 95, 184305 (2017).

    ADS  Google Scholar 

  72. 72.

    Zhao, Y. et al. Interlayer breathing and shear modes in few-trilayer MoS2 and WSe2. Nano Lett. 13, 1007–1015 (2013).

    ADS  CAS  Google Scholar 

  73. 73.

    Huang, S. et al. Low-frequency interlayer Raman modes to probe interface of twisted bilayer MoS2. Nano Lett. 16, 1435–1444 (2016).

    ADS  CAS  Google Scholar 

  74. 74.

    Zhang, X. et al. Raman spectroscopy of shear and layer breathing modes in multilayer MoS2. Phys. Rev. B 87, 115413 (2013).

    ADS  Google Scholar 

Download references

Acknowledgements

We thank D. F. Ogletree, J. Jureller, A. J. Mannix, J.-U. Lee, K.-H. Lee and M. Lee for their helpful discussions. We also acknowledge J.-H. Kang, A. Ye and C. Liang for their help with materials preparation. Primary funding for this work comes from Air Force Office of Scientific Research MURI projects (FA9550–18–1-0480 and FA9550–16-1-0031). Material growths done by F.M. and C.P. are partially supported by the University of Chicago MRSEC (NSF DMR-2011854) and Samsung Advanced Institute of Technology. This work makes use of the characterization facilities of the University of Chicago MRSEC (NSF DMR-2011854) and the Pritzker Nanofabrication Facility at the University of Chicago, which receives support from SHyNE Resource (NSF ECCS-1542205), a node of NSF’s NNCI network. TDTR measurements are supported by Office of Naval Research MURI grant N00014-16-1-2436 and are carried out in Frederick Seitz Materials Research Laboratory at the University of Illinois at Urbana-Champaign. F.E., E.F. and P.E. are funded by the Knut and Alice Wallenberg Foundation (2014.0226), the Swedish Research Council (2015-04153 and 2018-06482), and the FLAG-ERA JTC-2017 project MECHANIC funded by the Swedish Research Council (VR 2017-06819); they acknowledge the computer time allocations by the Swedish National Infrastructure for Computing at NSC (Linkӧping) and C3SE (Gothenburg). A.R. and electron microscopy at the Cornell Center for Materials Research are supported by the NSF MRSEC grant DMR-1719875. The Titan microscope was acquired with the NSF MRI grant DMR-1429155. F.M. acknowledges support by the NSF Graduate Research Fellowship Program under grant no. DGE-1746045. Y.Z. acknowledges support by the Camille and Henry Dreyfus Foundation, Inc., under the Dreyfus Environmental Postdoc award EP-16-094.

Author information

Affiliations

Authors

Contributions

S.E.K., J.S., and J.P. conceived the experiments. S.E.K. produced and performed the structural characterizations of the r-TMD films, and ran the finite-element simulations. F.M. and C.P. grew the TMD monolayers. A. Rai and D.G.C. performed the TDTR measurements. F.E., E.F. and P.E. conducted the atomistic simulations. P.P. fabricated the Au electrodes for the electromigration experiments. A. Ray and D.A.M. performed the STEM imaging. S.E.K. performed the Raman measurements and analysed the data with the help of D.G.C. and Y.Z. S.E.K., F.M. and J.P. wrote the manuscript. J.P., D.G.C. and P.E. oversaw the project. All authors discussed the manuscript and provided feedback.

Corresponding authors

Correspondence to Paul Erhart, David G. Cahill or Jiwoong Park.

Ethics declarations

Competing interests

The authors declare no competing interests.

Additional information

Peer review information Nature thanks Davide Donadio, M. Muñoz Rojo and Mariusz Zdrojek for their contribution to the peer review of this work.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Extended data figures and tables

Extended Data Fig. 1 Cross-sectional TEM images of N = 20 and 10 r-MoS2 on AlOx.

Each set of N images are taken from the same sample at different locations. The N = 10 films are coated with Al that was electron beam evaporated onto the surface.

Extended Data Fig. 2 GIWAXS data of N = 10 r-MoS2.

The peak position corresponds to a 2θ value of 14˚, which translates to an interlayer spacing of 6.4 Å (scattering direction).

Extended Data Fig. 3 Additional TDTR measurements and details.

a, Microscope image of an N = 10 r-MoS2 film coated with a square grid of Al pads. b, 4 × 4 TDTR map of RTDTR of an N = 10 r-MoS2 film. c, Histogram of RTDTR array measurements. d, TDTR measurements of N ≤ 10 r-TMD films coated with Au or Al. The error bars denote s.d.; number of TDTR measurements per film sample, n = 3 for Au samples; n = 3–5 for the Al samples. e, Picosecond acoustics of a MoS2 monolayer on thick sapphire substrate, coated with an Al transducer layer. The y axis Vin is the in-phase signal of the lock-in amplifier. The red arrows indicate the acoustic waves reflected at the Al/MoS2 interface.

Extended Data Fig. 4 Low-frequency Raman modes of r-MoS2.

a, Raman spectra reflecting the breathing modes (BM) of r-MoS2 (blue) and the shear mode (SM) for MoS2. b, The low-frequency Raman peak positions of r-MoS2 and exfoliated MoS2. The filled squares indicate the BM peak positions of r-MoS2. The open squares indicate the BM peak positions of exfoliated MoS2, and the open circles indicate the SM peak positions of exfoliated MoS272.

Extended Data Fig. 5 Raman thermometry on r-MoS2 films.

a, Δω–Pabs curves of representative N = 2 r-MoS2 films at different pressures. The Pabs values along the x axis are normalized to account for the slight differences in beam spot sizes (Δr = 20%). The results at 15 torr and 4 mtorr signify no effect of reducing the pressure to below 15 torr. Δω was approximately five-fold smaller at atmospheric pressure due to the extra heat loss channel by air. b, Δω–Pabs curves of r-MoS2 films made up of D = 400 nm (grain size) monolayers. c, Optical absorption of suspended r-MoS2 films, which follows the trend \(A=1-{(1-{A}_{0})}^{N}\), where A0 is the monolayer absorptance. From the fit, A0 = 0.08 ± 0.003. d, ω–T calibration measurements of suspended r-MoS2 films (D = 1 μm), with the N = 2 and N = 4 data as the representative curves. e, ω–T slopes versus layer number for all films, with D = 400 nm or 1 μm.

Extended Data Fig. 6 κ(T) and ρ of r-MoS2.

a, κ(T) of r-MoS2, with κ|| measured using Raman thermometry of N = 4 r-MoS2, and κ measured via TDTR as reported in Fig. 2c. The error bars are the propagated uncertainties from the calculation of the conductance value for each N. The error bars denote the propagated uncertainty of the calculations from the input parameters. We observed that κ|| decreased with T, alluding to phonon-mediated heat transport and attesting to the long-range crystallinity of the r-MoS2 films in-plane. This was in contrast with the κ(T) behaviour (Fig. 2c), which showed a slightly increasing trend. b, Catalogue of experimentally measured anisotropy ratios at room temperature versus slow-axis thermal conductivity (κs) of thermally anisotropic materials from literature, by category.

Extended Data Fig. 7 r-MoS2 efficacy as a heat spreader.

a, Finite element simulations of the linear temperature profiles of Au electrodes covered with MoS2 and r-MoS2. b, SiNx as heat spreaders for Au electrodes. Electrical properties of 10-nm-thick, 100-nm-wide and 10-μm-long Au electrodes before and after 16 nm SiNx film deposition onto the electrodes using plasma-enhanced chemical vapour deposition. In contrast to r-MoS2, the direct deposition of an ultrathin inorganic film such as SiNx with a comparable κ to κ|| of r-MoS2 negatively affects the performance of the Au electrodes.

Extended Data Fig. 8 Optimization of the MD simulations for κ calculations.

a, Optimization of the driving force of the system, where the grey zone denotes the error. b, Effect of thermal expansion on κ.

Extended Data Table 1 Comparison of literature thermal boundary resistances with experimentally measured values
Extended Data Table 2 Parameters used in the calculations of the thermal conductivity using GPUMD

Source data

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Kim, S.E., Mujid, F., Rai, A. et al. Extremely anisotropic van der Waals thermal conductors. Nature 597, 660–665 (2021). https://doi.org/10.1038/s41586-021-03867-8

Download citation

Comments

By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.

Search

Quick links