Main

During the last few years, signatures of viscous charge flow in the so-called Fermi-liquid hydrodynamic regime were observed in two-dimensional (2D) electron systems, especially graphene, using electrical device measurements7,8,9,11,12 and scanning probe microscopy10,13,22. A second hydrodynamic regime, which has no analogue in classical fluids, can occur very close to the Dirac point. When the Fermi temperature (TF = EF/kB, where EF is the Fermi energy and kB is the Boltzmann constant) becomes small compared to the electron temperature Te, the system becomes a quantum-critical fluid3,6,14,15,17. In this Dirac-fluid regime, the non-relativistic description of the viscous fluid is replaced by its ultra-relativistic counterpart, which accounts for the presence of both particles and holes, as well as for their linear energy dispersion. In line with theoretical predictions in this regime15, electrical measurements at cryogenic temperatures indicated a deviation from the Wiedemann–Franz law19 and from the Mott relation20, and a terahertz-probe study revealed the quantum-critical carrier scattering rate21.

Here, we follow electronic heat flow in the diffusive and hydrodynamic regimes at room temperature, and demonstrate a controlled Fermi-liquid to Dirac-fluid crossover, with a strongly enhanced thermal diffusivity close to the Dirac point. These observations are enabled by ultrafast spatiotemporal thermoelectric microscopy, a technique inspired by all-optical spatiotemporal diffusivity measurements23,24,25, with the crucial difference that the observable is the thermoelectric current, which is directly, and exclusively, sensitive to electronic heat26. We use a hexagonal boron nitride (hBN)-encapsulated graphene device that is both a Hall bar for electrical measurements and a split-gate thermoelectric detector (Fig. 1a). Since we use ultrashort laser pulses, with an approximate instrument response time (ΔtIRF where IRF means instrument response function) of 200 fs, to generate electronic heat, we are able to examine the system before momentum relaxation occurs, as we measure a momentum relaxation time, τmr, around 350 fs (Extended Data Fig. 1). In this temporal regime before momentum is relaxed, we enter the hydrodynamic window, because the electron–electron scattering time τee is <100 fs (ref. 27), that is τee < ΔτIRF < τmr. This is a different approach compared to most previous studies, where hydrodynamic effects were observed by using small system dimensions L to eliminate effects of momentum relaxation, that is vFτee < L < vFτmr (refs. 7,8,9,10,11,12,13,19,22) (vF = 106 m s−1 is the Fermi velocity). Our approach furthermore exploits elevated carrier temperatures, which greatly increases the accessibility of the Dirac-fluid regime, as for increasing carrier temperatures the crossover occurs increasingly far away from the Dirac point14,17 (Fig. 1b). As we will show, during the hydrodynamic window substantially more efficient heat spreading occurs in the Dirac-fluid regime than in the Fermi-liquid regime and in the diffusive regime (Fig. 1c,d).

Fig. 1: Spatiotemporal thermoelectric microscopy and heat spreading regimes.
figure 1

a, Concept of the experiment, where a graphene Hall bar and thermoelectric device is illuminated by two femtosecond heat-generating pulses with a relative temporal offset ∆t and a symmetric spatial offset ∆x with respect to the pn junction where electronic heat generates a thermoelectric current. The junction is created by applying +∆U to one gate and −∆U to the other. We isolate the differential thermoelectric current corresponding to light-induced electronic heat from both pulses that has travelled to the junction, where the heat adds up in a non-linear fashion. b, Phase diagram of the Dirac-fluid regime, calculated following ref. 14. For increasing Te, the Dirac-fluid regime occurs increasingly far away from the Dirac point. c,d, Illustration of light-triggered spreading of electronic heat in the Fermi-liquid regime (c) and Dirac-fluid regime (d). In both cases, for Δt > τmr, diffusive transport dominates (straight blue lines), while in the hydrodynamic window, with Δt < τmr, extremely efficient heat transport occurs in the Dirac-fluid regime (wavy red lines). e,f, Sketch of the spatial broadening of the heat spots for low (e) and high (f) diffusivity, D, indicating more interacting heat at the junction region, hence higher ∆ITE signal for higher D.

