Extremely anisotropic van der Waals thermal conductors

The densification of integrated circuits requires thermal management strategies and high thermal conductivity materials1–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.

The densification of integrated circuits requires thermal management strategies and high thermal conductivity materials [1][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 axis 4,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 MoS 2 , 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 MoS 2 (57 ± 3 mW m −1 K −1 ) and WS 2 (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 MoS 2 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.
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 nanostructuring 6 , such as fabricating inorganic superlattices [7][8][9][10][11] or designing symmetry-breaking crystal architectures in a single material 12 . 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, graphite 1 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, WSe 2 ) 14-17 and heterostructures 18 . 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 methods 19,20 (see Methods).
The transmission electron microscopy (TEM) diffraction (Fig. 1c, left) and darkfield (inset) images from a representative MoS 2 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 interfaces 20 with interlayer rotation at every stacked interface. The TEM diffraction pattern of N = 10 r-MoS 2 (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-MoS 2 ( 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 MoS 2 multilayers 21 . Both the growth and stacking steps are scalable, as shown by the optical images of N = 1 and N = 10 r-MoS 2 films (~1 cm 2 ) 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-MoS 2 or r-WS 2 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 (−V in /V out 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-MoS 2 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 R TDTR , the total thermal resistance between the Al transducer layer and sapphire across the r-TMD film for different N. Figure 2b shows R TDTR versus N for r-MoS 2 (N ≤ 22; solid circles) and r-WS 2 (N ≤ 10; open circles) measured under ambient conditions. We make two observations. First, R TDTR monotonically increases with N. Second, R TDTR 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 MoS 2 (as thick as 240 nm) 22,23 . A single parameter κ ⊥ characterizes the thermal resistance across r-MoS 2 (or r-WS 2 ) using the equation R TDTR = R 0 + Nd/κ ⊥ , where Nd is the total film thickness, and R 0 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-MoS 2 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 deposition 24,25 . We measure κ ⊥ = 57 ± 3 mW m −1 K −1 for r-MoS 2 and κ ⊥ = 41 ± 3 mW m −1 K −1 for r-WS 2 , which are similar to the lowest value ever observed in a fully dense solid 15 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 MoS 2 (2-5 W m −1 K −1 ) 26,27 or WS 2 (~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-MoS 2 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 MoS 2 (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-MoS 2 and bulk MoS 2 (see Methods and Extended Data Table 2) 28-30 . Figure 2c shows κ || and κ ⊥ of r-MoS 2 (solid circles) and bulk MoS 2 (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 MoS 2 to 0.16 ± 0.04 W m −1 K −1 in r-MoS 2 at 300 K, and also does not decrease with T, suggesting a transition away from the phonon-limited thermal transport mechanism observed in bulk MoS 2 .

Article
Further analysis of the vibrational spectrum of r-MoS 2 enables us to break down the reduction in κ ⊥ in terms of the changes in the group velocities (v g ) and lifetimes (τ), which are the two factors that determine the thermal conductivity according to Boltzmann transport theory. Figure 2d shows that the v g of the through-plane longitudinal acoustic (LA) mode in r-MoS 2 remains similar to that of bulk MoS 2 (dashed lines), but the transverse acoustic (TA) modes in r-MoS 2 undergo extreme softening with their v g s practically vanishing 31 . This implies a loss of resistance with respect to lateral shear, consistent with the low-frequency Raman spectra of r-MoS 2 films (see Methods and Extended Data Fig. 4) and previous calculations 32, 33 . In addition, the τ of both the LA and the TA modes ( Fig. 2e) in r-MoS 2 are more than one order of magnitude smaller than in bulk MoS 2 , with the LA lifetimes being close to the period of the LA mode vibration (dashed line). From these results, the median mean free path ∼ l = v g τ 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 >> ∼ 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-MoS 2 . 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-MoS 2 film ( Fig. 3a; hole diameter of 5 μm, at 15 torr), which increases its temperature (ΔT) locally upon absorbing laser power P abs . Δ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. Figure 3c plots Δω versus P abs for r-MoS 2 with different N (2 to 5). The slope of the linear fit (|d(Δω)/dP abs |), 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-MoS 2 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 ( 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 reports 23,[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-MoS 2 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 phononbased solids 15,27,43,44 (for a full comparison, see Extended Data Fig. 6b). Compared to a bulk MoS 2 crystal (ρ ≈ 20) 27 or disordered layered WSe 2 (ρ ≈ 30) 15,44 , r-MoS 2 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-MoS 2 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-nmthick r-MoS 2 film draped over a nanoscale Au electrode (15 nm tall, 100 nm wide) on a 50 nm SiO 2 /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-MoS 2 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-MoS 2 films provides thermal insulation in the through-plane direction, with much lower MoS 2 surface ∆T values that are only one-third of the value of the bare Au electrode. While single-crystal MoS 2 displays similar properties, the insulation effect is stronger in r-MoS 2 (Extended Data Fig. 7a). This implies that heat is efficiently directed away from the hot Au electrode laterally through r-MoS 2 but not to the surface of r-MoS 2 , making the surface of the entire device significantly cooler.
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-MoS 2 (~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 fail 45 . 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-MoS 2 can carry a larger current without breaking. The histogram of critical current I c (maximum current a Au electrode sustains for at least 20 s) measured from 20 electrodes (10 bare and 10 with r-MoS 2 ) reveals a ~50% increase in the median I c values (Fig. 4e). These results demonstrate our r-MoS 2 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 I c and maximum power before breaking is in good agreement with our simulation in Fig. 4c. Furthermore, we note that the r-MoS 2 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 Article 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 graphene 46 . Interlayer rotations can be combined with other parameters (such as pressure or interlayer spacing 47,48 ) and advanced structures (superlattices and heterostructures 18 ) to realize highly tunable ρ, allowing for the customization of thermal transport properties with an unprecedented level of directional and spatial control.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-021-03867-8. Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
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/.

Sample preparation
Large-area, polycrystalline transition metal dichalcogenide (TMD) (MoS 2 and WS 2 ) monolayers were grown on SiO 2 /Si substrates in a hot-walled tube furnace via metal-organic chemical vapour deposition adapted from a previously reported protocol 19 . 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 N 2 to 15 torr) were used as the metal precursors for the MoS 2 and WS 2 growths, respectively. (C 2 H 5 ) 2 S was used as the chalcogen source. All precursors were kept at room temperature. N 2 and H 2 were used as carrier gases. Typical growth times were 15-20 h for MoS 2 at a growth temperature of 525 °C. Typical growth times for WS 2 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 method 20 .
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/H 2 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 SiN x transmission electron microscopy (TEM) grids were cleaned in a N 2 /H 2 plasma at 100 °C and 180 mtorr for 3 min, followed by the transfer of stacked MoS 2 films onto the TEM grids. During the transfer process, the PMMA-coated r-MoS 2 was suspended on holey thermal release tape before contacting the TEM grid. The extra PMMA-MoS 2 not on the TEM grid was cut away at 180 °C so the PMMA layer was softened. PMMA was removed from r-MoS 2 on the TEM grid via annealing the film in 400/100 SCCM Ar/H 2 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-MoS 2 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 Imager 49 . 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 MoS 2 (bright bands), followed by an AlO x 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/e 2 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 outof-phase signals from the lock-in was fitted to a thermal diffusion model. The full details of the TDTR measurement can be found elsewhere 50,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 literature 52 . 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-MoS 2 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 G 1 . We used the bulk value of the volumetric heat capacity of 1.89 J K −1 cm −3 for the r-MoS 2 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 G 1 . To obtain κ TMD from G 1 , we perform TDTR on various N-layer TMD films, then perform a linear fit on the effective thermal resistance (R TDTR , equal to 1/G 1 ) 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 (R 0 ) of the top and bottom interfaces. In Extended Data Table 1, our R 0 values match the values reported in literature 22,27,53 . We note that, although R 0 changes depending on the chemical nature of the metal-TMD interface, the slope of the R TDTR -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/e 2 beam radius of ~3.2 μm, and so we deliberately chose a larger f and a 1/e 2 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 reports 34, 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-MoS 2 A 1g 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 A 1g peak since this out-of-plane vibrational mode is less sensitive to in-plane strain 56 . The ω-T calibration measurements were performed at atmospheric pressure and with low laser powers. The stage was purged with dry N 2 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-MoS 2 films, where a linear fit was performed to obtain the temperature-dependent Raman coefficients. This process was repeated for r-MoS 2 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-MoS 2 films with N = 2-5 and D = 1 μm in Fig. 3c and for r-MoS 2 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 < 0 ω P d d 2 2 . 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 (r 0 ) 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 r 0 = 0.71 ± 0.09 μm. The laser powers were measured using a Thorlabs standard silicon photodiode power sensor. The r-MoS 2 absorbance A = Absorbed light intensity 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-MoS 2 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 reports 20 . We use the value A measured at room temperature for our Raman analysis.
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 r R , ≤ susp and T r R , ≥ supp , T 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-MoS 2 and SiN x .
We solved for the expression of T(r) and obtained an expression for the average temperature measured by the Raman shift κ || was obtained by substituting the experimentally measured value for T 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 where ω is the Raman frequency of A 1g 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 κ || . Variable pressure Raman thermometry measurements. Previous Raman thermometry measurements on graphene films 57 and carbon nanotubes 58 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-MoS 2 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-MoS 2 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 T a 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 T a . We plot a κ || versus T curve, where the x axis is T = T a (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-MoS 2 heat spreader experiments (electromigration of Au nanoelectrodes)
All Au electrodes were fabricated on Si substrates with 50 nm dry SiO 2 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 O 2 plasma for 30 s to remove any resist residue and to promote adhesion of the r-MoS 2 film to the Au electrodes and the SiO 2 surface. A PMMA coated N = 16 r-MoS 2 film was transferred onto the electrodes using the same process as the stacking method as outlined above. The PMMA on the r-MoS 2 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 I c 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 SiN x 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 SiH 4 and N 2 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 literature 59 , which was implemented in Python using the atomic simulation environment package 60 . The structure models were subsequently relaxed using an analytic bond-order potential 61 and implemented in the LAMMPS package 62 .
The main r-MoS 2 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-MoS 2 and yields thermal conductivities and phonon dispersions of bulk MoS 2 that agree with previous experimental observations and Boltzmann transport calculations based on density functional theory 38 , 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 simulations 63 and implemented in the graphics processing units molecular dynamics (GPUMD) package 29 . We also included the effects of thermal expansion in the simulations (Extended Data Fig. 8b). The calculated κ ⟂ values of r-MoS 2 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 MoS 2 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 package 64 . 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 package 65 , which also provided us with the thermal conductivity as obtained from a direct solution of the Boltzmann transport equation 66 . 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 MoS 2 to that calculated using MD simulations to verify the accuracy of our MD simulations for calculating the phonon dispersion of r-MoS 2 . For both bulk and r-MoS 2 , 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 dynasor 67 . 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 paper 67 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 SiO 2 / 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 SiO 2 layer on a Si substrate. We layer a 10-nm-thick MoS 2 film onto the Au electrode and the SiO 2 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-MoS 2 /Au, r-MoS 2 /SiO 2 , Au/SiO 2 , SiO 2 /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-MoS 2 films, along with the spectrum for MoS 2 , are shown in Extended Data Fig. 4a. From the layer-dependence of the peak positions, we assigned these to be