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.
figure 1

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.

Full size image

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.
figure 2

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.

Full size image

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.
figure 3

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.

Full size image

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.
figure 4

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.

Full size image

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 RTDTR–N 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.