Full size image

Our technique works by using two ultrafast laser pulses that produce localized spots of electronic heat within tens of femtoseconds27. These spots are characterized by an increased carrier temperature Te > Tl, with Tl the lattice temperature (300 K). The degree of spatial spreading of these electronic heat spots as a function of time is governed by the diffusivity D. We control the relative spatial and temporal displacement, Δx and Δt, of the two pulses with sub-100-nm spatial precision and roughly 200 fs temporal resolution. Each laser pulse is incident on opposite sides of a pn junction at a distance ∆x/2 from the junction. This pn junction is created by applying opposite voltages ±ΔU with respect to the Dirac point voltages to the two backgates that form a split-gate structure. The two photo-generated electronic heat spots spread out spatially and part of the heat can reach the pn junction after a certain amount of time, generating a thermoelectric current at the junction through the Seebeck gradient26. The small region of the pn junction thus serves as a local probe of the electron temperature. While each of the heat spots can create thermoelectric current independently, we obtain spatiotemporal information by examining exclusively the signal that corresponds to heat generated by one of the pulses interacting with heat generated by the other pulse—the interacting heat current, ΔITE. Since the thermoelectric photocurrent scales sublinearly with incident power, we can isolate this interacting heat current ΔITE by modulating each laser beam at a different frequency, f1 and f2, and demodulating the thermoelectric current at the difference frequency f1 − f2. As illustrated in Fig. 1e,f, the higher the diffusivity D, the more interacting heat current ΔITE remains for increasing Δx and Δt.

Figure 2a shows the measured interacting heat current ∆ITE as a function of Δx and Δt. As expected, the largest ∆ITE occurs for the largest spatiotemporal overlap at the pn junction (Δx = Δt = 0). For increasing |∆t|, we find that the normalized signal extends further spatially, indicating the occurrence of heat spreading (Fig. 2b). This spatial spread is quantified via the second moment <Δx2>, which quantifies the width of the profile at different time delays (Methods). Similar to recent all-optical spatiotemporal microscopy24,25, we obtain spatial information beyond the diffraction limit by precise spatial sampling of diffraction-limited profiles. The experimentally obtained spatial spread as a function of Δt (Fig. 2c) is very similar to the calculated results (Fig. 2d), obtained by simulating the experiment with a given diffusivity, D (Methods and Supplementary Note 1). The white lines indicate the values of the spatial spread <Δx2> for different Δt. We also compare the simulated spatial spread <Δx2> versus Δt (blue dashed line in Fig. 2e) with the theoretical expectation according to the heat diffusion equation, <Δx2> = <Δx2>focus + 2DΔt (dash-dotted line in Fig. 2e). Here, D is the same diffusivity that was used as input for the simulation, and <Δx2>focus is the minimum second moment from the two overlapping pulses (Supplementary Note 2 and Supplementary Figs. 1–4). The initial slope is the same for both the simulated heat spreading and the theoretical spreading following the heat diffusion equation.

Fig. 2: Spatiotemporal tracking of heat spreading.
figure 2

a, The experimental spatiotemporal differential thermoelectric current ∆ITE as a function of Δx and Δt. b, Normalized profiles (norm.), showing a larger spatial spread for larger |∆t|. c,d, Experimental (c) and simulated (d) normalized ∆ITE for each ∆t, showing spatial broadening due to thermal transport as a function of ∆t. The white line indicates the spatial spread <∆x2>. e, Spatial spread <∆x2> of ∆ITE, as a function of ∆t for three different Fermi energies (symbols), with simulation results using as input the diffusivities from electrical mobility measurements (purple solid lines), with offset due to ultrafast heat spreading around time zero. Simulation (blue dashed line) and theoretical heat equation (black dash-dotted line) results with the same input diffusivity and no ultrafast spreading around time zero. Heat spreading with ultrahigh diffusivity in the Dirac-fluid regime (red line), which lasts for a few hundred fs, explains the time zero offset.

Full size image

