Letters https://doi.org/10.1038/s41565-021-00957-6 Observation of giant and tunable thermal diffusivity of a Dirac fluid at room temperature Alexander Block 1,2, Alessandro Principi3, Niels C. H. Hesp 1, Aron W. Cummings2, Matz Liebel 1 , Kenji Watanabe 4, Takashi Taniguchi 5, Stephan Roche 2,6, Frank H. L. Koppens 1,6, Niek F. van Hulst 1,6 and Klaas-Jan Tielrooij 2 ✉ Conducting materials typically exhibit either diffusive or enhanced thermal diffusivity close to the Dirac point. These ballistic charge transport. When electron–electron interac- observations are enabled by ultrafast spatiotemporal thermoelec- tions dominate, a hydrodynamic regime with viscous charge tric microscopy, a technique inspired by all-optical spatiotempo- flow emerges1–13. More stringent conditions eventually yield ral diffusivity measurements23–25, with the crucial difference that a quantum-critical Dirac-fluid regime, where electronic heat the observable is the thermoelectric current, which is directly, and can flow more efficiently than charge14–22. However, observing exclusively, sensitive to electronic heat26. We use a hexagonal boron and controlling the flow of electronic heat in the hydrodynamic nitride (hBN)-encapsulated graphene device that is both a Hall bar regime at room temperature has so far remained elusive. Here for electrical measurements and a split-gate thermoelectric detec- we observe heat transport in graphene in the diffusive and tor (Fig. 1a). Since we use ultrashort laser pulses, with an approxi- hydrodynamic regimes, and report a controllable transition to mate instrument response time (ΔtIRF where IRF means instrument the Dirac-fluid regime at room temperature, using carrier tem- response function) of 200 fs, to generate electronic heat, we are able perature and carrier density as control knobs. We introduce to examine the system before momentum relaxation occurs, as we the technique of spatiotemporal thermoelectric microscopy measure a momentum relaxation time, τmr, around 350 fs (Extended with femtosecond temporal and nanometre spatial resolu- Data Fig. 1). In this temporal regime before momentum is relaxed, tion, which allows for tracking electronic heat spreading. In we enter the hydrodynamic window, because the electron–elec- the diffusive regime, we find a thermal diffusivity of roughly tron scattering time τee is <100 fs (ref. 27), that is τee < ΔτIRF < τmr. 2,000 cm2 s−1, consistent with charge transport. Moreover, This is a different approach compared to most previous studies, within the hydrodynamic time window before momentum where hydrodynamic effects were observed by using small system relaxation, we observe heat spreading corresponding to a dimensions L to eliminate effects of momentum relaxation, that is giant diffusivity up to 70,000 cm2 s−1, indicative of a Dirac vFτee < L < vFτmr (refs. 7–13,19,22) (vF = 106 m s−1 is the Fermi velocity). fluid. Our results offer the possibility of further exploration Our approach furthermore exploits elevated carrier temperatures, of these interesting physical phenomena and their potential which greatly increases the accessibility of the Dirac-fluid regime, applications in nanoscale thermal management. as for increasing carrier temperatures the crossover occurs increas- During the last few years, signatures of viscous charge flow in ingly far away from the Dirac point14,17 (Fig. 1b). As we will show, the so-called Fermi-liquid hydrodynamic regime were observed in during the hydrodynamic window substantially more efficient heat two-dimensional (2D) electron systems, especially graphene, using spreading occurs in the Dirac-fluid regime than in the Fermi-liquid electrical device measurements7–9,11,12 and scanning probe micros- regime and in the diffusive regime (Fig. 1c,d). copy10,13,22. A second hydrodynamic regime, which has no analogue Our technique works by using two ultrafast laser pulses that in classical fluids, can occur very close to the Dirac point. When the produce localized spots of electronic heat within tens of femtosec- Fermi temperature (TF = EF/kB, where EF is the Fermi energy and kB onds27. These spots are characterized by an increased carrier tem- is the Boltzmann constant) becomes small compared to the electron perature Te > Tl, with Tl the lattice temperature (300 K). The degree temperature Te, the system becomes a quantum-critical fluid3,6,14,15,17. of spatial spreading of these electronic heat spots as a function of In this Dirac-fluid regime, the non-relativistic description of the time is governed by the diffusivity D. We control the relative spa- viscous fluid is replaced by its ultra-relativistic counterpart, which tial and temporal displacement, Δx and Δt, of the two pulses with accounts for the presence of both particles and holes, as well as for sub-100-nm spatial precision and roughly 200 fs temporal resolu- their linear energy dispersion. In line with theoretical predictions tion. Each laser pulse is incident on opposite sides of a pn junction in this regime15, electrical measurements at cryogenic tempera- at a distance ∆x/2 from the junction. This pn junction is created tures indicated a deviation from the Wiedemann–Franz law19 and by applying opposite voltages ±ΔU with respect to the Dirac point from the Mott relation20, and a terahertz-probe study revealed the voltages to the two backgates that form a split-gate structure. The quantum-critical carrier scattering rate21. two photo-generated electronic heat spots spread out spatially and Here, we follow electronic heat flow in the diffusive and hydro- part of the heat can reach the pn junction after a certain amount of dynamic regimes at room temperature, and demonstrate a con- time, generating a thermoelectric current at the junction through trolled Fermi-liquid to Dirac-fluid crossover, with a strongly the Seebeck gradient26. The small region of the pn junction thus 1 ICFO (Institut de Ciències Fotòniques), The Barcelona Institute of Science and Technology, Castelldefels, Spain. 2Catalan Institute of Nanoscience and Nanotechnology (ICN2), BIST and CSIC, Bellaterra, Spain. 3School of Physics and Astronomy, University of Manchester, Manchester, UK. 4Research Center for Functional Materials, National Institute for Materials Science, Tsukuba, Japan. 5International Center for Materials Nanoarchitectonics, National Institute for Materials Science, Tsukuba, Japan. 6ICREA - Institució Catalana de Recerca i Estudis Avançats, Barcelona, Spain. ✉e-mail: klaas.tielrooij@icn2.cat Nature Nanotechnology | VOL 16 | November 2021 | 1195–1200 | www.nature.com/naturenanotechnology 1195 Letters Nature Nanotechnology a b Carrier density n (1012 cm–2) x y Δt –2 –1 0 1 2 ITE 2 Electron temp. Te (103 K) Dirac fluid F T = e T 1 Δx Fermi liquid 300 K 0 2 1 0 1 2 –ΔU +ΔU Fermi temperature TF (103 K) c e f Fermi liquid 0 fs 200 fs 400 fs 600 fs Low D High D t0 Hydrodynamic window Diffusive regime d Dirac fluid t1 ~ΔITE t2 Hydrodynamic window Diffusive regime Δx Fig. 1 | Spatiotemporal thermoelectric microscopy and heat spreading regimes. 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. serves as a local probe of the electron temperature. While each of Here, D is the same diffusivity that was used as input for the simu- the heat spots can create thermoelectric current independently, we lation, and <Δx2>focus is the minimum second moment from the obtain spatiotemporal information by examining exclusively the two overlapping pulses (Supplementary Note 2 and Supplementary signal that corresponds to heat generated by one of the pulses inter- Figs. 1–4). The initial slope is the same for both the simulated heat acting with heat generated by the other pulse—the interacting heat spreading and the theoretical spreading following the heat diffusion current, ΔITE. Since the thermoelectric photocurrent scales sublin- equation. early with incident power, we can isolate this interacting heat cur- We first discuss the experimental results in the diffusive rent ΔITE by modulating each laser beam at a different frequency, regime, where |Δt| > τmr. For three different gate voltage combina- f1 and f2, and demodulating the thermoelectric current at the dif- tions, corresponding to Fermi energies between 75 and 190 meV ference frequency f1 − f2. As illustrated in Fig. 1e,f, the higher the (TF = 900 − 2,200 K), we extract the spatial spread as a function of diffusivity D, the more interacting heat current ΔITE remains for time delay (symbols in Fig. 2e), and compare it with the results from increasing Δx and Δt. simulations (solid lines). For these simulations, we have used the dif- Figure 2a shows the measured interacting heat current ∆ITE as fusivity values that we obtain directly from electrical measurements a function of Δx and Δt. As expected, the largest ∆ITE occurs for of charge mobility μ on the same device (Extended Data Fig. 1), and the largest spatiotemporal overlap at the pn junction (Δx = Δt = 0). the relation between mobility and diffusivity: D = μEF/2e (Methods). For increasing |∆t|, we find that the normalized signal extends fur- We find excellent agreement, if we account for short-lived ultrafast ther spatially, indicating the occurrence of heat spreading (Fig. 2b). heat spreading around Δt = 0, which leads to a larger-than-expected This spatial spread is quantified via the second moment <Δx2>, initial spread at time zero <Δx2>min, as we will explain below. The which quantifies the width of the profile at different time delays agreement between the measured heat spread for |Δt| > τmr and the (Methods). Similar to recent all-optical spatiotemporal micros- one calculated using the measured charge mobilities shows that copy24,25, we obtain spatial information beyond the diffraction limit electronic heat and charge flow together, as expected in the diffu- by precise spatial sampling of diffraction-limited profiles. The sive regime. Furthermore, it confirms that our technique is a reliable experimentally obtained spatial spread as a function of Δt (Fig. 2c) method for obtaining thermal diffusivities in a quantitative manner. is very similar to the calculated results (Fig. 2d), obtained by We now turn to the non-diffusive regime, by exploring the simulating the experiment with a given diffusivity, D (Methods and behaviour in the hydrodynamic window, where |Δt| < τmr. The Supplementary Note 1). The white lines indicate the values of the experimentally obtained spatial spreads start at a minimum value spatial spread <Δx2> for different Δt. We also compare the simu- <Δx2>min larger than 2 μm2, rather than starting at an expected lated spatial spread <Δx2> versus Δt (blue dashed line in Fig. 2e) <Δx2>focus = 0.56 μm2. A second device of hBN-encapsulated gra- with the theoretical expectation according to the heat diffusion phene with similar mobility reproduces this larger-than-expected equation, <Δx2> = <Δx2>focus + 2DΔt (dash-dotted line in Fig. 2e). spatial spread at time zero (Supplementary Note 3 and Extended 1196 Nature Nanotechnology | VOL 16 | November 2021 | 1195–1200 | www.nature.com/naturenanotechnology Nature Nanotechnology Letters a ΔITE (nA) b 1.0 –2 Δt = 0 ps ∣Δt∣ = 3 ps 0.4 ΔITE (norm.) Δx (µm) 0 0.2 0.5 √<Δx 2> 2 0 0 –5 0 5 –4 0 4 Pulse delay Δt (ps) Spatial displacement Δx (µm) c d –2 1 –2 1 ΔITE (norm.) ΔITE (norm.) Δx (µm) Δx (µm) 0 0 2 0 2 0 –5 0 5 –5 0 5 Pulse delay Δt (ps) Pulse delay Δt (ps) e Diffusive regime Hydrodynamic window Diffusive regime 6 D=1 5 meV .6 × 1 –1 s EF = 7 0 3 cm 2 2 Spatial spread <Δx 2> (µm2) s –1 5 cm 30 me V 4 D=2 .1 × 1 =1 0 EF = 1 0 3 cm 2 D s –1 V 90 me D=2 .3 × 1 EF = 1 03 cm 2 s –1 Hydrodynamic 2 τms ΔtIRF focus 0 –4 –2 0 2 4 Pulse delay Δt (ps) Fig. 2 | Spatiotemporal tracking of heat spreading. 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. Data Fig. 2). We exclude the possibility of an experimental artefact of D = (<Δx2>min − <Δx2>focus)/2ΔtIRF ≅ 70,000 cm2 s−1 for the low- such as an underestimation of the laser spot size, since we repeated est measured EF of 75 meV. Simulations of heat spreading with an the measurements while scanning through the laser focus, and mea- input diffusivity of 100,000 cm2 s−1 are indeed consistent with the sured the focus size (Supplementary Figs. 1–4). Furthermore, we experimentally observed spread in the hydrodynamic window (the observe that the offset depends on the Fermi energy, while keep- red line in Fig. 2e). ing all other experimental parameters fixed. Finally, we measured We attribute this observation of highly efficient initial heat a third device with a lower charge mobility and shorter hydrody- spreading to the presence of the quantum-critical electron-hole namic time window: τmr < 100 fs, which is smaller than ΔtIRF. This plasma. We can exclude that the observed initial spreading is the device exhibits systematically less heat spreading around time zero result of ballistic transport, as we calculate that the ballistic contri- (Supplementary Note 4 and Extended Data Fig. 3), consistent with bution to initial heat spreading would give only <Δx2>ball = 0.68 μm2 its smaller hydrodynamic time window. We therefore attribute (Supplementary Note 2 and Extended Data Fig. 4). Besides, ballistic the large experimentally observed minimum <Δx2>min in Fig. 2e transport has a very weak dependence (<10%) on carrier density in to ultrafast initial heat spreading that occurs before momentum this range, as the Fermi velocity does not change appreciably for the relaxation takes place, Δt ≲ 350 fs (see the schematic illustration Fermi energies considered here28. The reason for the high diffusiv- of spatiotemporal heat spreading in Fig. 1d). The dynamics of this ity in the Dirac-fluid regime is that the hot electrons and hot holes initial heat spreading are washed out by the finite time resolution that coexist in this regime move in the same direction under a ther- ΔtIRF, and manifests as a large minimum <Δx2>min at time zero. mal gradient, with inter-particle scattering events conserving total The observed initial spatial spread suggests a thermal diffusivity momentum19. We note that typical transport measurements probe Nature Nanotechnology | VOL 16 | November 2021 | 1195–1200 | www.nature.com/naturenanotechnology 1197 Letters Nature Nanotechnology a b c d ΔITE (nA) ΔITE (nA) ΔITE (nA) ΔITE (nA) 0 0.2 0 0.15 –0.2 0 –0.1 0 P = 4 µW P = 20 µW P = 20 µW P = 4 µW 2 Δy (µm) 0 –2 ΔU = –8.5 V ΔU = –1.5 V ΔU = 1.5 V ΔU = 8.5 V –2 0 2 –2 0 2 –2 0 2 –2 0 2 Δx (µm) Δx (µm) Δx (µm) Δx (µm) Te < TF Te > TF Te < TF e Te (K) P (µW) Dirac fluid σx2 (µm2) 2,000 20 1.4 =T F Te 1,600 12 0.9 1,200 Fermi liquid Fermi liquid 4 0.4 800 f Te (K) P (µW) Dirac fluid σy2 (µm2) 2,000 20 1.8 1,600 12 1.4 1,200 Fermi liquid Fermi liquid 1.0 4 800 –8 –4 0 4 8 ΔU (V) 175 150 125 100 75 75 100 125 150 175 EF (meV) TF (K) 2,000 1,600 1,200 800 800 1,200 1,600 2,000 g Dirac fluid D (cm2 s–1) 2,000 2 2 µm 2 µm µm 2 105 2 1 µm 1,500 0.5 5 Te (K) 1,000 104 Fermi liquid Fermi liquid 500 103 h D (103 cm2 s–1) 2,000 4 1,500 3 Te (K) 1,000 2 500 1 1,800 1,400 1,000 600 600 1,000 1,400 1,800 TF (K) TF (K) the sum of the momentum-conserving thermal resistivity (due to The ability of our technique to interrogate the system during the carrier–carrier scattering) and the momentum-non-conserving 350 fs before any momentum-non-conserving scattering occurs, thermal resistivity (due to carrier–impurity and carrier–pho- means that this contribution to the overall resistivity is negli- non scattering), where the latter dominates at room temperature. gible. Therefore, we probe exclusively the momentum-conserving 1198 Nature Nanotechnology | VOL 16 | November 2021 | 1195–1200 | www.nature.com/naturenanotechnology Nature Nanotechnology Letters Fig. 3 | Fermi-liquid to Dirac-fluid crossover. 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 σ 2calc (Methods). thermal conductivity, which diverges towards the Dirac point and conductivity of 18,000–40,000 W m−1 K−1. This is in agreement towards infinite electron temperature. with ref. 15, where values up to 100,000 W m−1 K−1 were predicted To provide further evidence of hydrodynamic heat transport, theoretically for large Te/TF. The thermal conductivity we obtain is we demonstrate the ability to control the crossover between the about three orders of magnitude larger than the one obtained in the Fermi-liquid and quantum-critical Dirac-fluid regimes via the ratio Dirac-fluid regime at cryogenic temperatures19. Our results show Te/TF, by independently varying Te via the incident laser power and that in the Dirac-fluid window the electronic contribution to heat TF via the applied gate voltage. A larger ratio results in less Coulomb transport can be much larger than the phononic contribution with screening and correspondingly stronger hydrodynamic effects due a conductivity of >2,000 W m−1 K−1 (ref. 30), which is already excep- to electron–electron interactions. If Te is larger than TF, electrons tionally high and can also be enhanced hydrodynamically, as shown and holes coexist and the Dirac-fluid regime becomes accessible recently31. Thus, the Dirac electron-hole plasma can contribute (Fig. 1b). We perform spatial scans in the hydrodynamic window strongly to thermal transport, extracting heat from hot spots much at a temporal delay of Δt = 0, in a geometry with one laser pulse faster than predicted by classical limits. impinging on the junction, while scanning the other pulse across In conclusion, our results show that the—until recently unreach- (x axis) and along (y axis) the junction region. Figure 3a–d shows able—physical phenomena associated with the Dirac fluid do not four representative spatial ∆ITE maps with varying Te/TF, yet simi- only offer an exciting playground for interesting physical phe- lar signal magnitudes. Clearly, the signal is broader for larger Te/TF, nomena, but also hold great promise for applications, for example indicating faster thermal transport. We repeat these measurements in thermal management of nanoscale devices. We note that the for a range of Te and TF values and quantify the initial heat spreading quantum-critical behaviour can be switched on and off using a using Gaussian functions, with widths σx and σy, to describe ∆ITE at modest gate voltage and in systems prepared by standard fabrica- Δt = 0 as a function of Δx or Δy (Fig. 3e,f and Supplementary Fig. 5). tion techniques. Finally, we believe that the optoelectronic tech- As expected for a crossover from the diffusive Fermi-liquid regime nique we have introduced—with the potential of increased spatial to the hydrodynamic Dirac-fluid regime, both spatial spreads σx and accuracy and temporal resolution—will be a valuable tool to reach σy increase substantially for increasing ratio Te/TF. These spreads a better understanding of the thermal behaviour of a broad range correspond to a diffusivity up to 40,000 cm2 s−1 (Methods), similar of quantum materials, with great promise for new technological to the 70,000 cm2 s−1 we found earlier. applications. We compare our experimental results to Boltzmann transport calculations following refs. 6,18, including carrier interactions and Online content long-range impurity scattering. We model impurities as Thomas– Any methods, additional references, Nature Research report- Fermi screened Coulomb scatterers of density 0.24 × 1012 cm−2. ing summaries, source data, extended data, supplementary infor- Figure 3g shows the calculated thermal diffusivity D as a function mation, acknowledgements, peer review information; details of of TF and Te, when considering only the hydrodynamic term due author contributions and competing interests; and statements of to electron–electron interactions, relevant in the hydrodynamic data and code availability are available at https://doi.org/10.1038/ window where |Δt| < τmr. A higher electron temperature, or lower s41565-021-00957-6. Fermi temperature, leads to strongly increased diffusivity, signal- ling a crossover from the Fermi-liquid to the Dirac-fluid regime. Received: 1 January 2021; Accepted: 12 July 2021; This is the same qualitative trend as for the experimental data taken Published online: 23 August 2021 at Δt = 0 in Fig. 3e–f, where a larger initial width originates from a larger diffusivity, thus supporting our interpretation of a hydrody- References namic crossover. 1. Polini, M. & Geim, A. K. Viscous electron fluids. Phys. Today 73, 28 (2020). 2. Müller, M., Schmalian, J. & Fritz, L. Graphene: a nearly perfect fluid. A more quantitative comparison shows that the calculated D in Phys. Rev. Lett. 103, 025301 (2009). the diffusive regime is around 2,000 cm2 s−1 (Fig. 3h), in quantita- 3. Foster, M. S. & Aleiner, I. L. Slow imbalance relaxation and thermoelectric tive agreement with the experiment in the diffusive regime. The transport in graphene. Phys. Rev. B. 79, 085415 (2009). obtained thermal diffusivity in the hydrodynamic window close to 4. Narozhny, B. N., Gornyi, I. V., Titov, M., Schütt, M. & Mirlin, A. D. the Dirac point reaches values above 100,000 cm2 s−1, even higher Hydrodynamics in graphene: linear-response transport. Phys. Rev. B. 91, 035414 (2015). than our experimental estimates of 35,000–70,000 cm2 s−1. Using the 5. Levitov, L. & Falkovich, G. Electron viscosity, current vortices and negative calculated diffusivities, we estimate the spatial spread at time zero nonlocal resistance in graphene. Nat. Phys. 12, 672–676 (2016). σcalc (Methods), as shown in Fig. 3g. These are similar to the experi- 6. Zarenia, M., Principi, A. & Vignale, G. Disorder-enabled hydrodynamics of mentally obtained ones, thus confirming our conclusion of highly charge and heat transport in monolayer graphene. 2D Mater. 6, 035024 (2019). efficient heat spreading in the Dirac-fluid regime at room tempera- 7. Moll, P. J. W., Kushwaha, P., Nandi, N., Schmidt, B. & Mackenzie, A. P. Evidence for hydrodynamic electron flow in PdCoO2. Science 351, ture, with a diffusivity that is almost two orders of magnitude larger 1061–1064 (2016). than in the diffusive regime. We note that the theoretical calcula- 8. Bandurin, D. A. et al. Negative local resistance caused by viscous electron tions predict that even higher diffusivities are attainable. backflow in graphene. Science 351, 1055–1058 (2016). Finally, we discuss the three-dimensional (3D) thermal conduc- 9. Krishna Kumar, R. et al. Superballistic flow of viscous electron fluid through tivity, to assess the ability to transport useful amounts of heat. We graphene constrictions. Nat. Phys. 13, 1182–1185 (2017). 10. Braem, B. A. et al. Scanning gate microscopy in a viscous electron fluid. find roughly 100 W m−1 K−1 in the diffusive regime (Methods), in Phys. Rev. B. 98, 241304 (2018). agreement with ab initio calculations29. In the Dirac-fluid regime, 11. Gooth, J. et al. Thermal and electrical signatures of a hydrodynamic electron with an electron temperature of roughly 1,000 K, we obtain a thermal fluid in tungsten diphosphide. Nat. Commun. 9, 4093 (2018). Nature Nanotechnology | VOL 16 | November 2021 | 1195–1200 | www.nature.com/naturenanotechnology 1199 Letters Nature Nanotechnology 12. Berdyugin, A. I. et al. Measuring hall viscosity of graphene’s electron fluid. 25. Block, A. et al. Tracking ultrafast hot-electron diffusion in space and time by Science 364, 162–165 (2019). ultrafast thermomodulation microscopy. Sci. Adv. 5, eaav8965 (2019). 13. Sulpizio, J. A. et al. Visualizing Poiseuille flow of hydrodynamic electrons. 26. Gabor, N. M. et al. Hot carrier-assisted intrinsic photoresponse in graphene. Nature 576, 75–79 (2019). Science 334, 648–652 (2011). 14. Sheehy, D. E. & Schmalian, J. Quantum critical scaling in graphene. Phys. Rev. 27. Brida, D. et al. Ultrafast collinear scattering and carrier multiplication in Lett. 99, 226803 (2007). graphene. Nat. Commun. 4, 1987 (2013). 15. Xie, H. Y. & Foster, M. S. Transport coefficients of graphene: interplay of 28. Stauber, T. et al. Interacting electrons in graphene: Fermi velocity impurity scattering, Coulomb interaction, and optical phonons. Phys. Rev. B. renormalization and optical response. Phys. Rev. Lett. 118, 266801 (2017). 93, 195103 (2016). 29. Kim, T. Y., Park, C. H. & Marzari, N. The electronic thermal conductivity of 16. Lucas, A., Crossno, J., Fong, K. C., Kim, P. & Sachdev, S. Transport in graphene. Nano Lett. 16, 2439–2443 (2016). inhomogeneous quantum critical fluids and in the Dirac fluid in graphene. 30. Balandin, A. A. et al. Superior thermal conductivity of single-layer graphene. Phys. Rev. B. 93, 075426 (2016). Nano Lett. 8, 902–907 (2008). 17. Lucas, A. & Fong, K. C. Hydrodynamics of electrons in graphene. J. Phys. 31. Machida, Y., Matsumoto, N., Isono, T. & Behnia, K. Phonon hydrodynamics Condens. Matter 30, 053001 (2018). and ultrahigh–room-temperature thermal conductivity in thin graphite. 18. Zarenia, M., Smith, T. B., Principi, A. & Vignale, G. Breakdown of the Science 367, 309–312 (2020). Wiedemann–Franz law in AB-stacked bilayer graphene. Phys. Rev. B 99, Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in 161407 (2019). published maps and institutional affiliations. 19. Crossno, J. et al. Observation of the Dirac fluid and the breakdown of the Wiedemann–Franz law in graphene. Science 351, 1058–1061 (2016). Open Access This article is licensed under a Creative Commons 20. Ghahari, F. et al. Enhanced thermoelectric power in graphene: Attribution 4.0 International License, which permits use, sharing, adap- violation of the Mott relation by inelastic scattering. Phys. Rev. Lett. 116, tation, distribution and reproduction in any medium or format, as long 136802 (2016). as you give appropriate credit to the original author(s) and the source, provide a link to 21. Gallagher, P. et al. Quantum-critical conductivity of the Dirac fluid in the Creative Commons license, and indicate if changes were made. The images or other graphene. Science 364, 158–162 (2019). third party material in this article are included in the article’s Creative Commons license, 22. Ku, M. J. H. et al. Imaging viscous flow of the Dirac fluid in graphene. Nature unless indicated otherwise in a credit line to the material. If material is not included in 583, 537–541 (2020). the article’s Creative Commons license and your intended use is not permitted by statu- 23. Ruzicka, B. A. et al. Hot carrier diffusion in graphene. Phys. Rev. B. 82, tory regulation or exceeds the permitted use, you will need to obtain permission directly 195414 (2010). from the copyright holder. To view a copy of this license, visit http://creativecommons. 24. Zhu, T., Snaider, J. M., Yuan, L. & Huang, L. Ultrafast dynamic microscopy of org/licenses/by/4.0/. carrier and exciton transport. Annu. Rev. Phys. Chem. 70, 219–244 (2019). © The Author(s) 2021 1200 Nature Nanotechnology | VOL 16 | November 2021 | 1195–1200 | www.nature.com/naturenanotechnology Nature Nanotechnology Letters Methods with only one pulse at a time, analogous to the experimental difference-frequency Fabrication of split-gate thermoelectric device. The split-gate device with demodulation. Hall geometry consists of exfoliated, single layer graphene encapsulated by hBN, prepared using standard exfoliation and dry transfer techniques. The Quantifying the spatial spread. The following analysis is conducted both on the hBN-graphene-hBN stack is placed on a predefined split-gate structure made of experimental and the simulated data of ∆ITE(∆x, ∆t) for ‘symmetric experiments’ graphene, grown by chemical vapour deposition, where the gap between the two with optical pulses incident at a distance ∆x on each side of the pn-junction (Figs. 1 gates is roughly 100 nm, created via electron-beam lithography and reactive ion and 2). For each ∆t of the datasets ∆ITE(∆x, ∆t) we calculate the width of the etching (RIE). The top hBN and graphene are etched into a Hall bar shape with signal via the second moment, which for an ideal Gaussian profile is equal to the laser lithography and RIE, keeping the split-gate intact and not etching completely squared Gaussian width σ2. The second moment is calculated from the pixels ∆xi through the bottom hBN. Finally, the Ti/Au side contacts are created by a further (i = 1,…, N) via step of lithography, RIE and metal evaporation. The fabrication steps are shown in ∑ |Δx −Δx|2 ΔITE (Δxi ,Δt) Supplementary Fig. 6. < Δx2 > (Δt) = i ∑i ΔITE (Δxi , Δt) , with the mean i ∑ Δxi ΔITE (Δxi ,Δt) Spatiotemporal thermoelectric current microscopy setup. Our setup enables Δx = i∑ ΔITE (xi ,t) . i 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 We note that the minimum second moment at the focus <∆x2>focus of (the pn junction), while we vary the spatial displacement of our two laser pulses 0.56 µm2 comes from simulating the symmetric experiment, using as input the with respect to this junction and vary the temporal delay between the two measured Gaussian beam width at the focus σfocus2 = 0.14 µm2 (Supplementary ultrashort pulses. This means that we are following in space and time the diffusion Note 2). For the ‘asymmetric experiments’ with one optical pulse always incident of light-induced electronic heat from the location of light incidence to the pn on the pn junction (data for Fig. 3), we always consider the spatial profile only at junction. It is the thermoelectric effect at the pn junction, governed by the Seebeck time zero. Here we find that Gaussian fits with a background give the most reliable coefficient, that generates our observable signal, the thermoelectric current. We results. The entire set of data is shown in Supplementary Fig. 5. For each note that although the value of the Seebeck coefficient itself changes when changing dataset ∆ITE(∆x) or ∆ITE(∆y)(taken )at ∆t = 0, we perform Gaussian fits using EF, and when entering the hydrodynamic regime20, this only affects the magnitude 2 the function f(Δx) = a exp − Δx + b, where the Gaussian squared width of the thermoelectric current, rather than how electronic heat is diffusing outside 2σ 2 the pn junction, which is what we are following with our spatiotemporal technique. σ2 indicates the thermal spreading. Here, the minimum simulated Gaussian A sketch of the setup is shown in Supplementary Fig. 7. A Ti:sapphire oscillator widths are (σx2)focus = 0.34 µm2 and (σy2)focus = 0.44 µm2 (Supplementary Note 2). (886-nm centre wavelength, 76-MHz repetition rate), is split into two beam The experimentally obtained widths from this dataset as function of gate voltage paths. Both beams are modulated with optical choppers, at frequencies f1 = 741 and optical power are also shown in Supplementary Fig. 5, showing an increase and f2 = 529 Hz. The relative time delay between the two pulses is controlled with power, that is a larger Te, and an increase towards the Dirac point, that by a mechanical delay line. The spatial offset of one beam with respect to the is a smaller TF. We estimate the theoretical Gaussian widths in Fig. 3g using other is controlled with a mirror galvanometer, while the position of the sample σ 2calc = (σx2)focus + 2D ΔtIRF, where D are the calculated diffusivities. 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. Electrical measurements. We characterize our device electrically with four-probe We collect the thermoelectric (TE) photocurrent between the source and drain measurements (Extended Data Fig. 1), finding a charge mobility μ of 30,000– contacts through the graphene sheet on either side of the junction via lock-in 50,000 cm2 Vs−1, depending on carrier density. The measured mobilities correspond amplification. By demodulating the current signal at the difference frequency of to a momentum relaxation time 𝜏mr of 300–500 fs. These relaxation times are the two modulation frequencies, f2 − f1 = 211.7 Hz, we isolate the signal caused longer than the temporal resolution (the IRF) of our measurement technique, by the interaction of both heating sources, which we call the interacting heat ΔtIRF ≅ 200 fs, thus allowing us to probe our system before and after momentum current ∆ITE. The temporal resolution of the setup of 200 fs is determined in the relaxation occurs, that is in the non-diffusive and diffusive regime. We use these sample plane of the microscope (Supplementary Fig. 8). The spatial resolution measured charge mobilities to calculate the expected thermal diffusivity via the e ∂ne/h defined by our spot sizes is below 1 μm, whereas the accuracy with which we can Einstein relation33,34 μe/h = ne/h ∂EF De/h, where e is the elementary charge, EF is observe electronic heat spreading is determined by the signal-to-noise ratio, and is the Fermi energy and ne/h is the electron/hole carrier density. For highly doped E2 estimated to be below 100 nm. graphene ( EF ≫ kB T ) the carrier density expression ne/h = π h̄2Fv2 , leads to the F simple relation: De/h = E2eF μe/h. We note that we obtain the identical result by Estimating Fermi temperature controlled by gate voltage. During photocurrent calculating D from the ratio of the 2D thermal conductivity κe,2D and the electronic measurements, the gate voltage Ux is applied to one side (x is ‘A’) or the other heat capacity Ce and using the Wiedemann–Franz law: κ e,2D /σ = π 2 /3 · (kB /e)2 Te, side (x is ‘B’) side of the split gate. We always apply a symmetric voltage around where kB is the Boltzmann constant and e the elementary charge, together with the experimentally determined Dirac point voltage UDP DP x : UA = UA + ΔU and the conductivity σ = neμ and the following heat capacity for graphene (valid UB = UDP B − ΔU . The gate electrode and the graphene form a capacitor with 2πε k2 T for Te < TF): Ce = 3h̄F2 vB2 e . Given the measured mobilities, we expect thermal the dielectric hexagonal boron nitride (hBN), with a thickness of thBN = 70 nm, F diffusivities around 2,000 cm2 s−1 for our sample. and a relative permittivity of ϵhBN = 3.56. The carrier density n is calculated via n = ϵe0tϵhBN hBN ΔU , where ϵ0 is the vacuum permittivity. We calculate the Fermi energy Thermal diffusivity and conductivity of the Dirac fluid. We estimate the EF and the Fermi temperature TF via E2F = πħ2 v2F · n and TF = EkBF , where kB is the enhanced thermal diffusivity of the Dirac fluid by comparing the measured Boltzmann constant. 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 Estimating carrier temperature controlled by laser power. The thermoelectric for the symmetric scan (Fig. 2), and 29,000 and 39,000 cm2 s−1 for the x and photovoltage is assumed to be proportional to the time-averaged increase of the y directions of the asymmetric scan (Fig. 3), respectively, where <Δx2> is electronic temperature Te above the ambient temperature T0, as in ref. 32. The replaced with (σx2) and (σy2), respectively. The same calculation for a second sublinear dependence of the thermoelectric current ITE on optical power for the device (Supplementary Note 3 and Extended Data Fig. 2) gives a diffusivity of device under study here for illumination with a single pulsed laser (λ = 886 nm) 100,000 cm2 s−1. The 3D thermal conductivity κ3D of the Dirac fluid is calculated is shown in Supplementary Fig. 9. With a linear temperature scaling of the from the diffusivity D and the electronic heat capacity Ce, via κ3D = DCe/d, electronic heat capacity for graphene away from the Dirac point, Ce(T) = γT, we where d is the thickness of graphene, 0.3 nm. For the Dirac fluid, we have Q0 +ΔQ Te Te > TF, and therefore use the ‘undoped’ electronic heat capacity35 π18(h̄v ζ (3) 3 2 F) 2 kB Te , integrate the heat energy per unit area dQ = CedT, that is, ∫ dQ = ∫ γTdT . where ζ (3) ≈ 1.202. With the above estimate D = 35,000–70,000 cm2 s−1 Q0 T0 and Te = 1,000 K, we obtain the 3D thermal conductivity κ3D = 18,000– With the incident power P proportional to the absorbed heat energy per unit area 40,000 W mK−1. This corresponds to a 2D κ2D > 5 μW K−1. This value is orders of ∆Q, we √find that the peak Te as a function of the laser power P scales as in ref. , 32 magnitude larger than the value found in ref. 19. The reason for this is that our Te = 2 T20 + bP . Here, the parameter b is defined via bP = 2∆Q/γ, and is used to electron temperature is more than ten times higher, and therefore the electronic convert incident power to peak electron temperature (Supplementary Fig. 9). 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 Simulation of the experiment. A detailed description of the simulation can be Dirac-fluid regime with its diverging thermal diffusivity. found in Supplementary Note 1 and Supplementary Fig. 10. In brief, we calculate the spatiotemporal evolution of electronic heat generated by the two optical Dirac-fluid crossover temperature. Following the treatment in ref. 14, we find the pulses in the graphene sheet via the heat equation with a finite difference method. crossover temperature from Fermi liquid to Dirac fluid, as a function of Fermi We define Gaussian heating pulses and calculate their temperature rise via the temperature as experimentally measured non-linear power scaling. We extract the differential TE T0 ( ( )) current contribution as a function of ∆x and ∆t by the difference of the heating at Tcross (TF ) = TF 1 + λ ln , the pn junction region in the presence of both pulses with respect to simulations TF Nature Nanotechnology | www.nature.com/naturenanotechnology Letters Nature Nanotechnology where λ = e2 /16ϵ0 ϵr vF h̄ ≈ 0.55/ϵr for graphene √ with the dielectric environment (N.F.v.H.) and the CERCA program (ICN2 and ICFO); Spanish MINECO under ϵr ≈ 3.56 for hBN. The temperature T0 = 23h̄v F 3/4 k a π B 0 ≈ 8.4 × 104 K , with the grant nos. SEV-2017-0706 (ICN2) and CEX-2019-000910-S (ICFO); the International −10 inter-atomic distance a0 = 1.42 × 10 m. The resulting crossover temperature PhD fellowship program ‘la Caixa’ (A.B.); Leverhulme Trust grant no. RPG-2019-363 is shown in Fig. 1b. We note that the relatively high refractive index of the hBN (A.P.); the Elemental Strategy Initiative conducted by the MEXT, Japan, grant no. encapsulant makes the Dirac fluid more easily accessible, as it lowers the crossover JPMXP0112101001 (K.W. and T.T.), JSPS KAKENHI grant no. JP20H00354 (K.W. and temperature compared to vacuum, by a factor of about two for the range of TF T.T.) and the CREST (grant no. JPMJCR15F3) and JST (K.W. and T.T.); Fundació Privada values studied here. Cellex (ICFO) and Fundació Mir-Puig (ICFO). Data availability Author contributions The data that support the findings of this study are available from the K.-J.T. conceived and supervised the project. A.B. developed the experimental setup corresponding author on reasonable request. under supervision of M.L. and N.F.v.H. and input from K.-J.T. N.C.H.H. performed sample fabrication, with material input from K.W. and T.T., under the supervision of F.H.L.K. and K.-J.T. A.B. performed the measurements under the supervision of K.-J.T. References A.P. performed the theoretical hydrodynamic transport calculations. A.B. and K.-J.T. 32. Tielrooij, K. J. et al. Out-of-plane heat transfer in van der Waals stacks through interpreted and analysed the data, with input from M.L., A.W.C. and N.F.v.H. A.B. electron-hyperbolic phonon coupling. Nat. Nanotechnol. 13, 41–46 (2018). developed the model that simulated the experiment with input from K.-J.T. A.W.C. and 33. Zebrev, G. I. I. in Physics and Applications of Graphene—Theory S.R. developed the ballistic transport model. A.B. and K.-J.T. wrote the paper, with input (ed. Mikhailov, S.) 475–498 (IntechOpen, 2011). from all authors. 34. Rengel, R. & Martín, M. J. Diffusion coefficient, correlation function, and power spectral density of velocity fluctuations in monolayer graphene. J. Appl. Phys. 114, 143702 (2013). Competing interests 35. Lui, C. H., Mak, K. F., Shan, J. & Heinz, T. F. Ultrafast photoluminescence The authors declare no competing interests. from graphene. Phys. Rev. Lett. 105, 127404 (2010). Additional information Acknowledgements Extended data are available for this paper at https://doi.org/10.1038/ We thank M. Polini and P. Piskunow for fruitful discussions, and H. Agarwal and K. s41565-021-00957-6. Soundarapandian for help with sample fabrication. We acknowledge the following Supplementary information The online version contains supplementary material funding sources: European Union’s Horizon 2020 research and innovation programme available at https://doi.org/10.1038/s41565-021-00957-6. under grant nos. 804349 (K.-J.T.), 873028 (A.P.), 785219 (F.H.L.K. and S.R.), 881603 Correspondence and requests for materials should be addressed to K.-J.T. (F.H.L.K. and S.R.) and 670949 (N.F.v.H.); Spanish MCIU/AEI under grant nos. RYC- 2017-22330 (K.-J.T.), PID2019-111673GB-I00 (K.-J.T.), BES-2016-078727 (M.L.), Peer review information Nature Nanotechnology thanks Giulio Cerullo and the other, RTI2018-099957-J-I00 (M.L.) and PGC2018-096875-B-I00 (M.L. and N.F.v.H.); the anonymous, reviewer(s) for their contribution to the peer review of this work. Government of Catalonia under grant nos. SGR1656 (F.H.L.K.) and 2017SGR1369 Reprints and permissions information is available at www.nature.com/reprints. Nature Nanotechnology | www.nature.com/naturenanotechnology Nature Nanotechnology Letters Extended Data Fig. 1 | Electrical mobility measurement. Momentum relaxation time from four-probe measurements and corresponding (calculated) heat diffusivity (solid line). Four-probe measurements were performed by applying 1 V to a MΩ series resistor, such that a current of 1 μA flows between the two outer contacts of the device (see inset). We then measure the voltage drop between two lateral contacts of the graphene device. This yields the μEF sheet conductance σ as a function of gate voltage. We then use σ = neμ, in order to extract the mobility µ and use τ mr = to obtain the momentum ev2F relaxation time. Here, n is the carrier density, EF is the Fermi energy, vF ≈ 10 m/s is the Fermi velocity, and e is the elementary charge. We measured up to 6 EF = 150 meV, and extrapolated the data to higher Fermi energies. The three symbols indicate the Fermi energy and predicted diffusivity corresponding to the diffusivity measurements in Fig. 2e. Nature Nanotechnology | www.nature.com/naturenanotechnology Letters Nature Nanotechnology Extended Data Fig. 2 | See next page for caption. Nature Nanotechnology | www.nature.com/naturenanotechnology Nature Nanotechnology Letters Extended Data Fig. 2 | Spatiotemporal results from a second device. a, Asymmetric spatiotemporal ∆ITE maps for three different gate voltages. b, Extracted width as a function of ∆t. A lower Fermi level leads to a higher time-zero width, in accordance with hydrodynamic transport, as presented for the main device in the manuscript. (c) ∆ITE maps, taken at ∆t = 0, as a function of beam offset (∆x, ∆y), as well as sample height (z). (d) extracted line profiles for the two dimensions. e, Resulting signal width σ2 for both dimensions as extracted from Gaussian fits at each z-position. The same measurements are presented in Supplementary Fig. 4 for the main device. These experiments were performed with two beams of wavelength 443 nm and 886 nm, respectively. Nature Nanotechnology | www.nature.com/naturenanotechnology Letters Nature Nanotechnology Extended Data Fig. 3 | Spatiotemporal results from a third device. a, Microscope image of the device. This sample has split-gates made from graphite, with a gap that is 200 nm. It then has a 30 nm thick layer of SiO2, and then we transferred a graphene flake grown by chemical vapour deposition (CVD) on top of the split-gate structure using an hBN flake. b, Gate-dependent current measurement, which gives an estimated mobility of this half-encapsulated CVD graphene sample of around 8,500 cm2/Vs (solid lines). c-d, Asymmetric spatiotemporal ∆ITE maps at time zero for different gate voltages and laser powers without (c) and with normalization (d). e-f, Comparison between first and third device. 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 power and gate voltage, for both the first device (hBN-encapsulated with high mobility, presented in the main text, blue-purple colours, ‘hBN’) and the third device (on SiO2 with low mobility, yellow-red colors, ‘SiO2’). The low-mobility sample with shorter hydrodynamic time window shows systematically less heat spreading around time zero, in agreement with our picture of hydrodynamic heat spreading during the hydrodynamic time window. Nature Nanotechnology | www.nature.com/naturenanotechnology Nature Nanotechnology Letters Extended Data Fig. 4 | Ballistic spreading simulation. a, Time-dependent output distributions for quantum mechanical calculations for a single electron (top row) and an ensemble of independent electrons (bottom row). b, Resulting width σball2 for ballistic transport for the Monte Carlo method with varying Fermi velocities, as well as the quantum calculation for an ensemble of independent electrons. Both calculations essentially agree and the spread within 0.25 ps leads to final widths of below 0.25 µm2 for realistic values of the Fermi velocity. Nature Nanotechnology | www.nature.com/naturenanotechnology