Experimental details
MoSe2 and MoS2 flakes of varying thicknesses were mechanically exfoliated from bulk crystals onto a viscoelastic polydimethylsiloxane stamp and subsequently transferred using a dry-transfer technique onto titanium (5 nm)/gold (50 nm)-coated silicon nitride substrates containing circular holes with a diameter of 15 μm (Norcada, NTPR005D-C15), resulting in suspended flakes. The flakes are monocrystalline and exhibit clean surfaces and low strain36. For the measurements, we used samples suspended over the 15-μm holes, with similar results also obtained for the 10-μm holes (Supplementary Fig. 2). The spatiotemporal pump–probe thermometry setup (Extended Data Fig. 2 and ref. 21) uses a FLINT laser from LIGHT CONVERSION, operating at a frequency of 76 MHz, which produces pulses centred at 1,030 nm with <200-fs temporal resolution. The majority of the laser power drives an optical parametric oscillator that generates a variable signal output ranging from 1,320 to 2,000 nm. We use the second or third harmonic of this signal output as probe, tuning it to the exciton resonances—800 nm for MoSe2 and 615 nm for MoS2. Utilizing a scanning mirror system (Optics in Motion OIM101), we control the probe’s direction, whereas the pump beam (at 515 nm, the second harmonic of the fundamental laser source) goes through a Newport DL255 delay line and an electro-optic modulator that controls the pump modulation frequency fmod (240 kHz). Both beams are combined using a dichroic mirror and then focused to submicrometre spots via a microscope objective lens with a numerical aperture of 0.65. The pump and probe spot sizes are σpu ≈ 0.27 μm and σpr ≈ 0.35 μm, respectively. We detect the reflected probe beam using a silicon photodiode.
The spatiotemporal pump–probe thermometry measurements, described in detail in ref. 21, rely on the fact that within a nanosecond37, energy from the initial pump-induced electronic excitation transfers to the lattice, where phonons transport the energy as lattice heat towards the metallic heat sink at the edge of the suspended crystal, which takes place on a microsecond timescale (Extended Data Fig. 3). In the case of purely diffusive heat transport following Fourier’s law, this leads to a spatial temperature profile that exclusively depends on thermal diffusivity, the accurately known geometry of the system, and two frequencies related to the pulsed photoexcitation: the laser repetition rate frep and the pump modulation frequency fmod. Measuring the spatial profile, thus, provides direct access to thermal diffusivity. We record the spatial profile at a slightly negative pump–probe delay, corresponding to a delay time of 1/frep ≈ 13 ns, where the transient reflectivity results from the pump-induced lattice temperature change, without any contribution from the photoexcited electronic species. The transient reflectivity is sensitive to the lattice temperature ΔT via the temperature-dependent exciton linewidth38. We modulate the intensity of the pump beam at a frequency fmod = 240 kHz and demodulate using a lock-in amplifier (Zurich MFLI), obtaining in-phase and out-of-phase signals. We perform all measurements in a vacuum (~10−5 mbar) and at room temperature.
To ensure that we can accurately correlate the transient reflectivity signal to a lattice temperature, we need to operate in the linear regime, where an increase in pump fluence gives a proportional increase in temperature and, therefore, transient reflectivity. To confirm that we are in the linear regime in our experiments, we first conduct temporal dynamics measurements by varying the pump fluence. In addition, we perform spatial scans on the suspended region by varying the pump fluence. Extended Data Fig. 6a–d shows the results obtained for an almost eight-layer MoSe2 sample. Both in-phase and out-phase signals exhibit a linear increase with pump fluence, as evident from their respective normalized plots. This linear dependence of the signal supports the use of the transient reflectivity profile as a representation of the spatial temperature profile. Moreover, we vary the probe fluence (and fix the pump power) on the suspended trilayer MoSe2 sample, yielding similar observations (Extended Data Fig. 6e–h). On the basis of these measurements, we keep the pump and probe fluences in the range of 6–20 μJ cm−2, to ensure that the system is in a regime in which the transient reflectivity scales linearly with pump power (Extended Data Figs. 5, 6 and 10). In this case, the typical temperature increase is on the order of 10 K.
Given recent reports that layered semiconductors can exhibit exotic electron–hole plasma39 phases due to strong electron interactions, we conducted photoluminescence (PL) measurements on a monolayer MoSe2 sample to confirm that we are operating within the linear electronic regime. This ensures that our observations are not influenced by nonlinear electron–hole plasma effects. Extended Data Fig. 10a displays the PL spectra with increasing laser fluence (pump, 515 nm). Consistently, we observe a linear increase in spectral average (area under PL spectrum) with laser fluence (Extended Data Fig. 10b), suggesting that we are well below the Mott transition limit and ruling out the presence of such exotic phases. In particular, we conduct these PL measurements under almost identical configurations as the spatiotemporal thermometry measurements, with the only difference being the collection of PL emission by a fibre cable to the spectrometer (Extended Data Fig. 10c).
Phenomenological diffusive Fourier model
To obtain an effective diffusivity from the transient reflectivity evolution observed in our experiments, we use Fourier’s diffusive heat law to model the propagation of heat from the localized pump-induced hotspot to the heat sink:
$$\frac{\partial T}{\partial t}=D{\nabla }^{2}T+\frac{Q}{{C}_{{\rm{v}}}},$$
(3)
where T is the lattice temperature, t is time, D is the diffusivity, Q is the power injected by the pump laser and Cv is the heat capacity. Initially, we start with a Gaussian temperature profile that represents the localized heating from a single pump pulse. As time progresses, this temperature profile diffuses outwards and decays, as heat flows towards the heat sink. Due to the relatively short time interval between consecutive laser pulses (13 ns, corresponding to the 76-MHz repetition rate, frep), the temperature profile does not fully return to equilibrium before the next pump pulse arrives.
At t = 13 ns (or \(t=\frac{1}{{f}_{{\rm{rep}}}}\)), we inject a Gaussian energy pulse into the existing, partially decayed, temperature distribution. This cycle repeats until the end of the pump-on window. The number of these Gaussian pulses (N) that contribute to the accumulated heat during the pump-on window is determined by the laser repetition rate and the pump modulation frequency (fmod), given by \(N=\frac{{f}_{{\rm{rep}}}}{2{f}_{{\rm{mod}}}}\) for a 50% duty cycle (~160 pulses in our case). When the pump is blocked during the pump-off window, the accumulated heat fully dissipates. The heating cycle then resumes with the next pump-on window.
We perform our simulations in two dimensions with radial symmetry and calculate the temperature rise ΔT(x, y) over time using the forward-time centred-space finite difference method. We set the boundary condition at r = 7.5 μm (the edge of the suspended region) as a perfect heat sink, ΔT = 0 K, to represent the gold-coated substrate. To obtain the in-phase and out-of-phase convoluted transient reflectivity profiles from the model predictions, we perform the lock-in operations described in the ‘Describing in-phase and out-of-phase signals’ section.
DFT input for the mesoscopic model
The ab initio calculations are conducted using density functional theory as implemented in the SIESTA program40,41, as described in detail in ref. 26. Briefly, the cell and atomic positions are optimized using a conjugate gradient method, ensuring that the maximum forces on atoms are kept below 10−5 eV Å−1. For low-dimensional structures, a vacuum thickness of 17 Å is used to avoid periodic interactions along the stacking direction. Unlike in ref. 24, we do not normalize the thermal conductivity taking into account this additional volume. All calculations use the Perdew–Burke–Ernzerhof generalized gradient approximation42,43 with van der Waals functionals, as reparameterized in ref. 44. A k-mesh of 20 × 20 × 1 is used for monolayers, bilayers and trilayers, and 20 × 20 × 20 for other systems, with norm-conserving pseudopotentials45 sourced from the PseudoDojo library46.
To compute the interatomic force constants (IFCs), we generate supercells of size 10 × 10 × 10 for bulk systems and 10 × 10 × 1 for the other structures. Atomic displacements are thermally initialized using a Bose–Einstein distribution, simulating phonons at 300 K, using TDEP routines47. With these thermally populated supercells, we calculate both second- and third-order IFCs, with real-space force cut-offs of 5 Å for the third-order IFCs and 8 Å for the second-order IFCs. Thermal properties calculations are subsequently conducted on a q-point grid of 128 × 128 × 1 for monolayers, 64 × 64 × 1 for bilayers and trilayers, and 24 × 24 × 24 for bulk systems. The thermal properties calculated for MoSe2 and MoS2 are presented in Tables 1 and 2, respectively, and serve as input parameters for our mesoscopic model. Specifically, we directly obtain the heat capacity Cv, and use the full solution of the Boltzmann transport equation to obtain the thermal conductivity κ. Furthermore, the mode-resolved lifetimes τλ, group velocities vλ and phase velocities cλ serve as input parameters to obtain the heat flux relaxation time τ and the non-local length ℓ used in the hydrodynamic transport equation that we describe later28.
Hydro-thermoelastic model
The hydrodynamic and thermoelastic model predictions of the convoluted thermal profiles are obtained by solving the hydro-thermoelastic transport equation (equation (2)), along with the energy balance equation including thermomechanical energy exchange, and Newton’s law for the elastic field evolution including a linear thermal expansion term (Supplementary Note 1). Numerical solutions are obtained using finite elements29. In Tables 1 and 2, we present the intrinsic thermoelastic parameter values appearing in the model equations for MoSe2 and MoS2, respectively. In the absence of hydrodynamic and thermoelastic effects, the hydro-thermoelastic equation reduces to Fourier’s law. By combining this simplified transport equation with the energy balance equation, one recovers the thermal diffusion equation described in the ‘Phenomenological diffusive Fourier model’ section.
We impose radial symmetry and restrict the deformation to the in-plane directions. The temperature and displacement in the edges of the sample are fixed. The excitation is modelled as a Gaussian power density in the energy balance equation with a characteristic size of L = 700 nm, which is assumed to be larger than the optical spot size to account for fast diffusion of the optically excited hot electrons and excitons before they deposit energy into the lattice48. We note that the resulting temperature profiles are qualitatively insensitive to the exact size of the thermal source. In the time domain, the excitation power density follows a square function with a frequency of 240 kHz. Thus, we do not explicitly simulate the energy pulses during the pump-on stage, as described earlier. This simplification greatly reduces the computational cost by preventing the propagation of acoustic waves, and we verified that it does not significantly modify the average temperature and elastic fields predicted by the model (Supplementary Fig. 3). Finally, to obtain the in-phase and out-of-phase convoluted thermal profiles from the model predictions, we perform the lock-in operations described in the ‘Describing in-phase and out-of-phase signals’ section. More details are provided in Supplementary Note 1.
Describing in-phase and out-of-phase signals
To model the photoresponse, specifically the in-phase and out-of-phase signals, we use the results of the simulations from both the phenomenological diffusive Fourier model and the mesoscopic hydrodynamic and thermoelastic models, and simulate the lock-in detection process. As an example, Extended Data Fig. 4a shows the temperature dynamics obtained from the Fourier model over multiple modulation cycles. First, we examine the case with a thermal diffusivity of 0.20 cm2 s−1 at a modulation frequency fmod of 10 kHz. The temperature rises during the pump-on period, reaching its peak by the end of this window, and then decays during the pump-off period, returning to zero in approximately 2 μs during the pump-off window (Extended Data Fig. 4a).
To simulate the lock-in operation, we multiply this temperature signal by a sine function (to obtain the in-phase signal) and a cosine function (to obtain the out-of-phase signal), both at the modulation frequency of the pump (Extended Data Fig. 4b,c). We obtain the in-phase and out-of-phase signals by averaging these products over time. Repeating this process for temperature dynamics at different spatial positions (from the hotspot to the heat sink) allows us to reconstruct the in-phase and out-of-phase spatial profiles (Extended Data Fig. 4d). The cosine-weighted temperature dynamics average out to nearly zero, whereas the sine-weighted dynamics remain finite, indicating a predominantly in-phase response and a negligible out-of-phase contribution.
Next, we consider a higher modulation frequency, fmod = 240 kHz, and repeat the Fourier simulation with D = 0.20 cm2 s−1. Extended Data Fig. 4e shows the resulting temperature dynamics at various spatial positions. The temperature rises during the pump-on period, reaches its maximum near the end of this window, and then decays during the pump-off period, returning nearly to zero before the onset of the next pump-on cycle. Unlike the 10-kHz case, a finite-temperature signal persists during the pump-off window. This key difference leads to a finite out-of-phase contribution when averaging the cosine-weighted temperature dynamics (Extended Data Fig. 4f–h).
Finally, we consider the case of a lower diffusivity, 0.02 cm2 s−1, and repeat the Fourier simulation at fmod = 240 kHz. Extended Data Fig. 4i illustrates the temperature dynamics for this scenario at various spatial positions. Two key differences emerge. First, as shown in Extended Data Fig. 3, the temperature does not fully decay during the pump-off period, unlike in the case of higher diffusivity. Second, around 2 μm from the heat source, the temperature rises slowly during the pump-on window and decays equally slowly during the pump-off window. The lower diffusivity causes heat to propagate more slowly, resulting in delayed temperature rise, with the maximum temperature at 2 μm occurring during the pump-off phase.
This delayed heating causes the temperature dynamics at this location to be π out of phase, leading to a sign change in the in-phase signal. Note that the out-of-phase signal is π/2 out of phase with respect to the in-phase signal. Mathematically, these negative in-phase signals occur because when the temperature dynamics are multiplied by the sine wave, a substantial portion of the sine’s negative part overlaps with the temperature during the pump-off period. Consequently, after averaging, certain regions in the in-phase spatial profile exhibit sign changes (Extended Data Fig. 4j–l). In summary, the observed sign change in the in-phase signal in our experiments does not indicate negative differential temperatures but is a result of the lock-in operation. This phenomenon becomes apparent when there is a slow, highly viscous flow of heat, as seen in ultrathin layered semiconductors (see the main text) with low thermal diffusivity.
DFT details for LBTE calculations
As an independent approach, we perform LBTE calculations that include non-diffusive effects, except for thermoelastic effects. We follow the strategy presented in ref. 49 to solve the LBTE:
$$\frac{\partial {g}_{n}}{\partial t}+{{\bf{v}}}_{n}\cdot \nabla {g}_{n}=Q{p}_{n}+\mathop{\sum }\limits_{j}{\omega }_{n}{W}_{n,j}\frac{1}{{\omega }_{j}}({c}_{j}\Delta T-{g}_{j}),$$
(4)
where gn is the deviational phonon energy density per mode, Q is the macroscopic volumetric heat generation rate and the values pn describe the distribution of heating among the phonon modes. The continuous integral of the collision term in the Boltzmann transport equation has been discretized as matrix Wn,j, which is a general scattering matrix of dimensions M × M describing the scattering rate between phonon states n and j, acting on the difference between the equilibrium and non-equilibrium distribution functions. The heat generation rate expressed as Q = δ(t)χ(x), where δ(t) is the Dirac delta, corresponds to impulse laser heating with a spatial intensity profile χ(x). The temperature T is defined49 as the value for which the equilibrium energy density of phonons matches the non-equilibrium energy density. Linearization gives the temperature rise ΔT as the ratio of the non-equilibrium energy density of phonons divided by the heat capacity Cv:
$$\Delta T=\frac{1}{{C}_{{\rm{v}}}}\mathop{\sum }\limits_{n}{g}_{n}.$$
(5)
To obtain the necessary inputs for the LBTE and ballistic predictions, we performed DFT calculations using the QUANTUM ESPRESSO package. Several pseudopotentials and exchange–correlation functionals were tested including the full relativistic norm-conserving pseudopotential, ultrasoft pseudopotentials and standard solid-state pseudopotentials with the Perdew–Burke–Ernzerhof generalized gradient approximation and PBEsol. The DFT calculation parameters used in this work are as follows: a 20 × 20 × 1, 20 × 20 × 1, 16 × 16 × 16 and 16 × 16 × 16 Monkhorst–Pack k-mesh with a kinetic energy cut-off of 110, 120, 100 and 100 Ry and convergence criteria of 1 × 10−15, 1 × 10−16, 1 × 10−14, 1 × 10−14 Ry is used for monolayer MoS2, monolayer MoSe2, bulk MoS2 and bulk MoSe2, respectively. 8 × 8 × 1, 8 × 8 × 1, 4 × 4 × 4 and 4 × 4 × 4 second-order force constant supercells and 4 × 4 × 1, 4 × 4 × 1, 4 × 4 × 4 and 4 × 4 × 4 third-order force constants up to the fifth nearest neighbour were used, where only wavefunctions at the gamma point were considered for extraction of the third-order force constants for monolayer MoS2, monolayer MoSe2, bulk MoS2 and bulk MoSe2, respectively. DFT and Boltzmann transport equation input and output files are available via GitHub at https://github.com/jamalabouhaibeh/Calculations. The results for the ballistic regime are described in Supplementary Note 2 and Supplementary Figs. 4 and 5.