We first discuss the experimental results in the diffusive regime, where |Δt| > τmr. For three different gate voltage combinations, corresponding to Fermi energies between 75 and 190 meV (TF = 900 − 2,200 K), we extract the spatial spread as a function of time delay (symbols in Fig. 2e), and compare it with the results from simulations (solid lines). For these simulations, we have used the diffusivity values that we obtain directly from electrical measurements of charge mobility μ on the same device (Extended Data Fig. 1), and the relation between mobility and diffusivity: D = μEF/2e (Methods). We find excellent agreement, if we account for short-lived ultrafast heat spreading around Δt = 0, which leads to a larger-than-expected initial spread at time zero <Δx2>min, as we will explain below. The agreement between the measured heat spread for |Δt| > τmr and the one calculated using the measured charge mobilities shows that electronic heat and charge flow together, as expected in the diffusive regime. Furthermore, it confirms that our technique is a reliable method for obtaining thermal diffusivities in a quantitative manner.

We now turn to the non-diffusive regime, by exploring the behaviour in the hydrodynamic window, where |Δt| < τmr. The experimentally obtained spatial spreads start at a minimum value <Δx2>min larger than 2 μm2, rather than starting at an expected <Δx2>focus = 0.56 μm2. A second device of hBN-encapsulated graphene with similar mobility reproduces this larger-than-expected spatial spread at time zero (Supplementary Note 3 and Extended Data Fig. 2). We exclude the possibility of an experimental artefact such as an underestimation of the laser spot size, since we repeated the measurements while scanning through the laser focus, and measured the focus size (Supplementary Figs. 1–4). Furthermore, we observe that the offset depends on the Fermi energy, while keeping all other experimental parameters fixed. Finally, we measured a third device with a lower charge mobility and shorter hydrodynamic time window: τmr < 100 fs, which is smaller than ΔtIRF. This device exhibits systematically less heat spreading around time zero (Supplementary Note 4 and Extended Data Fig. 3), consistent with its smaller hydrodynamic time window. We therefore attribute the large experimentally observed minimum <Δx2>min in Fig. 2e to ultrafast initial heat spreading that occurs before momentum relaxation takes place, Δt ≲ 350 fs (see the schematic illustration of spatiotemporal heat spreading in Fig. 1d). The dynamics of this initial heat spreading are washed out by the finite time resolution ΔtIRF, and manifests as a large minimum <Δx2>min at time zero. The observed initial spatial spread suggests a thermal diffusivity of D = (<Δx2>min − <Δx2>focus)/2ΔtIRF ≅ 70,000 cm2 s−1 for the lowest measured EF of 75 meV. Simulations of heat spreading with an input diffusivity of 100,000 cm2 s−1 are indeed consistent with the experimentally observed spread in the hydrodynamic window (the red line in Fig. 2e).

We attribute this observation of highly efficient initial heat spreading to the presence of the quantum-critical electron-hole plasma. We can exclude that the observed initial spreading is the result of ballistic transport, as we calculate that the ballistic contribution to initial heat spreading would give only <Δx2>ball = 0.68 μm2 (Supplementary Note 2 and Extended Data Fig. 4). Besides, ballistic transport has a very weak dependence (<10%) on carrier density in this range, as the Fermi velocity does not change appreciably for the Fermi energies considered here28. The reason for the high diffusivity in the Dirac-fluid regime is that the hot electrons and hot holes that coexist in this regime move in the same direction under a thermal gradient, with inter-particle scattering events conserving total momentum19. We note that typical transport measurements probe the sum of the momentum-conserving thermal resistivity (due to carrier–carrier scattering) and the momentum-non-conserving thermal resistivity (due to carrier–impurity and carrier–phonon scattering), where the latter dominates at room temperature. The ability of our technique to interrogate the system during the 350 fs before any momentum-non-conserving scattering occurs, means that this contribution to the overall resistivity is negligible. Therefore, we probe exclusively the momentum-conserving thermal conductivity, which diverges towards the Dirac point and towards infinite electron temperature.

To provide further evidence of hydrodynamic heat transport, we demonstrate the ability to control the crossover between the Fermi-liquid and quantum-critical Dirac-fluid regimes via the ratio Te/TF, by independently varying Te via the incident laser power and TF via the applied gate voltage. A larger ratio results in less Coulomb screening and correspondingly stronger hydrodynamic effects due to electron–electron interactions. If Te is larger than TF, electrons and holes coexist and the Dirac-fluid regime becomes accessible (Fig. 1b). We perform spatial scans in the hydrodynamic window at a temporal delay of Δt = 0, in a geometry with one laser pulse impinging on the junction, while scanning the other pulse across (x axis) and along (y axis) the junction region. Figure 3a–d shows four representative spatial ∆ITE maps with varying Te/TF, yet similar signal magnitudes. Clearly, the signal is broader for larger Te/TF, indicating faster thermal transport. We repeat these measurements for a range of Te and TF values and quantify the initial heat spreading using Gaussian functions, with widths σx and σy, to describe ∆ITE at Δt = 0 as a function of Δx or Δy (Fig. 3e,f and Supplementary Fig. 5). As expected for a crossover from the diffusive Fermi-liquid regime to the hydrodynamic Dirac-fluid regime, both spatial spreads σx and σy increase substantially for increasing ratio Te/TF. These spreads correspond to a diffusivity up to 40,000 cm2 s−1 (Methods), similar to the 70,000 cm2 s−1 we found earlier.

Fig. 3: Fermi-liquid to Dirac-fluid crossover.
figure 3

a–d, Time-zero spatial maps of ∆ITE for low optical power P and high gate voltage ∆U with an np junction (a) and with a pn junction (d); and for high optical power and low gate voltage with an np junction (b), and with a pn junction (c). For larger ratio Te/TF (that is, larger P/∆U) the spatial extent is clearly larger. e,f, Time zero Gaussian widths for spatial scans with one pulse on the junction and the second one scanning across (e) and along (f) the graphene pn junction, as a function of P and ∆U. The strong dependence on Te/TF demonstrates our ability to transition under control into the Dirac-fluid regime with strongly increased thermal diffusivity. g,h, Calculation of the thermal diffusivity following refs. 6,18 with only electron–electron interactions (g) and only long-range Coulomb scattering (h). The contours in g are the calculated time zero spreads \(\sigma _{{{{\mathrm{calc}}}}}^2\) (Methods).

Full size image

We compare our experimental results to Boltzmann transport calculations following refs. 6,18, including carrier interactions and long-range impurity scattering. We model impurities as Thomas–Fermi screened Coulomb scatterers of density 0.24 × 1012 cm−2. Figure 3g shows the calculated thermal diffusivity D as a function of TF and Te, when considering only the hydrodynamic term due to electron–electron interactions, relevant in the hydrodynamic window where |Δt| < τmr. A higher electron temperature, or lower Fermi temperature, leads to strongly increased diffusivity, signalling a crossover from the Fermi-liquid to the Dirac-fluid regime. This is the same qualitative trend as for the experimental data taken at Δt = 0 in Fig. 3e–f, where a larger initial width originates from a larger diffusivity, thus supporting our interpretation of a hydrodynamic crossover.

A more quantitative comparison shows that the calculated D in the diffusive regime is around 2,000 cm2 s−1 (Fig. 3h), in quantitative agreement with the experiment in the diffusive regime. The obtained thermal diffusivity in the hydrodynamic window close to the Dirac point reaches values above 100,000 cm2 s−1, even higher than our experimental estimates of 35,000–70,000 cm2 s−1. Using the calculated diffusivities, we estimate the spatial spread at time zero σcalc (Methods), as shown in Fig. 3g. These are similar to the experimentally obtained ones, thus confirming our conclusion of highly efficient heat spreading in the Dirac-fluid regime at room temperature, with a diffusivity that is almost two orders of magnitude larger than in the diffusive regime. We note that the theoretical calculations predict that even higher diffusivities are attainable.

Finally, we discuss the three-dimensional (3D) thermal conductivity, to assess the ability to transport useful amounts of heat. We find roughly 100 W m−1 K−1 in the diffusive regime (Methods), in agreement with ab initio calculations29. In the Dirac-fluid regime, with an electron temperature of roughly 1,000 K, we obtain a thermal conductivity of 18,000–40,000 W m−1 K−1. This is in agreement with ref. 15, where values up to 100,000 W m−1 K−1 were predicted theoretically for large Te/TF. The thermal conductivity we obtain is about three orders of magnitude larger than the one obtained in the Dirac-fluid regime at cryogenic temperatures19. Our results show that in the Dirac-fluid window the electronic contribution to heat transport can be much larger than the phononic contribution with a conductivity of >2,000 W m−1 K−1 (ref. 30), which is already exceptionally high and can also be enhanced hydrodynamically, as shown recently31. Thus, the Dirac electron-hole plasma can contribute strongly to thermal transport, extracting heat from hot spots much faster than predicted by classical limits.

In conclusion, our results show that the—until recently unreachable—physical phenomena associated with the Dirac fluid do not only offer an exciting playground for interesting physical phenomena, but also hold great promise for applications, for example in thermal management of nanoscale devices. We note that the quantum-critical behaviour can be switched on and off using a modest gate voltage and in systems prepared by standard fabrication techniques. Finally, we believe that the optoelectronic technique we have introduced—with the potential of increased spatial accuracy and temporal resolution—will be a valuable tool to reach a better understanding of the thermal behaviour of a broad range of quantum materials, with great promise for new technological applications.

Methods

Fabrication of split-gate thermoelectric device

The split-gate device with Hall geometry consists of exfoliated, single layer graphene encapsulated by hBN, prepared using standard exfoliation and dry transfer techniques. The hBN-graphene-hBN stack is placed on a predefined split-gate structure made of graphene, grown by chemical vapour deposition, where the gap between the two gates is roughly 100 nm, created via electron-beam lithography and reactive ion etching (RIE). The top hBN and graphene are etched into a Hall bar shape with laser lithography and RIE, keeping the split-gate intact and not etching completely through the bottom hBN. Finally, the Ti/Au side contacts are created by a further step of lithography, RIE and metal evaporation. The fabrication steps are shown in Supplementary Fig. 6.

Spatiotemporal thermoelectric current microscopy setup

Our setup enables us to follow electronic heat spreading in space and time, because we use the thermoelectric signal generated by electronic heat interacting at a fixed location (the pn junction), while we vary the spatial displacement of our two laser pulses with respect to this junction and vary the temporal delay between the two ultrashort pulses. This means that we are following in space and time the diffusion of light-induced electronic heat from the location of light incidence to the pn junction. It is the thermoelectric effect at the pn junction, governed by the Seebeck coefficient, that generates our observable signal, the thermoelectric current. We note that although the value of the Seebeck coefficient itself changes when changing EF, and when entering the hydrodynamic regime20, this only affects the magnitude of the thermoelectric current, rather than how electronic heat is diffusing outside the pn junction, which is what we are following with our spatiotemporal technique.

A sketch of the setup is shown in Supplementary Fig. 7. A Ti:sapphire oscillator (886-nm centre wavelength, 76-MHz repetition rate), is split into two beam paths. Both beams are modulated with optical choppers, at frequencies f1 = 741 and f2 = 529 Hz. The relative time delay between the two pulses is controlled by a mechanical delay line. The spatial offset of one beam with respect to the other is controlled with a mirror galvanometer, while the position of the sample with respect to the beams is controlled with a piezo scanning stage. The beams are focused onto the sample with a ×40 0.6 numerical aperture objective lens. We collect the thermoelectric (TE) photocurrent between the source and drain contacts through the graphene sheet on either side of the junction via lock-in amplification. By demodulating the current signal at the difference frequency of the two modulation frequencies, f2 − f1 = 211.7 Hz, we isolate the signal caused by the interaction of both heating sources, which we call the interacting heat current ∆ITE. The temporal resolution of the setup of 200 fs is determined in the sample plane of the microscope (Supplementary Fig. 8). The spatial resolution defined by our spot sizes is below 1 μm, whereas the accuracy with which we can observe electronic heat spreading is determined by the signal-to-noise ratio, and is estimated to be below 100 nm.

Estimating Fermi temperature controlled by gate voltage

During photocurrent measurements, the gate voltage Ux is applied to one side (x is ‘A’) or the other side (x is ‘B’) side of the split gate. We always apply a symmetric voltage around the experimentally determined Dirac point voltage \(U_{{{\mathrm{x}}}}^{{{{\mathrm{DP}}}}}\): \(U_{{{\mathrm{A}}}} = U_{{{\mathrm{A}}}}^{{{{\mathrm{DP}}}}} + {\Delta} U\) and \(U_{{{\mathrm{B}}}} = U_{{{\mathrm{B}}}}^{{{{\mathrm{DP}}}}} - {\Delta} U\). The gate electrode and the graphene form a capacitor with the dielectric hexagonal boron nitride (hBN), with a thickness of thBN = 70 nm, and a relative permittivity of \({\it{\epsilon }}_{{{{\mathrm{hBN}}}}} = 3.56\). The carrier density n is calculated via \(n = \frac{{{\it{\epsilon }}_0{\it{\epsilon }}_{{{{\mathrm{hBN}}}}}}}{{e\,t_{{{{\mathrm{hBN}}}}}}}{\Delta} U\), where \(\epsilon_{0}\) is the vacuum permittivity. We calculate the Fermi energy EF and the Fermi temperature TF via \(E_{{{\mathrm{F}}}}^2\) = πħ2\(v_{{{\mathrm{F}}}}^2 \cdot n\) and \(T_{{{\mathrm{F}}}} = \frac{{E_{{{\mathrm{F}}}}}}{{k_{{{\mathrm{B}}}}}},\) where kB is the Boltzmann constant.

Estimating carrier temperature controlled by laser power

The thermoelectric photovoltage is assumed to be proportional to the time-averaged increase of the electronic temperature Te above the ambient temperature T0, as in ref. 32. The sublinear dependence of the thermoelectric current ITE on optical power for the device under study here for illumination with a single pulsed laser (λ = 886 nm) is shown in Supplementary Fig. 9. With a linear temperature scaling of the electronic heat capacity for graphene away from the Dirac point, Ce(T) = γT, we integrate the heat energy per unit area dQ = CedT, that is, \(\mathop {\smallint }\limits_{Q_0}^{Q_0 + {\Delta} Q} {{{\mathrm{d}}}}Q = \mathop {\smallint }\limits_{T_0}^{T_e} \gamma T{{{\mathrm{d}}}}T\). With the incident power P proportional to the absorbed heat energy per unit area ∆Q, we find that the peak Te as a function of the laser power P scales as in ref. 32, \(T_e = \root {2} \of {{T_0^2 + bP}}\). Here, the parameter b is defined via bP = 2∆Q/γ, and is used to convert incident power to peak electron temperature (Supplementary Fig. 9).

Simulation of the experiment

A detailed description of the simulation can be found in Supplementary Note 1 and Supplementary Fig. 10. In brief, we calculate the spatiotemporal evolution of electronic heat generated by the two optical pulses in the graphene sheet via the heat equation with a finite difference method. We define Gaussian heating pulses and calculate their temperature rise via the experimentally measured non-linear power scaling. We extract the differential TE current contribution as a function of ∆x and ∆t by the difference of the heating at the pn junction region in the presence of both pulses with respect to simulations with only one pulse at a time, analogous to the experimental difference-frequency demodulation.

Quantifying the spatial spread

The following analysis is conducted both on the experimental and the simulated data of ∆ITE(∆x, ∆t) for ‘symmetric experiments’ with optical pulses incident at a distance ∆x on each side of the pn-junction (Figs. 1 and 2). For each ∆t of the datasets ∆ITE(∆x, ∆t) we calculate the width of the signal via the second moment, which for an ideal Gaussian profile is equal to the squared Gaussian width σ2. The second moment is calculated from the pixels ∆xi (i = 1,…, N) via

$$\begin{array}{rcl}< {\Delta} x^2 > ({\Delta} t) & = & \frac{{\mathop {\sum }\nolimits_i \left| {{\Delta} x_i - \overline {{\Delta} x} } \right|^2{\Delta} I_{{{{\mathrm{TE}}}}}({\Delta} x_i,{\Delta} t)}}{{\mathop {\sum }\nolimits_i {\Delta} I_{{{{\mathrm{TE}}}}}({\Delta} x_i,\,{\Delta} t)}},\,{{{\mathrm{with}}}}\,{{{\mathrm{the}}}}\,{{{\mathrm{mean}}}} \\\overline {{\Delta} x} & = & \frac{{\mathop {\sum }\nolimits_i {\Delta} x_i{\Delta} I_{{{{\mathrm{TE}}}}}({\Delta} x_i,{\Delta} t)}}{{\mathop {\sum }\nolimits_i {\Delta} I_{{{{\mathrm{TE}}}}}(x_i,t)}}.\end{array}$$

We note that the minimum second moment at the focus <∆x2>focus of 0.56 µm2 comes from simulating the symmetric experiment, using as input the measured Gaussian beam width at the focus σfocus2 = 0.14 µm2 (Supplementary Note 2). For the ‘asymmetric experiments’ with one optical pulse always incident on the pn junction (data for Fig. 3), we always consider the spatial profile only at time zero. Here we find that Gaussian fits with a background give the most reliable results. The entire set of data is shown in Supplementary Fig. 5. For each dataset ∆ITE(∆x) or ∆ITE(∆y) taken at ∆t = 0, we perform Gaussian fits using the function \({{{\mathrm{f}}}}({\Delta} x) = a\exp \left( { - \frac{{{\Delta} x^2}}{{2\sigma ^2}}} \right) + b\), where the Gaussian squared width σ2 indicates the thermal spreading. Here, the minimum simulated Gaussian widths are (σx2)focus = 0.34 µm2 and (σy2)focus = 0.44 µm2 (Supplementary Note 2). The experimentally obtained widths from this dataset as function of gate voltage and optical power are also shown in Supplementary Fig. 5, showing an increase with power, that is a larger Te, and an increase towards the Dirac point, that is a smaller TF. We estimate the theoretical Gaussian widths in Fig. 3g using \(\sigma _{{{{\mathrm{calc}}}}}^2\) = (σx2)focus + 2D ΔtIRF, where D are the calculated diffusivities.

Electrical measurements

We characterize our device electrically with four-probe measurements (Extended Data Fig. 1), finding a charge mobility μ of 30,000–50,000 cm2 Vs−1, depending on carrier density. The measured mobilities correspond to a momentum relaxation time 𝜏mr of 300–500 fs. These relaxation times are longer than the temporal resolution (the IRF) of our measurement technique, ΔtIRF ≅ 200 fs, thus allowing us to probe our system before and after momentum relaxation occurs, that is in the non-diffusive and diffusive regime. We use these measured charge mobilities to calculate the expected thermal diffusivity via the Einstein relation33,34 \(\mu _{{{{\mathrm{e}}}}/{{{\mathrm{h}}}}} = \frac{e}{{n_{{{{\mathrm{e}}}}/{{{\mathrm{h}}}}}}}\frac{{\partial n_{{{{\mathrm{e}}}}/{{{\mathrm{h}}}}}}}{{\partial E_{{{\mathrm{F}}}}}}D_{{{{\mathrm{e}}}}/{{{\mathrm{h}}}}}\), where e is the elementary charge, EF is the Fermi energy and ne/h is the electron/hole carrier density. For highly doped graphene (\(E_{{{\mathrm{F}}}} \gg k_{{{\mathrm{B}}}}T\)) the carrier density expression \(n_{{{{\mathrm{e}}}}/{{{\mathrm{h}}}}} = \frac{{E_{{{\mathrm{F}}}}^2}}{{\pi \hbar ^2v_{{{\mathrm{F}}}}^2}}\), leads to the simple relation:\(D_{{{{\mathrm{e}}}}/{{{\mathrm{h}}}}} = \frac{{E_{{{\mathrm{F}}}}}}{{2e}}\mu _{{{{\mathrm{e}}}}/{{{\mathrm{h}}}}}\). We note that we obtain the identical result by calculating D from the ratio of the 2D thermal conductivity κe,2D and the electronic heat capacity Ce and using the Wiedemann–Franz law: \(\kappa _{{{{\mathrm{e}}}},2{{{\mathrm{D}}}}}/\sigma = \pi ^2/3\cdot (k_{{{\mathrm{B}}}}/e)^2T_{{{\mathrm{e}}}}\), where kB is the Boltzmann constant and e the elementary charge, together with the conductivity σ = neμ and the following heat capacity for graphene (valid for Te < TF): \(C_{{{\mathrm{e}}}} = \frac{{2\pi \varepsilon _{{{\mathrm{F}}}}k_{{{\mathrm{B}}}}^2T_{{{\mathrm{e}}}}}}{{3\hbar ^2v_{{{\mathrm{F}}}}^2}}\). Given the measured mobilities, we expect thermal diffusivities around 2,000 cm2 s−1 for our sample.

Thermal diffusivity and conductivity of the Dirac fluid

We estimate the enhanced thermal diffusivity of the Dirac fluid by comparing the measured width at time zero <Δx2>min to the expected width <Δx2>focus explained above, via D = (<Δx2>min − <Δx2>focus)/2ΔtIRF. We find values up to 74,000 cm2 s−1 for the symmetric scan (Fig. 2), and 29,000 and 39,000 cm2 s−1 for the x and y directions of the asymmetric scan (Fig. 3), respectively, where <Δx2> is replaced with (σx2) and (σy2), respectively. The same calculation for a second device (Supplementary Note 3 and Extended Data Fig. 2) gives a diffusivity of 100,000 cm2 s−1. The 3D thermal conductivity κ3D of the Dirac fluid is calculated from the diffusivity D and the electronic heat capacity Ce, via κ3D = DCe/d, where d is the thickness of graphene, 0.3 nm. For the Dirac fluid, we have Te > TF, and therefore use the ‘undoped’ electronic heat capacity35 \(\frac{{18\,\zeta (3)}}{{\pi (\hbar v_{{{\mathrm{F}}}})^2}}k_{{{\mathrm{B}}}}^3T_{{{\mathrm{e}}}}^2\), where \(\zeta \left( 3 \right) \approx 1.202\). With the above estimate D = 35,000–70,000 cm2 s−1 and Te = 1,000 K, we obtain the 3D thermal conductivity κ3D = 18,000–40,000 W mK−1. This corresponds to a 2D κ2D > 5 μW K−1. This value is orders of magnitude larger than the value found in ref. 19. The reason for this is that our electron temperature is more than ten times higher, and therefore the electronic heat capacity is >100× higher. Furthermore, we reach a Te/TF > 3, while their maximum Te/TF was around two, which means that we are further in the Dirac-fluid regime with its diverging thermal diffusivity.

Dirac-fluid crossover temperature

Following the treatment in ref. 14, we find the crossover temperature from Fermi liquid to Dirac fluid, as a function of Fermi temperature as

$$T_{{{{\mathrm{cross}}}}}(T_{{{\mathrm{F}}}}) = T_{{{\mathrm{F}}}}\left( {1 + \lambda \ln \left( {\frac{{T_0}}{{T_{{{\mathrm{F}}}}}}} \right)} \right),$$

where \(\lambda = {{{\mathrm{e}}}}^2/16\epsilon _0{\it{\epsilon }}_{{{\mathrm{r}}}}v_{{{\mathrm{F}}}}\hbar \approx 0.55/{\it{\epsilon }}_{{{\mathrm{r}}}}\) for graphene with the dielectric environment \({\it{\epsilon }}_{{{\mathrm{r}}}} \approx 3.56\) for hBN. The temperature \(T_0 = \frac{{2\hbar v_{{{\mathrm{F}}}}\sqrt \pi }}{{3^{3/4}k_{{{\mathrm{B}}}}a_0}} \approx 8.4\times 10^4\,{{{\mathrm{K}}}}\), with the inter-atomic distance \(a_0 = 1.42\times 10^{ - 10}\,{{{\mathrm{m}}}}\). The resulting crossover temperature is shown in Fig. 1b. We note that the relatively high refractive index of the hBN encapsulant makes the Dirac fluid more easily accessible, as it lowers the crossover temperature compared to vacuum, by a factor of about two for the range of TF values studied here.