Skip to main content

Thank you for visiting You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.

Thermal dynamics and electronic temperature waves in layered correlated materials


Understanding the mechanism of heat transfer in nanoscale devices remains one of the greatest intellectual challenges in the field of thermal dynamics, by far the most relevant under an applicative standpoint. When thermal dynamics is confined to the nanoscale, the characteristic timescales become ultrafast, engendering the failure of the common description of energy propagation and paving the way to unconventional phenomena such as wave-like temperature propagation. Here, we explore layered strongly correlated materials as a platform to identify and control unconventional electronic heat transfer phenomena. We demonstrate that these systems can be tailored to sustain a wide spectrum of electronic heat transport regimes, ranging from ballistic, to hydrodynamic all the way to diffusive. Within the hydrodynamic regime, wave-like temperature oscillations are predicted up to room temperature. The interaction strength can be exploited as a knob to control the dynamics of temperature waves as well as the onset of different thermal transport regimes.


The capability to access ultrafast thermal dynamics recently gave access to striking phenomena that take place in materials at the nanoscale before complete local energy equilibration among heat carriers is achieved1,2,3,4,5. For instance, non-Fourier heat transport regimes have been reported for hot spots dimensions inferior to the phonon mean free-path6,7,8, in which energy is ballistically carried point to point, or have been engineered via nano-patterning of dielectric substrates9,10,11. As a consequence of the existence of two non-thermal populations, wave-like thermal transport, often referred to as second sound12,13, has been predicted in graphene, both in the frame of microscopic14,15,16,17 and macroscopic models18. Temperature wave-like phenomena have been recently observed at high temperatures in graphene19 and 2D materials20 on sub-nanosecond timescales and scheme for their coherent control have been proposed21. So far most of the effort has been devoted to phononic non-Fourier heat transport17,19,22,23,24,25,26, where, only recently, a theoretical framework, covering on equal footing Fourier diffusion, hydrodynamic propagation, and all regimes in between, has been proposed27. On the contrary, electronic non-Fourier heat transport remains relatively unexplored18,20,28,29.

In this work, we propose layered correlated materials (LCM) as the ideal platform to access the entire spectrum of unconventional electronic heat transport regimes. In quantum correlated materials, the strong electronic interactions give rise to emerging many-body properties, such as collective and decoupled diffusion of energy and charge30,31. Tuning the interaction strength thus opens the possibility to investigate transport regimes with no counterpart in conventional weakly-interacting materials32,33.

We consider an impulsive excitation on the surface of a LCM characterized by a strong local Coulomb interaction U (see Fig. 1). The interaction U can drive fast local thermalization processes leading to the rapid build up of a hot intra-layer electronic temperature before relaxation via slower scattering paths takes place. At the same time, the interaction leads to heavier quasiparticles with enhanced effective mass m* and a reduced kinetic energy. As a consequence, energy propagation across the layers is expected to slow down for increasing U. Overall the interaction U, together with the anisotropy of the inter-layer and intra-layer hopping terms, may thus act as a tuning parameter to control the relative inter- and intra-layer energy exchange processes in LCM. Eventually, as the interaction increases, the two processes can effectively decouple, thus allowing for unconventional electronic heat transport regimes to occur on the ultrashort space and timescales.

Fig. 1: Setup.
figure 1

Cartoon of the layered correlated material impulsively excited on the top surface by ultrafast light pulses. We assume that the excitation drives a fast thermalization of the electronic population establishing an electronic temperature Thot on the topmost layers of the sample. The right panels display the electronic distributions at different depths after the equilibration step (see text).

We investigate the possibility for unconventional heat transport regimes by focusing on the impulsive thermal dynamics of the layered single-band Hubbard model, which represents a general framework for understanding the effects of electronic interactions in a large family of correlated materials. We show that on sub-picosecond timescales the electronic heat transfer is initially characterized by ballistic wavefront propagation, followed by an hydrodynamic regime, which eventually evolves into conventional Fourier heat transfer on longer timescales. In the hydrodynamic regime, we predict that LCM may sustain temperature wave oscillations at THz frequencies and up to ambient temperature.

In order to contextualize the present concepts within the frame of real systems and to connect with the realm of technologically relevant materials, we focus on the correlated metal SrVO3 (SVO). SVO is a paradigmatic representative of the wider class of correlated transition metal oxides (TMOs) and it has been proposed as a platform for a wealth of potential technological applications ranging from ideal electrode materials34, to Mott transistors35 and transparent conductors36. We argue that the degree of correlation of SVO, as measured by the interaction strength, is such that ballistic transport first, and wave-like thermal transport afterwords, are accessible on the sub-picosecond timescale. Our results, together with the possibility of heterostructuring TMO to atomic layer accuracy, promote these materials to ideal building blocks for nanothermal device architectures based on non-Fourier heat transport.

The present work rationalizes the microscopic interactions underlying unconventional electronic heat transfer phenomena in LCM. Our findings enlarge the functionalities of quantum materials32,33 to the realm of nanoscale heat transport37, beyond the case of radiative energy transfer38,39,40.

Results and discussion

Model and observables

The Hamiltonian of the layered Hubbard model reads:

$$H=\mathop{\sum }\limits_{n=1}^{L}{h}_{n}+\mathop{\sum }\limits_{n=1}^{L-1}{\tau }_{n,n+1}$$


$${h}_{n}=\mathop{\sum}\limits_{ < i,j > \sigma }{t}_{\parallel }{c}_{in\sigma }^{{{{\dagger}}} }{c}_{jn\sigma }^{}+U\mathop{\sum}\limits_{i}{n}_{in\uparrow }{n}_{in\downarrow }$$


$${\tau }_{n,n+1}=\mathop{\sum}\limits_{\sigma }{t}_{\perp }{c}_{in\sigma }^{{{{\dagger}}} }{c}_{in+1\sigma }^{}+h.c.$$

where \({c}_{in\sigma }^{{{{\dagger}}} }\) is a fermionic creation operator for an electron with spin σ at the site i belonging to the layer indexed by n, which ranges from 0 to L. t and t represent, respectively, the intra- and inter-plane hopping amplitudes. The sum of the in-plane hopping term runs over pairs of nearest neighbouring sites and we introduce the number operator \({n}_{in\sigma }={c}_{in\sigma }^{{{{\dagger}}} }{c}_{in\sigma }^{}\). We assume in-plane translational invariance so that we can introduce an in-plane momentum k = (kx, ky) and recast \({h}_{n}={\sum }_{{{{{{{{\bf{k}}}}}}}}\sigma }\epsilon ({{{{{{{\bf{k}}}}}}}}){c}_{{{{{{{{\bf{k}}}}}}}}n\sigma }^{{{{\dagger}}} }{c}_{{{{{{{{\bf{k}}}}}}}}n\sigma }^{}+U{\sum }_{i}{n}_{in\uparrow }{n}_{in\downarrow }\) with \(\epsilon ({{{{{{{\bf{k}}}}}}}})=-2{t}_{\parallel }(\cos ({k}_{x}a)+\cos ({k}_{y}a))\) and a the lattice spacing. We fix the chemical potential in order to have an average occupation of one electron per site (half-filling) corresponding to the perfect particle-hole symmetric case. As a consequence, the total number of electrons per layer is conserved during the dynamics. This choice allows us to model energy transport in the absence of mass and charge transport. We consider the model parameters t = t = 60 meV and U = 0.65 eV, which correspond to an interaction-driven mass renormalization m/m* 0.3, and a lattice spacing a = 5 Å. This value of the effective mass renormalization is consistent with experimental estimates for SrVO3.

We study the non-equilibrium thermal dynamics in the frame of model (1) by means of a time-dependent variational approach based on the generalized Gutzwiller approximation for layered systems41,42 (see Methods). This approach provides a versatile tool for describing in a non-perturbative way the dynamics in the Hubbard model, which is governed by the interplay between the hopping terms t, t and the local Coulomb interaction U. In this framework, the dynamics is described by using a variational density matrix \(\hat{\rho }(t)\) with the structure43,44 \(\hat{\rho }(t)\equiv {{{{{{{\mathcal{P}}}}}}}}(t){\hat{\rho }}_{* }(t){{{{{{{{\mathcal{P}}}}}}}}}^{{{{\dagger}}} }(t)\), where \({\hat{\rho }}_{* }\) is a density matrix which describes the dynamics of quasiparticles through an effective non-interacting Hamiltonian, and \({{{{{{{\mathcal{P}}}}}}}}={\prod }_{i}{{{{{{{{\mathcal{P}}}}}}}}}_{i}\) are projectors onto the local Hilbert space at site i which control the relative weights of the local many-body configurations. The mutual feedback between the dynamics of the localized degrees of freedom and the low-energy quasiparticle excitations results in an effective mass renormalization m* = m/Z, which is controlled by the interaction U through the quasiparticle weight Z(U). In the non-interacting limit Z(U = 0) = 1, whereas at finite interaction Z < 1 and decreases as a function of U. Eventually, for a critical interaction strength, Uc, the system undergoes a metal-to-insulator Mott transition, corresponding to a vanishing quasiparticle weight, i.e. Z(Uc) → 0. In this regime, quasiparticle excitations are completely suppressed and the dynamics becomes dominated by high-energy incoherent excitations at energies ~ U45. In this work we focus on the thermal dynamics of hot quasiparticles in the correlated metal regime, where Z is finite but significantly smaller than one.

The thermal dynamics is triggered by a sudden increase of the electronic temperature localized within the first few surface layers of the LCM as can be achieved, for instance, by excitation with a femtosecond light pulse28. We mimic the impulsive excitation by considering the two-step non-equilibrium protocol (see Methods): first, we initialize the system in a non-equilibrium state characterized by two-electronic populations at two different temperatures. In the second step, we consider the unitary evolution of the initialized state, as regulated by the Hamiltonian (1). The initialization steps consists in a preliminary dynamics starting from the zero temperature state in which each layer n is independently coupled with baths of non-interacting electrons at temperatures Tbath(n). All the layers are at a base temperature Tbath(n) = Tc0 except the five topmost which are set at Tbath(n) = ThTc0, for n = 1, …, 5, a temperature characterizing the hot electrons. In the notation Tc0, the subscript “c” stands for cold and “0” indicates the instant right after the impulsive excitation, whereas the subscript “h” in Th stands for hot. Throughout the paper, the temperatures used in the equilibration step are Tc0 35 K for the cold temperature in the bulk and Th = 10 × Tc0 for the hot temperature in the topmost layers.

The end of the equilibration step defines the instant of time t = 0. At t = 0, we initialize the density matrix with the density matrix obtained at the end of the equilibration step \(\hat{\rho }(t=0)={\hat{\rho }}_{{{{{{{{\rm{equ}}}}}}}}}\). We switch off all the couplings to the baths and let the initialized density matrix unitary evolve for t > 0. We refer the reader to the Methods section for further details on the non-equilibrium protocol.

At positive times, we study thermal transport by tracking the time evolution of the electronic temperature and of the heat flux at each layer. The electronic temperature is extracted from the time and layer-dependent quasiparticle non-equilibrium distribution functions \({N}_{n}^{{{{{{{{\rm{neq}}}}}}}}}(\epsilon ,t)\) (see Methods). For any instant of time t > 0, and for any layer index n, we find upon fitting that \({N}_{n}^{{{{{{{{\rm{neq}}}}}}}}}(\epsilon ,t)\) can be expressed as a superposition of two equilibrium Fermi distributions: (i) a hot distribution at the temperature Th, fixed by the initial equilibration temperature in the five topmost layers, and of weight ρhot(n, t); (ii) a cold distribution characterized by a time- and layer-dependent temperature T(n, t) and of weight ρcold(n, t) ≡ 1 − ρhot(n, t)

$${N}_{n}^{{{{{{{{\rm{neq}}}}}}}}}(\epsilon ,t)={f}_{{{{{{{{\rm{hot}}}}}}}}}{\rho }_{{{{{{{{\rm{hot}}}}}}}}}(n,t)+{f}_{{{{{{{{\rm{cold}}}}}}}}}(n,t)\left[1-{\rho }_{{{{{{{{\rm{hot}}}}}}}}}(n,t)\right],$$

with 0 ≤ ρhot(n, t) ≤ 1, fhot = f(ϵ, Th) and fcold(n, t) = f(ϵ, T(n, t)). For any instant of time t and layer index n, we fit Nneq(ϵ, t), computed via Eq. (19), with the expression given by Eq. (4), ρhot(n, t) and T(n, t) being the only two fitting parameters. While the temperature of the hot electrons is fixed at Th, the temperature of the remaining 1 − ρhot(n, t) fraction of electrons in the “cold” state T(n, t) changes in time. In Fig. 2 we show the decomposition, Eq. (4), for the layer n = 15 and two instants of time at t = 0.32 ps and t = 0.81 ps for which the non-equilibrium distribution functions are described by different populations of hot electrons, and different cold temperatures T(n, t).

Fig. 2: Non-equilibrium distribution functions decomposition.
figure 2

Non-equilibrium distribution functions, Eq. (4), decomposed in terms of the cold and hot Fermi-distribution functions. Data correspond to the non-equilibrium distribution functions for the layer n = 15 at two instants of time t = 0.32 ps (left panels) and t = 0.81 ps (right panels). Top panels a, b: ρcoldfcold (blue lines) and sum of the contributions ρcoldfcold + ρhotfhot (red line) compared with the non-equilibrium distribution functions \({N}_{n}^{{{{{{{{\rm{neq}}}}}}}}}\) (dashed black line), obtained by the numerical simulation. For simplicity we removed the layer and time index from the labels. The blue/red shaded regions helps to visualize the cold and hot contributions to the non-equilibrium distribution functions. Bottom panels c, d: ΔNcold, difference between the cold distribution function at the time t, and the initial one at temperature Tc0.

Remarkably, the effective electronic temperatures at t = 0.32 ps and t = 0.81 ps extracted from the fit of Eq. (4) to \({N}_{n}^{neq}(\epsilon ,t)\) result to be, respectively, larger and smaller than the initial Tc0. In order to demonstrate that these oscillations of T(n, t) are not an artifact related to the fitting procedure, we report in Fig. 2c, d direct evidence of temperature oscillations in the calculated distribution function. More specifically, we subtracted the contribution of the hot electrons (ρhot(n, t)fhot) from the total distribution function \({N}_{n}^{neq}(\epsilon ,t)\), as calculated from the unitary dynamics. The normalized cold distribution can be thus expressed as \({N}_{n}^{{{{{{\mathrm{cold}}}}}}}(\epsilon ,t)=[{N}_{n}^{neq}(\epsilon ,t)-{\rho }_{{{{{\mathrm{{hot}}}}}}}(n,t){f}_{{{{{{\mathrm{hot}}}}}}}]/[1-{\rho }_{{{{{{\mathrm{hot}}}}}}}(n,t)]\). The difference between \({N}_{n}^{{{{{{\mathrm{cold}}}}}}}(\epsilon ,t)\) and the initial Fermi distribution, i.e. \({{\Delta }}{N}_{n}^{{{{{{\mathrm{cold}}}}}}}={N}_{n}^{{{{{{\mathrm{cold}}}}}}}(\epsilon ,t)-f(\epsilon ,{T}_{c0})\), clearly shows a change of occupation in the vicinity of the Fermi level which resembles a broadening (narrowing) of a Fermi distribution at t = 0.32 (0.81) ps. These data directly demonstrate that, during the dynamics, the electronic distribution changes in a complex oscillatory way, which can be accounted for by the two-temperature effective model described by Eq. (4). This fact allows for a clear physical interpretation of the transient propagation of energy and a practical definition of a local time-dependent electronic temperature. Initially, the perturbation creates a population of hot electrons described by ρhot(n, t), which propagates across the layers. The interaction between the cold electrons on each layer and the hot electrons propagating through the system creates a modulation of the temperature of the cold electronic population, encoded in the spatio-temporal evolution of T(n, t). In the rest of the manuscript, we will refer to this quantity as the cold electronic temperature, bearing in mind that T(n, t = 0) = Tn, with Tn the bath temperatures defined in the initialization step.

We extract the heat flux by considering the time evolution of the energy density defined for each layer as

$${E}_{n}(t)\equiv {{{{{{{\rm{Tr}}}}}}}}[\hat{\rho }(t)({h}_{n}+{\tau }_{n,n+1})]$$

where \(\hat{\rho }(t)\) represents the time-evolved Gutzwiller density matrix. At positive times the density matrix evolves under unitary dynamics, so that energy is conserved. As a consequence, there exists a continuity equation that relates the time evolution of the local energy density En to the divergence of the current associated with the transport of energy, namely the heat flux. Therefore, by starting from the dynamics of the energy density, we compute the heat flux qn at layer n, and along the z-direction perpendicular to the planes, by applying the continuity equation

$$\frac{\partial {q}_{n}}{\partial z}+\frac{\partial {E}_{n}}{\partial t}=0$$

where the discrete spatial derivative is defined with respect to the inter-layer distance a, \(\frac{\partial {q}_{n}}{\partial z}=({q}_{n+1}-{q}_{n})/a\).

Ultrafast thermal dynamics

We now show how this model offers the possibility to access different regimes of non-conventional heat transport on the sub-picosecond timescale. Each regime will be then separately discussed in the following. Figure 3a reports the results for the time evolution of the hot population weight and the local electronic relative effective temperature variation, ΔT/Tc0(n, t) with ΔT = T(n, t) − Tc0, recorded on layer n = 15, which we take as representative of the inner region of the slab. For times 0 < t 150 fs, both ρhot(15, t) and T(15, t) remain fixed to the equilibrium values ρhot = 0 and T = Tc0. At t ~ 150 fs the perturbation reaches the n = 15 layer and the dynamics that follows can be neatly divided in three steps.

Fig. 3: Sub-picosecond thermal dynamics.
figure 3

a Dynamics of the hot-electron population (top) and relative temperature variation of the “cold” electronic population (bottom) recorded on layer n = 15. b Layer profiles of the temperature of the “cold” electronic population (red, left axis) and of the heat flux (blue, right axis) at different instant of times. The blurred yellow band highlights the wave packet of temperature oscillations that follows the ballistic front. c Dynamics of the heat flux (blue, right axis) and temperature gradient of the “cold” electronic population (brown, left axis) on the layer n = 15. Arrows indicates the three regimes of thermal transport discussed in the main text. Panels a and c share the same x-axis. For simplicity the labels are shown only in the x-axis of panel c.

  1. (i)

    In the time window 150–400 fs the dynamics is characterized by a significant increase of ρhot(15, t) highlighting the arrival of the propagating hot-electron population. On this timescale, the electronic relative temperature variation remains limited. This is indicative of a ballistic regime of energy transport in which the energy flows without inducing any heating in the underlying quasi-equilibrium distribution.

  2. (ii)

    For t 400 fs the hot-electron population displays a sharp drop and, concomitantly, we observe the activation of a fast oscillatory dynamics in the electronic temperature of the cold electrons. Initially the oscillations are centered around a value higher than the initial equilibrium temperature Tc0 indicating that the transit of the ballistic front of hot electrons induced the heating of the population of cold electrons on the layer.

  3. (iii)

    Eventually the system equilibrates for t 0.9 ps with the residual damped temperature oscillations converging to Tc0.

We gain further insight into the thermal dynamics by comparing the dynamics of the local cold electronic temperature with the heat flux qn(t) at layer n. Figure 3b reports the spatial profiles of the heat flux (right axis, blue trace) and of the local electronic temperature (left axis, red trace) at fixed instants of time. The broad feature at the forefront of the heat-flux profile indicates the propagation of a ballistic energy front accompanied by a small and more localized perturbation of the electronic temperature. At the back front of the ballistic heat flux, as indicated by the blurred yellow band in Fig. 3b, we observe the formation of a sharp sinusoidal feature in the spatial profile of the temperature. In the time domain, this sharp feature marks the separation between the first two dynamical regimes of the local temperature observed in panel a for the layer n = 15. The presence of this pronounced oscillation of the temperature spatial profile is accompanied by weaker temperature oscillations with smaller spatial periodicity in the layers behind the ballistic front.

To fully characterize the thermal dynamics regimes occurring after the ballistic front has transited, we further compare the dynamics of the heat flux with that of the temperature gradient T(t) = (Tn+1(t) − Tn(t))/a perpendicular to the layers. These quantities are shown in Fig. 3c for the n = 15 layer. In the time window 0.15–0.4 ps, the ballistic regime shows up as a sharp increase of the heat flux with no sizeable effect on the temperature gradient. In correspondence of the end of the ballistic regime, i.e. the sharp drop of the heat current, an oscillatory dynamics is activated for the temperature gradient. The oscillatory dynamics of T(t) is maintained in the 0.4–0.9 ps time window, along with a residual positive heat current on the layer. At t 0.9 ps the heat current displays damped oscillations centred around zero indicating the recovery of local thermal equilibrium. Remarkably, the equilibration is characterized by the synchronization between the dynamics of the temperature gradient and the heat flux. In this regime, we can define an instantaneous proportionality between the heat flux and temperature gradient, i.e. qn(t) T(t), indicating that the heat transfer process is well described by a Fourier-like heat transfer law.

At intermediate times (0.4 < t < 0.9 ps), before Fourier-like transport sets in, there is a residual positive flow of the heat current with an oscillatory dynamics of −T(t) that is not simply proportional to that of qn(t). This fact reveals the presence of a new heat transport regime, which bridges the ballistic regime established at the arrival of the perturbation (0.15 < t < 0.4 ps) and the Fourier-like transport setting in at long times after the perturbation has transited (t > 0.9 ps). This intermediate regime is characterized by a residual population of hot electrons on the layer and by an oscillatory dynamics of the temperature of the cold electron population. We identify this regime as a hydrodynamic transport of heat sustained by the exchange of energy between the two sub-populations of hot and cold electrons. By comparing the dynamics on the single layer (Fig. 3a, c) with the layer profiles at different times (Fig. 3b), we can observe that the emergence of the hydrodynamic regime coincides with transit of the sharp sinusoidal feature in the spatial profile of the temperature at the trailing edge of the heat-flux ballistic front. As it will be further discussed, this feature can be considered as a temperature wave-packet propagating through the system.

Summarizing, the sub-picosecond thermal dynamics of electrons displays three subsequent regimes of heat transport: (i) the ballistic propagation of energy at the front of the perturbation; (ii) the hydrodynamic regime at the trailing edge of the ballistic front. The former is characterized by a wave-like propagation of the electronic temperature; (iii) a Fourier-like heat transport driving the recovery of thermal equilibrium. The time and space extension of the three regimes are indicated by the arrows in the plots of the dynamics at fixed layer index (see Fig. 3c) and of spatial profiles at fixed time (see Fig. 3b). In the remaining of the paper we analyse in detail the different regimes and discuss the possibility to control their onset in layered correlated materials.

Ballistic energy transport

We now address the possibility of controlling the initial ballistic energy transport by tuning the microscopic parameters entering in the Hubbard model (1). In the ballistic regime, the energy is mostly carried by the population of hot electrons at temperature Th. The energy propagates through hopping processes of the hot electrons excited in the first layers. Layered correlated materials thus offer two complementary ways to control the inter-layer coupling and, in turn, the velocity of propagation of the ballistic front, namely tuning either the anisotropy of the system, t/t, or the strength of the interaction, U. The increase of the latter drives a reduction of the quasiparticle weight Z, which leads to a larger effective mass for the inter-layer motion and a smaller effective hopping, \({t}_{\perp }^{* }=Z{t}_{\perp }\).

We show these effects in Fig. 4 where we report the spatio-temporal dynamics of the hot-electron population ρhot(n, t) obtained for different values of anisotropy (horizontal gray arrow) and relative interaction strength (vertical red arrow). Increasing either one, the propagation velocity of the wavefront is diminished. In the inset we plot the velocity of ballistic propagation vb as a function of U for t/t = 1. vb is defined as the slope of the white dashed line in Fig. 4. The correlation-induced renormalization of \({t}_{\perp }^{* }\) strongly suppresses the energy propagation along the z-direction. For the sake of applications, we note that, in nanosystems with sizes of the order of the ballistic mean free-path, the thermal conductivity becomes a size-dependent property46,47,48,49. Nanoengineering of LCM, combined with proper tuning of U and t/t, thus offers a mean to control, on the picosecond timescale, the velocity of ballistic heat pulses and, therefore, the thermal conductivity of nanodevices.

Fig. 4: Control of ballistic energy propagation.
figure 4

The four panels matrix displays the ballistic dynamics of the hot-electron population for varying values of the correlation strength U and anisotropy t/t (t = 60 meV). The color scale represents the amplitude of ρhot(n, t) (yellow: maximum; black: minimum). The inset displays the speed vb of the ballistic wavefront for different values of U at t = t. The ballistic wavefronts are highlighted by the dashed lines in the main panels.

Hydrodynamic energy transport: temperature waves

The results reported in Fig. 3 demonstrate that a purely electronic hydrodynamic transport regime can be achieved in our correlated system on much faster timescales than the more conventional phononic counterpart13,14,19. Similarly to the phononic case, this hydrodynamic regime manifests itself by a wave-like propagation of temperature oscillations, which emerge after the ballistic front has transited (see arrows in Fig. 3c). In this section, we will quantitatively describe the characteristics of temperature wave-like propagation, as it emerges from our microscopic model, and compare with a macroscopic model for the description of the hydrodynamic regime characterized by the emergence of electronic temperature waves.

In the microscopic model, we characterize the hydrodynamic regime by tracking the position of the minimum of the wave packet \({X}_{\min }\). We observe that \({X}_{\min }\) linearly increases in time (See the inset of Fig. 5a), allowing an estimate of the wave-packet group velocity from the simple relation \({X}_{\min }={v}_{g}t\), with vg ~ 30 nm/ps of the same order of magnitude of the ballistic energy wavefront velocity. A similar result is obtained when tracking the time-dependent maximum of the wave packet, \({X}_{\max }\). This result suggests that we can approximately describe the wave packet as a superposition of weakly dispersive waves with frequencies νk = vgk/2π.

Fig. 5: Spectral analysis in k-space of the electronic temperature waves.
figure 5

Top panel: spatial Fourier transform of the layer profile of the electronic temperature at three different times, t = 0.505 ps (diamonds), t = 0.735 ps (squares) and t = 0.945 ps (circles). The inset shows the position of the minimum of the temperature oscillation, indicated in the bottom panel, as a function of time. The vertical dashed line indicated by k* shows the average of the position of the three peaks. Bottom panel: Portions of the temperature profiles at different times used to compute the discrete Fourier transform. The arrows indicate the positions Xmin and Xmax.

In order to identify the barycentric wavevector of the propagating wave packet, in the top panel of Fig. 5 we report the spatial Fourier transform of the electronic temperature profile in the spatial window where the propagating packet is present, as highlighted in the three curves of the bottom panel, which correspond to three different times, t = 0.505 ps, t = 0.735 ps and t = 0.945 ps. The small number of layers included in the Fourier window produces spectrally broaden peaks with the maximum occurring at slightly different k values for different times. We estimate the peak wavevector by taking the average of the three peaks observed at the three chosen times, obtaining k* ~ 2.2 nm−1, corresponding to a wavelength λ ~ 2.85 nm. Inserting this result in the linear dispersion relation we obtain a frequency ν* ~ 10.5 THz. We notice that in the time domain, and at fixed layer index, this frequency corresponds to the inverse of the period of the large amplitude temperature oscillation originating after the transit of the ballistic energy wavefront, as shown in Fig. 3a.

We now compare the predictions of the microscopic model to a macroscopic description of electronic temperature waves based on the phenomenological Dual Phase Lag Model (DPLM)50. The DPLM builds on a modification of Fourier law by the introduction of a delay between the time at which the temperature gradient T is established, t + τT, and the time when the inter-layer heat-flux q sets in, t + τq

$$q(t+{\tau }_{q},z)=-{\kappa }_{T,el}\,{\nabla }_{\perp }T\left(t+{\tau }_{T},z\right).$$

This approach, although phenomenological, allows to easily accounts for both thermal wave damping, and dispersion. Recently, it proved effective in describing phononic temperature wave oscillations in graphite18.

The expansion of Eq. (7) to first order, and its combination with the local conservation of energy at time t, gives rise to a second order parabolic differential equation for the temperature variation ΔT(t, z) = T(t, z) − Tc0. We look for wave-like solutions of this differential equation starting from a temperature pulse triggered at initial time on the top side of the sample slab. Following ref. 18, the pulse can be described by a superposition of plane-waves of real-valued wave vectors k and complex frequencies ν. Underdamped plane-wave solutions for ΔT(t, z) are found if the condition τq > 2τT is met. These temperature waves are characterized by the complex-valued dispersion relation

$$\nu (k,R,\alpha )={\nu }_{1}(k,R,\alpha )+{{{{{{{\rm{i}}}}}}}}{\nu }_{2}(k,R,\alpha ),$$

where ν1,2(k, R, α) depend on the wavevector k and on the parameters \(R=\frac{{\tau }_{T}}{{\tau }_{q}}\) and thermal diffusivity \(\alpha =\frac{{\kappa }_{T,el}}{{C}_{el}}\), κT,el and Cel being the electronic thermal conductivity and specific heat, respectively. The analytical expressions for the real-valued ν1 and ν2, together with the quality factor defined as \(Q(k,R,\alpha )=\frac{{\nu }_{1}}{{\nu }_{2}}\), are reported in the Supplementary Eqs. (1)–(3). It is important to stress that the derivation of the dispersion for the DPLM, Eq. (8), directly follows from energy conservation and the assumption of a delay between the formation of a gradient of temperature and the heat flux. As such, there is no direct connection between the DPLM and the assumptions used for studying the thermal dynamics in the microscopic model. We further notice that, in principle, the quantities R and α do depend on the electronic temperature T. However, since we are looking for relative variation ΔT(t)/Tc0 1 (see Fig. 3a), the temperature dependence may be limited to the initial base temperature Tc0, i.e. R = R(Tc0) and α = α(Tc0).

In order to reveal under which conditions temperature waves are sustained, we exploit the dispersion ν1(k) and its quality factor Q, upon inserting in Eq. (8) the parameters relevant for SrVO3. In particular, we extract the thermodynamic quantities from the literature, whereas we use the sub-picosecond dynamics of the microscopic model to estimate the delay time τq. We first identify the time for setting a variation in the temperature gradient, τT, as the electronic thermalization time. The local thermalization time in SVO is estimated to be as short as ~5 fs on the basis of angle-resolved photoemission spectroscopy51 and optical conductivity36 data (see Supplementary Fig. 1). This timescale is compatible with the attribution in the microscopic model of an instantaneous local temperature on the sub-picosecond timescale, as described in Fig. 3. We thus set τT = 5 fs. On the other hand, the heat-flux dynamics in Fig. 3c shows that the synchronization between T and q starts at ~900 fs, i.e. 500 fs after the ballistic wavefront has transited through the 15th layer. We can thus assume τq 500 fs. Based on these assumptions, we obtain R = τT/τq ~ 0.01, which is well below the threshold R < 0.5 for the observation of a wave-like behaviour18. While the electronic scattering time is expected to weakly depend on the temperature, the temperature dependence of τq is tested by studying the dynamics of the single-band Hubbard model at different base temperatures Tc0. As shown in Supplementary Fig. 2, τq is almost independent of Tc0, thus allowing to assume a temperature independent value of R. The temperature dependence of the wave frequencies is instead retained through the thermal diffusivity α. Specifically, for the case of SVO, Cel = γT with γ = 2.4 × 102 Jm−3 K−252. As for κT,el(T) we retrieve it from the temperature-dependent electrical conductivity, σ(T), of SVO single crystals52 upon application of the Wiedemann-Franz-Lorentz relation: κT,el = LσT, L = 2.44  10−8 WΩK−2 being the Lorentz number. The temperature-dependent κT,el ranges from 10 Wm−1K−1 at 300 K to 20 Wm−1K−1 at 35 K.

With this parameters at hand, in Fig. 6 we show the dispersion relation for the temperature oscillation frequency ν1 (top panel) and the corresponding Q-factor (bottom panel) as a function of wavelength λ = 2π/k and base temperature Tc0. The temperature wave frequency ν ~ 10.5 THz, obtained from the microscopic model at the base temperature Tc0 = 35 K, falls within the range of the allowed frequencies and is compatible with two possible wavelengths, λ ~ 6.5 nm and λ ~ 1.1 nm. These wavelengths correspond to Q-factors ~5 and 0.2, respectively, therefore only the longest wavelength is expected to be detectable. This wavelength falls pretty close to the estimate λ ~ 2.85 nm obtained from the microscopic single-band Hubbard model.

Fig. 6: Temperature wave dispersion in SVO3.
figure 6

Top panel: electronic temperature oscillation frequency vs oscillation’s wavelength, λ, at a base temperatures Tc0 = 35 K (blue line) and 300 K (red line). Inset: Structure of SVO3. Bottom panel: quality factor (colormap) as a function of the oscillation’s wavelength, λ, and the base temperature Tc0. Calculations are based on the dispersion relations, Supplementary Eqs. (1)–(3), upon insertion of input parameters from experiments, α and τT, and from non-equilibrium thermal dynamics results from the layered Hubbard model (see text), τq. In both panels a linear-log plot is adopted.

Given the quite general assumptions on the parameters of the microscopic model and the realistic values used in the phenomenological model, the above comparison shows an overall good agreement between the temperature waves dynamics obtained from the sub-picoseconds dynamics of the single-band Hubbard model and the predictions based on a macroscopic model. Such an agreement further confirms that LCM can sustain, in the hydrodynamic regime, temperature waves with wavelengths and periods fully compatible with state-of-the-art materials growth techniques and time-resolved spectroscopies. More in general, the frequencies and Q-factor values reported in Fig. 6 show that the manifestation of temperature waves in LCM can be observed up to temperatures as high as 300 K. This is the consequence of the fact that the energy scales controlling the electronic dynamics, i.e. t = t = 60 meV and U = 650 meV, correspond to temperatures of 700 K and 7000 K, respectively. At variance with the phononic case, the sub-picosecond electronic hydrodynamic regime is thus expected to be very robust against temperature, giving rise to the emergence of temperature wave-like oscillations in real materials at ambient conditions.

Similarly to the ballistic transport regime, the electronic correlations are key to control the wave-like temperature propagation. In Fig. 7 we report the temperature dynamics at layer n = 15 for different values of the interaction U. We observe that the smaller the interaction the smaller is the temperature oscillation amplitude triggered in the population of cold electrons by the transit of the hot-electron wavefront. The data further show that the temperature oscillation periods, indicated in Fig. 7 by the black arrows, decrease as the interaction U is decreased. In general, as already observed for the case of the ballistic energy propagation, the thermal dynamics of quasiparticles becomes slower as the interaction is increased. This may be traced back to the effect of the correlation-driven renormalization of the quasiparticle effective mass. Therefore, the electronic correlations strength, which controls the quasiparticle effective mass renormalization, can act as a control parameter for the frequency and amplitude of transient temperature waves in LCM.

Fig. 7: Control of temperature wave-like oscillations.
figure 7

Temperature T(t) of the “cold” electronic oscillation at the n = 15 layer for different values of U (the same values used in the inset of Fig. 4) Some of the data have been magnified and the curves shifted for graphical reasons. The horizontal arrows highlight the oscillation periods that match the frequency ν* extracted from spectral analysis of the temperature wave packet.

Recovery of Fourier-like heat transport

After the transit of the ballistic heat wavefront and of the temperature wave packet, the hydrodynamic regime gradually evolves into a more conventional dynamics (t > 0.9 ps in Fig. 3). Here, the non-equilibrium hot-electron population has already left the region of interest, giving rise to a free oscillatory equilibration dynamics of the temperature of “cold" electrons. The wavelength of the temperature oscillation is smaller than that of the temperature wave-packet propagating with speed vg (hydrodynamic regime), as may be seen in Fig. 3b. In the present regime, the oscillation frequency ν gives an intuitive physical picture of the temperature oscillations. Indeed, the oscillation frequency ν corresponds to an energy hν, which exactly matches the energy \(4{t}_{\perp }^{* }\) corresponding to the renormalized bandwidth in the direction perpendicular to the layers (Fig. 1). This fact clearly highlights that the oscillatory dynamics of the temperature is governed by the particle inter-layer hoppings between occupied and empty states due to unbalanced energy distributions on the different layers. In this simple picture, the heat flux left behind by the temperature wave packet freely oscillates with a frequency controlled by \({t}_{\perp }^{* }\), which is the only intrinsic energy scale of the Hubbard Hamiltonian playing a role on the hundreds femtoseconds timescale. The oscillating qn(t) thus acts as the source for the temperature gradient, which instantaneously follows the temperature variation, i.e. without any delay, as expressed in Fourier law \(q\left(t,z\right)=-{\kappa }_{T,el}\,{\nabla }_{\perp }T\left(t,z\right)\). For instance, for U = 0.65 eV one has \({t}_{\perp }^{* }\simeq 0.33{t}_{\perp }=20\,{{{{{{{\rm{meV}}}}}}}}\) and the oscillation periods reads \(h/4{t}_{\perp }^{* }\simeq 50\,{{{{{{{\rm{fs}}}}}}}}\).

Interestingly, we can estimate the electronic thermal conductivity by evaluating the ratio between the oscillation amplitude of the heat flux and that of the temperature gradient, i.e. κT,el = q/T where the symbol represents the oscillation amplitude. As an example of the moderately interacting regime we consider U = 0.45 eV, resulting in \(\frac{U}{{U}_{c}} \sim 0.55\), that is quite far from the Mott transition critical point, U/Uc ~ 1. In doing so we estimate κT,el 440 Wm−1 K−1, which is in the order of typical zero-frequency electronic thermal conductivity of conventional metals. On the other hand, when we increase U, the interactions drive a larger temperature gradient (see Supplementary Fig. 3), which in turn results in a very small value of κT,el53. For instance, when U/Uc = 0.7–0.8, the estimated thermal conductivity is in the range 40–2.5 Wm−1 K−1, a value of the same order of the zero-frequency conductivity reported for SVO52. Thus, despite its simplicity, the thermal dynamics of the layered Hubbard model predicts the correct order of magnitude of Fourier-like thermal conductivity of materials for a very wide range of correlation strengths.

In conclusion, we have proposed layered correlated materials as a platform enabling to access a rich variety of heat transport regimes. We have considered a layered Hubbard model and studied the thermal dynamics triggered by the creation of a non-equilibrium population of hot electrons on the top of the system. The transient dynamics undergoes, on ultrashort space and timescales, a crossover between ballistic energy transport and electronic temperature wave-like oscillations in the hydrodynamic regime. Eventually, the Fourier-like heat transport regime is recovered on the picosecond timescale. Specifically, transition metal oxides thin films and heterostructures, with typical thicknesses and periodicities in the few nanometers range, are here predicted to sustain electronic temperature wave-like oscillations in a parameter space fully compatible with state-of-the-art time-resolved calorimetry techniques54. We stress that, in contrast to the phononic case for which scattering is detrimental to the formation of temperature waves, in the present electronic case the correlation strengths are the key ingredients to observe this phenomenon. As a results, electronic temperature waves can exists for a much broader range of temperatures, and are predicted to persist up to room temperature.

The outreach of our results ranges beyond LCM. Among the most interesting applications we foresee is nanoengineering of superlattices made out of correlated materials, allowing for coherent control of temperature waves in nanodevices. For instance, LCM can be grown in heterostructures with control of the physical properties at the level of single atomic layers55 and with the possibility of engineering artificial periodicities to select high-Q modes of temperature waves. The recent introduction of the temperonic crystal21, i.e. a periodically modulated structure, which behaves like a crystal for temperature waves, provides a new tool to coherently control temperature pulses in correlated heterostructures. Furthermore, strong correlations, and their control via the inter-layer twist angle, have recently been reported in graphene superlattices56,57,58. The present work paves the way to the control of electronic ballistic propagation and to the engineering of nanodevices exploiting the wave-like nature of the electronic heat transfer on the sub-picosecond timescale.


Gutzwiller variational dynamics

In this section we describe the general scheme of the Gutzwiller variational dynamics used in this work to address the thermal dynamics in the layered Hubbard model. We refer the reader to existing literature for an in-depth review about the method41,42,44. We start from the definition of a variational density matrix of the form

$$\hat{\rho }(t)\equiv {{{{{{{\mathcal{P}}}}}}}}(t){\hat{\rho }}_{* }(t){{{{{{{{\mathcal{P}}}}}}}}}^{{{{\dagger}}} }(t),$$

where \({\hat{\rho }}_{* }(t)\) represents the density matrix of an effective system of non-interacting quasiparticles, and \({{{{{{{\mathcal{P}}}}}}}}(t)={\prod }_{in}{{{{{{{{\mathcal{P}}}}}}}}}_{in}(t)\) is a projector, defined by operators \({{{{{{{{\mathcal{P}}}}}}}}}_{in}\) acting on the local Hilbert spaces defined at the site with in-plane index i and layer index n. The projectors \({{{{{{{{\mathcal{P}}}}}}}}}_{in}(t)\) control the relative weights of the local many-body configurations. We obtain dynamical equations of motions by using the Dirac-Frenkel time-variational principle applied to the Gutzwiller ansatz (9). This variational dynamics has been originally introduced for the description of the dynamics of pure quantum states41, whereas extension to the dynamics of mixed states has been presented in ref. 44. By representing the local projectors by a set of matrices \({\hat{{{\Phi }}}}_{in}\) defined on the local Hilbert space, the variational dynamics is defined by the set of coupled equation of motions

$$i\hslash {\partial }_{t}{\hat{\rho }}_{* }(t)=\left[{H}_{* }[\hat{{{\Phi }}}(t)],{\hat{\rho }}_{* }(t)\right]$$
$$i\hslash {\partial }_{t}{\hat{{{\Phi }}}}_{in}(t)=\frac{\delta }{\delta {{{\Phi }}}_{in}^{{{{\dagger}}} }(t)}{{{{{{{\rm{Tr}}}}}}}}\left[{\hat{\rho }}_{* }(t){H}_{* }[\hat{{{\Phi }}}(t)]\right]+{H}_{in}{\hat{{{\Phi }}}}_{in}(t).$$

In the above equations \({H}_{* }[\hat{{{\Phi }}}(t)]\) is a single-particle Hamiltoinan obtained by replacing the fermionic operators in the hopping Hamiltonian H0, \({c}_{in\sigma }^{{{{\dagger}}} }\to {R}_{in}[\hat{{{\Phi }}}]{c}_{in\sigma }^{{{{\dagger}}} }\), with \({R}_{in}[\hat{{{\Phi }}}]\) hopping renormalization parameters determined by the local projectors \({\hat{{{\Phi }}}}_{in}\). Hin is the matrix representation of the Hubbard interaction onto the local Hilbert space. The density matrix \({\hat{\rho }}_{* }(t)\) describes the dynamics of effective quasiparticles with renormalized hopping or, equivalently, enhanced mass. The quasiparticle renormalization is defined as \({Z}_{i}=| {R}_{i}[\hat{{{\Phi }}}]{| }^{2}\ \le \ 1,\) and it is controlled by the local projectors \({\hat{{{\Phi }}}}_{i}\). The local projectors \({\hat{{{\Phi }}}}_{i}\) describe the dynamics of the local degrees of freedom associated with the interaction term. In the single-band case at half-filling, the local degrees of freedom reduce to the excitations of double occupied sites (doublons) at energies ~U. The two dynamics are coupled so that the method is able to capture a non-trivial feedback between the delocalized (quasipartcles) and localized (doublons) degrees of freedom. It is important to mention that Eqs. ((10)–(11)) represent an exact solution of the variational problem only in the limit of a lattice of infinite coordination. For lattices of finite coordination Eqs. ((10)–(11)) represent an approximation to the variation problem, known as the Gutzwiller approximation. In this work, we consider full in-plane translation invariance, so that the matrices Φin depend only on the layer index n. The non-equilibrium dynamics is completely determined by the solution of the set dynamical Eqs. ((10)–(11)) with initial conditions set by the non-equilibration protocol described below.

Non-equilibrium protocol

The non-equilibrium protocol defines the initial condition for the unitary dynamics described in this work. We first solve the variational problem at T = 0, by looking at stationary solutions of Eqs. ((10)–(11)) for a pure quantum state \({\hat{\rho }}_{* }=\left|{{{\Psi }}}_{* }(t)\right\rangle \left\langle {{{\Psi }}}_{* }(t)\right|\). Therefore, we set the temperatures of the two-electronic populations by running an equilibration dynamics in which each layer is coupled to an external bath at temperature Tn. To this extent we supplement the unitary dynamics for the quasiparticles (10) with a dissipative term \({{{{{{{{\mathcal{L}}}}}}}}}_{{{{{{{{\rm{bath}}}}}}}}}\)

$$i\frac{\partial {\rho }_{* }(t)}{\partial t}=\left[{H}_{* }(t),{\rho }_{* }(t)\right]+{{{{{{{{\mathcal{L}}}}}}}}}_{{{{{{\mathrm{bath}}}}}}}[{\rho }_{* }].$$

which is obtained by considering, on each layer, a non-interacting bath defined on the same square lattice of the Hubbard layer. We assume all the baths to be identical, so that the bath Hamiltonian reads \({H}_{{{{{{{{\rm{bath}}}}}}}}}={\sum }_{{{{{{{{\bf{k}}}}}}}}n\sigma }{\epsilon }_{bath}({{{{{{{\bf{k}}}}}}}}){d}_{{{{{{{{\bf{k}}}}}}}}n\sigma }^{{{{\dagger}}} }{d}_{{{{{{{{\bf{k}}}}}}}}n\sigma }^{}\). The interaction with the electrons in the layers is described by local single-particle hopping processes \({H}_{{{{{{{{\rm{bath-sys}}}}}}}}}=v{\sum }_{in\sigma }{c}_{in\sigma }^{{{{\dagger}}} }{d}_{in\sigma }^{}+h.c.=v{\sum }_{{{{{{{{\bf{k}}}}}}}}n\sigma }{c}_{{{{{{{{\bf{k}}}}}}}}n\sigma }^{{{{\dagger}}} }{d}_{{{{{{{{\bf{k}}}}}}}}n\sigma }^{}+h.c.\). By integrating out the baths in the so-called Born-Markov approximation, and by assuming that the baths belonging to different layers are completely decoupled, i.e. \(\langle {d}_{{{{{{{{\bf{k}}}}}}}}n\sigma }^{{{{\dagger}}} }{d}_{{{{{{{{\bf{k}}}}}}}}n^{\prime} \sigma }^{}\rangle \propto {\delta }_{nn^{\prime} }\), we obtain a standard Lindbladt form of the dissipative term containing both single-particle losses, \({{{{{{{{\mathcal{L}}}}}}}}}^{-}\), and gains, \({{{{{{{{\mathcal{L}}}}}}}}}^{+}\),

$${{{{{{{{\mathcal{L}}}}}}}}}_{{{{{{\mathrm{bath}}}}}}}[{\rho }_{* }]={{{{{{{{\mathcal{L}}}}}}}}}^{+}[{\rho }_{* }]+{{{{{{{{\mathcal{L}}}}}}}}}^{-}[{\rho }_{* }],$$


$${{{{{{{{\mathcal{L}}}}}}}}}^{-}[{\rho }_{* }]=\mathop{\sum}\limits_{{{{{{{{\bf{k}}}}}}}}n\sigma }{{{\Gamma }}}_{{{{{{{{\bf{k}}}}}}}}n}^{-}(t)\left[2{c}_{{{{{{{{\bf{k}}}}}}}}n\sigma }^{}{\rho }_{* }(t){c}_{{{{{{{{\bf{k}}}}}}}}n\sigma }^{{{{\dagger}}} }-\{{c}_{{{{{{{{\bf{k}}}}}}}}n\sigma }^{{{{\dagger}}} }{c}_{{{{{{{{\bf{k}}}}}}}}n\sigma }^{},{\rho }_{* }(t)\}\right]$$
$${{{{{{{{\mathcal{L}}}}}}}}}^{+}[{\rho }_{* }]=\mathop{\sum}\limits_{{{{{{{{\bf{k}}}}}}}}n\sigma }{{{\Gamma }}}_{{{{{{{{\bf{k}}}}}}}}n}^{+}(t)\left[2{c}_{{{{{{{{\bf{k}}}}}}}}n\sigma }^{{{{\dagger}}} }{\rho }_{* }(t){c}_{{{{{{{{\bf{k}}}}}}}}n\sigma }^{}-\{{c}_{{{{{{{{\bf{k}}}}}}}}n\sigma }^{}{c}_{{{{{{{{\bf{k}}}}}}}}n\sigma }^{{{{\dagger}}} },{\rho }_{* }(t)\}\right].$$

The couplings \({{{\Gamma }}}_{{{{{{{{\bf{k}}}}}}}}}^{\pm }(t)\) are defined by the distribution functions of the baths as

$${{{\Gamma }}}_{{{{{{{{\bf{k}}}}}}}}n}^{+}\equiv r(t){{\Gamma }}{N}_{n}^{{{{{{\mathrm{bath}}}}}}}({{{{{{{\bf{k}}}}}}}})\,\qquad {{{\Gamma }}}_{{{{{{{{\bf{k}}}}}}}}n}^{-}\equiv r(t){{\Gamma }}(1-{N}_{n}^{{{{{{\mathrm{bath}}}}}}}({{{{{{{\bf{k}}}}}}}})),$$

where Γ is a constant proportionals to the couplings between layer and baths Γ v2 and r(t) = θ(t − t)θ(− t) is a box function function highlighting that the coupling is switched-on at a time t < 0 and switched-off at time t = 0 when the equilibration is reached. \({N}_{n}^{{{{{{\mathrm{bath}}}}}}}({{{{{{{\bf{k}}}}}}}})=f({\epsilon }_{{{{{{\mathrm{bath}}}}}}}({{{{{{{\bf{k}}}}}}}}),{T}_{{{{{{\mathrm{bath}}}}}}}(n))\) is the Fermi-distribution of the nth bath at temperature Tbath(n). The structure of the couplings reflects the fact that the probability of jumping from the bath to the layer (gain) in a given momentum state k is proportional to the occupation number in the corresponding state in the bath. On the contrary, the probability to jump from the layer to the bath (loss), is proportional to the probability of finding an empty state in the bath. The Eq. (12) is solved together with the coupled equation for the local degrees of freedom (11), so that the coupling of quasiparticles with the baths has an indirect effect also onto the local degrees of freedom and, in turn, on the quasiparticle renormalization Z. The quasiparticle renormalization corresponding to the mass enhancement parameter mentioned in the main text m/m*  0.3 corresponds to the value obtained at the end of the equilibration step. The absolute value t corresponds to the time needed to reach equilibration, being larger the smaller the couplings Γ. In practice, it is found that the stationary solution does not depend on the strength of the coupling, so that the equilibration can be reached considering an equilibration dynamics much shorter than the picosecond dynamics discussed in the main text.

When equilibration is reached, the density matrix \(\hat{\rho }(t)\) approaches a stationary value \({\hat{\rho }}_{{{{{{{{\rm{equ}}}}}}}}}\equiv {{{{{{{{\mathcal{P}}}}}}}}}_{{{{{{{{\rm{equ}}}}}}}}}{\hat{\rho }}_{* }^{{{{{{{{\rm{equ}}}}}}}}}{{{{{{{{\mathcal{P}}}}}}}}}_{{{{{{{{\rm{equ}}}}}}}}}^{{{{\dagger}}} }\), independent of time. In the equilibrated state, the quasiparticle occupation numbers \({N}_{n}^{{{{{{{{\rm{equ}}}}}}}}}({{{{{{{\bf{k}}}}}}}})\equiv {{{{{{{\rm{Tr}}}}}}}}[{\hat{\rho }}_{* }^{{{{{{{{\rm{equ}}}}}}}}}{c}_{{{{{{{{\bf{k}}}}}}}}n\sigma }^{{{{\dagger}}} }{c}_{{{{{{{{\bf{k}}}}}}}}n\sigma }^{}]\) become equal to the occupation numbers of the baths. By changing variable from the quasi-momentum k to energy using the bath dispersions ϵbath(k), the equilibration condition, in terms of the quasiparticle energy distribution function, reads

$${N}_{n}^{{{{{{{{\rm{equ}}}}}}}}}(\epsilon )\equiv {N}_{n}^{{{{{{{{\rm{equ}}}}}}}}}({\epsilon }_{{{{{{\mathrm{bath}}}}}}}({{{{{{{\bf{k}}}}}}}}))=f(\epsilon ,{T}_{n}).$$

In Fig. 1 we show two energy distribution functions for a surface (n ≤ 5), and a bulk (n > 5) layer, respectively, obtained after the equilibration step.

Non-equilibrium distribution functions

Here, we describe the definition of the non-equilibrium distribution functions \({N}_{n}^{{{{{{{{\rm{neq}}}}}}}}}(\epsilon ,t)\). We consider the time evolution of the layer-dependent quasiparticle occupation numbers

$${N}_{n}({{{{{{{\bf{k}}}}}}}},t)\equiv {{{{{{{\rm{Tr}}}}}}}}\left[{\hat{\rho }}_{* }(t){c}_{{{{{{{{\bf{k}}}}}}}}n\sigma }^{{{{\dagger}}} }{c}_{{{{{{{{\bf{k}}}}}}}}n\sigma }^{}\right].$$

In order to track the evolution of the occupation numbers with respect to the initial condition (17), we compare Nn(k, t) with the occupation numbers that would be obtained after an equilibration step with a bath at fixed temperature. To this extent, we define non-equilibrium distribution functions from the layer-dependent occupation numbers, Eq. (18), in the same way as done in Eq. (17) for the equilibrated distribution functions, namely

$${N}_{n}^{{{{{{{{\rm{neq}}}}}}}}}(\epsilon ,t)\equiv {N}_{n}^{{{{{{{{\rm{neq}}}}}}}}}({\epsilon }_{{{{{{\mathrm{bath}}}}}}}({{{{{{{\bf{k}}}}}}}}),t).$$

We mention that such a definition of a non-equilibrium distribution is needed as the Gutzwiller method does not give direct access to spectral quantities, but only to quantities integrated in frequency. However, the definition (19) is physically motivated by the need of having a meaningful comparison with the distribution functions after the equilibration step, Eq. (17) which are known, being entirely determined by the baths. In particular, at t = 0, the definition Eq. (19) correctly reproduces \({N}_{n}^{{{{{{{{\rm{neq}}}}}}}}}(\epsilon ,t=0)={N}_{n}^{{{{{{{{\rm{equ}}}}}}}}}(\epsilon )\), with \({N}_{n}^{{{{{{{{\rm{equ}}}}}}}}}(\epsilon )\) defined in Eq. (17), whereas at t > 0, when all the couplings with the baths are set to zero, the relation Eq. (17) is no more satisfied and \({N}_{n}^{{{{{{{{\rm{neq}}}}}}}}}(\epsilon ,t)\) can be expressed as a superposition of two equilibrium Fermi distributions, Eq. (4). This allows for the description of the thermal dynamics in terms of a propagating hot electronic population and the temperature dynamics of the cold electronic population.

Data availability

The data generated in this study, including the code used to generate the data, have been deposited in the YARETA database59.


  1. 1.

    Chen, G. Nanoscale Energy Transport and Conversion (Oxford University Press, 2005).

  2. 2.

    Volz, S. et al. Nanophononics: state of the art and perspectives. Eur. Phys. J. B 89, 15 (2016).

    ADS  Google Scholar 

  3. 3.

    Li, N. et al. Colloquium: Phononics: manipulating heat flow with electronic analogs and beyond. Rev. Mod. Phys. 84, 1045–1066 (2012).

    ADS  Google Scholar 

  4. 4.

    Luo, T. & Chen, G. Nanoscale heat transfer—from computation to experiment. Phys. Chem. Chem. Phys. 15, 3389–3412 (2013).

    CAS  PubMed  Google Scholar 

  5. 5.

    Cahill, D. G. et al. Nanoscale thermal transport. ii. 2003–2012. Appl. Phys. Rev. 1, 011305 (2014).

    ADS  Google Scholar 

  6. 6.

    Siemens, M. E. et al. Quasi-ballistic thermal transport from nanoscale interfaces observed using ultrafast coherent soft x-ray beams. Nat. Mater. 9, 26 (2010).

    ADS  CAS  PubMed  Google Scholar 

  7. 7.

    Minnich, J. et al. Thermal conductivity spectroscopy technique to measure phonon mean free paths. Phys. Rev. Lett. 107, 095901 (2011).

    ADS  CAS  PubMed  Google Scholar 

  8. 8.

    Johnson, J. A. et al. Direct measurement of room-temperature nondiffusive thermal transport over micron distances in a silicon membrane. Phys. Rev. Lett. 110, 025901 (2013).

    ADS  PubMed  Google Scholar 

  9. 9.

    Hoogeboom-Pot, K. M. et al. A new regime of nanoscale thermal transport: collective diffusion increases dissipation efficiency. Proc. Natl Acad. Sci. USA 112, 4846 (2015).

    ADS  CAS  PubMed  PubMed Central  Google Scholar 

  10. 10.

    Chen, X., Hua, C., Zhang, H., Ravichandran, N. K. & Minnich, A. J. Quasiballistic thermal transport from nanoscale heaters and the role of the spatial frequency. Phys. Rev. Appl. 10, 054068 (2018).

    ADS  CAS  Google Scholar 

  11. 11.

    Frazer, T. D. et al. Engineering nanoscale thermal transport: Size- and spacing-dependent cooling of nanostructures. Phys. Rev. Appl. 11, 024042 (2019).

    ADS  CAS  Google Scholar 

  12. 12.

    Guyer, R. A. & Krumhansl, J. A. Thermal conductivity, second sound, and phonon hydrodynamic phenomena in nonmetallic crystals. Phys. Rev. 148, 778–788 (1966).

    ADS  CAS  Google Scholar 

  13. 13.

    Beck, H., Meier, P. & Thellun, A. Phonon hydrodynamics in solids. Phys. Stat. Sol. a 24, 11–63 (1974).

    ADS  CAS  Google Scholar 

  14. 14.

    Lee, S., Broido, D., Esfarjani, K. & Chen, G. Hydrodynamic phonon transport in suspended graphene. Nat. Commun. 6, 6290 (2015).

    ADS  CAS  PubMed  Google Scholar 

  15. 15.

    Ding, Z. et al. Phonon hydrodynamic heat conduction and knudsen minimum in graphite. Nano Lett. 18, 638–649 (2018).

    ADS  CAS  PubMed  Google Scholar 

  16. 16.

    Cepellotti, A. et al. Phonon hydrodynamics in two-dimensional materials. Nat. Commun. 6, 1–7 (2015).

    Google Scholar 

  17. 17.

    Li, X. & Lee, S. Crossover of ballistic, hydrodynamic, and diffusive phonon transport in suspended graphene. Phys. Rev. B 99, 085202 (2019).

    ADS  CAS  Google Scholar 

  18. 18.

    Gandolfi, M., Benetti, G., Glorieux, C., Giannetti, C. & Banfi, F. Accessing temperature waves: a dispersion relation perspective. Int. J. Heat. Mass Transf. 143, 118553 (2019).

    Google Scholar 

  19. 19.

    Huberman, S. et al. Observation of second sound in graphite at temperatures above 100 K. Science 364, 375 (2019).

    ADS  CAS  PubMed  Google Scholar 

  20. 20.

    Zhang, Y. et al. Coherent modulation of the electron temperature and electron–phonon couplings in a 2D material. Proc. Natl Acad. Sci. USA 117, 8788–8793 (2020).

    CAS  PubMed  PubMed Central  Google Scholar 

  21. 21.

    Gandolfi, M., Giannetti, C. & Banfi, F. Temperonic crystal: a superlattice for temperature waves in graphene. Phys. Rev. Lett. 125, 265901 (2020).

    ADS  CAS  PubMed  Google Scholar 

  22. 22.

    Cepellotti, A. & Marzari, N. Transport waves as crystal excitations. Phys. Rev. Mater. 1, 045406 (2017).

    Google Scholar 

  23. 23.

    Torres, P., Alvarez, F. X., Cartoixà, X. & Rurali, R. Thermal conductivity and phonon hydrodynamics in transition metal dichalcogenides from first-principles. 2D Mater. 6, 035002 (2019).

    CAS  Google Scholar 

  24. 24.

    Machida, Y., Matsumoto, N., Isono, T. & Behnia, K. Phonon hydrodynamics and ultrahigh–room-temperature thermal conductivity in thin graphite. Science 367, 309–312 (2020).

    ADS  CAS  PubMed  Google Scholar 

  25. 25.

    Beardo, A. et al. Observation of second sound in a rapidly varying temperature field in Ge. Sci. Adv. 7, eabg4677 (2021).

    ADS  CAS  PubMed  PubMed Central  Google Scholar 

  26. 26.

    Melis, C., Fugallo, G. & Colombo, L. Room temperature second sound in cumulene. Phys. Chem. Chem. Phys.23, 15275–15281 (2021).

    CAS  PubMed  Google Scholar 

  27. 27.

    Simoncelli, M., Marzari, N. & Cepellotti, A. Generalization of Fourier’s law into viscous heat equations. Phys. Rev. X 10, 011019 (2020).

    CAS  Google Scholar 

  28. 28.

    Gandolfi, M. et al. Emergent ultrafast phenomena in correlated oxides and heterostructures. Phys. Scr. 92, 034004 (2017).

    ADS  Google Scholar 

  29. 29.

    Sobolev, S. Local nonequilibrium electron transport in metals after femtosecond laser pulses: a multi-temperature hyperbolic model. Nanoscale Microscale Thermophys. Eng. 1–13 (2021).

  30. 30.

    Hartnoll, S. A. Theory of universal incoherent metallic transport. Nat. Phys. 11, 54–61 (2015).

    CAS  Google Scholar 

  31. 31.

    Lee, S. et al. Anomalously low electronic thermal conductivity in metallic vanadium dioxide. Science 355, 371–374 (2017).

    ADS  CAS  PubMed  Google Scholar 

  32. 32.

    Tokura, Y., Kawasaki, M. & Nagaosa, N. Emergent functions of quantum materials. Nat. Phys. 13, 1056–1068 (2017).

    CAS  Google Scholar 

  33. 33.

    Basov, D. N., Averitt, N. D. & Hsieh, D. Towards properties on demand in quantum materials. Nat. Mat. 16, 1077–1088 (2017).

    CAS  Google Scholar 

  34. 34.

    Moyer, J. A., Eaton, C. & Engel-Herbert, R. Highly conductive srvo3 as a bottom electrode for functional perovskite oxides. Adv. Mater. 25, 3578–3582 (2013).

    CAS  PubMed  Google Scholar 

  35. 35.

    Zhong, Z. et al. Electronics with correlated oxides: SrVO3/SrTiO3 as a mott transistor. Phys. Rev. Lett. 114, 246401 (2015).

    ADS  PubMed  Google Scholar 

  36. 36.

    Zhang, L. et al. Correlated metals as transparent conductors. Nat. Mater. 15, 204–210 (2015).

    ADS  PubMed  Google Scholar 

  37. 37.

    Ordonez-Miranda, J., Ezzahri, Y., Joulain, K., Drevillon, J. & Alvarado-Gil, J. J. Modeling of the electrical conductivity, thermal conductivity, and specific heat capacity of VO2. Phys. Rev. B 98, 075144 (2018).

    ADS  CAS  Google Scholar 

  38. 38.

    Cesarini, G. et al. Quantitative evaluation of emission properties and thermal hysteresis in the mid-infrared for a single thin film of vanadium dioxide on a silicon substrate. Int. J. Therm. Sci. 146, 106061 (2019).

    CAS  Google Scholar 

  39. 39.

    Ben-Abdallah, P. & Biehs, S.-A. Near-field thermal transistor. Phys. Rev. Lett. 112, 044301 (2014).

    ADS  PubMed  Google Scholar 

  40. 40.

    Ordonez-Miranda, J., Ezzahri, Y., Tiburcio-Moreno, J. A., Joulain, K. & Drevillon, J. Radiative thermal memristor. Phys. Rev. Lett. 123, 025901 (2019).

    ADS  CAS  PubMed  Google Scholar 

  41. 41.

    Fabrizio M. The Out-of-Equilibrium Time-Dependent Gutzwiller Approximation. In New Materials for Thermoelectric Applications: Theory and Experiment. NATO Science for Peace and Security Series B: Physics and Biophysics. (eds Zlatic, V. & Hewson, A.) (Springer, Dordrecht, 2013)

  42. 42.

    Mazza, G., Amaricci, A., Capone, M. & Fabrizio, M. Electronic transport and dynamics in correlated heterostructures. Phys. Rev. B 91, 195124 (2015).

    ADS  Google Scholar 

  43. 43.

    Sandri, M., Capone, M. & Fabrizio, M. Finite-temperature gutzwiller approximation and the phase diagram of a toy model for v2o3. Phys. Rev. B 87, 205108 (2013).

    ADS  Google Scholar 

  44. 44.

    Lanatà, N., Deng, X. & Kotliar, G. Finite-temperature gutzwiller approximation from the time-dependent variational principle. Phys. Rev. B 92, 081108 (2015).

    ADS  Google Scholar 

  45. 45.

    Georges, A., Kotliar, G., Krauth, W. & Rozenberg, M. J. Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions. Rev. Mod. Phys. 68, 13–125 (1996).

    ADS  MathSciNet  CAS  Google Scholar 

  46. 46.

    Mingo, N. & Broido, D. A. Carbon nanotube ballistic thermal conductance and its limits. Phys. Rev. Lett. 95, 096105 (2005).

    ADS  CAS  PubMed  Google Scholar 

  47. 47.

    Muñoz, E., Lu, J. & Yakobson, B. I. Ballistic thermal conductance of graphene ribbons. Nano Lett. 10, 1652–1656 (2010).

    ADS  PubMed  Google Scholar 

  48. 48.

    Bae, M.-H. et al. Ballistic to diffusive crossover of heat flow in graphene ribbons. Nat. Commun. 4, 1734 (2013).

    ADS  PubMed  Google Scholar 

  49. 49.

    Caddeo, C. et al. Thermal boundary resistance from transient nanocalorimetry: a multiscale modeling approach. Phys. Rev. B 95, 085306 (2017).

    ADS  Google Scholar 

  50. 50.

    Tzou, D. Y. Macro-to Microscale Heat Transfer: The Lagging Behavior (John Wiley & Sons, 2014).

  51. 51.

    Aizaki, S. et al. Self-energy on the low- to high-energy electronic structure of correlated metal SrVO3. Phys. Rev. Lett. 109, 056401 (2012).

    ADS  CAS  PubMed  Google Scholar 

  52. 52.

    Inoue, I. H., Goto, O., Makino, H., Hussey, N. E. & Ishikawa, M. Bandwidth control in a perovskite-type 3d1-correlated metal Ca1−xSrxVO3. I. Evolution of the electronic properties and effective mass. Phys. Rev. B 58, 4372–4383 (1998).

    ADS  CAS  Google Scholar 

  53. 53.

    Karrasch, C., Kennes, D. M. & Heidrich-Meisner, F. Thermal conductivity of the one-dimensional Fermi-Hubbard model. Phys. Rev. Lett. 117, 116401 (2016).

    ADS  CAS  PubMed  Google Scholar 

  54. 54.

    Giannetti, C. et al. Ultrafast optical spectroscopy of strongly correlated materials and high-temperature superconductors: a non-equilibrium approach. Adv. Phys. 65, 58 (2016).

    ADS  CAS  Google Scholar 

  55. 55.

    Hwang, H. Y. et al. Emergent phenomena at oxide interfaces. Nat. Mater. 11, 103–113 (2012).

    ADS  MathSciNet  CAS  PubMed  Google Scholar 

  56. 56.

    Cao, Y. et al. Correlated insulator behaviour at half-filling in magic-angle graphene superlattices. Nature 556, 80?84 (2018).

    PubMed  Google Scholar 

  57. 57.

    Seyler, K. L. et al. Signatures of moiré-trapped valley excitons in MoSe2/WSe2 heterobilayers. Nature 567, 66–70 (2019).

    ADS  CAS  PubMed  Google Scholar 

  58. 58.

    Lu, X. et al. Superconductors, orbital magnets and correlated states in magic-angle bilayer graphene. Nature 574, 653–657 (2019).

    ADS  CAS  PubMed  Google Scholar 

  59. 59.

    Mazza, G., Gandolfi, M., Capone, M., Banfi, F. & Giannetti, C. Dataset “Ultrafast thermal dynamics in strongly correlated systems". YARETA (2021)

Download references


G.M. acknowledges financial support from the Swiss National Science Foundation through an AMBIZIONE grant. Part of this work has been supported from the European Research Council (ERC-319286-QMAC). M.G. acknowledges financial support from the CNR Joint Laboratories program 2019–2021 and Project No. SAC.AD002.026 (OMEN). F.B. acknowledges financial support from Université de Lyon in the frame of the IDEXLYON Project—Programme Investissements d’ Avenir (ANR-16-IDEX-0005) and from Université Claude Bernard Lyon 1 through the BQR Accueil EC 2019 grant and from CNRS (Delegation CNRS 2021–2022). C.G. acknowledges support from Università Cattolica del Sacro Cuore through D.2.2 and D.3.1 grants. M.C. and C.G. acknowledge financial support from MIUR through the PRIN 2015 (Prot 2015C5SEJJ001) and PRIN 2017 (Prot. 20172H2SC4_005) programs.

Author information




G.M. developed the theoretical framework and performed the numerical study of the non-equilibrium dynamics in the layered Hubbard model. M.G. performed the temperature wave analysis within the Dual Phase Lag Model. G.M., M.G., M.C., F.B. and C.G. contributed to conceive the project, to discuss the results and to write the manuscript.

Corresponding authors

Correspondence to Giacomo Mazza, Francesco Banfi or Claudio Giannetti.

Ethics declarations

Competing interests

The authors declare no competing interests.

Additional information

Peer review information Nature Communications thanks Jong Han, Sebastian Volz the other, anonymous, reviewer for their contribution to the peer review of this work.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Mazza, G., Gandolfi, M., Capone, M. et al. Thermal dynamics and electronic temperature waves in layered correlated materials. Nat Commun 12, 6904 (2021).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:


By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.


Quick links

Nature Briefing

Sign up for the Nature Briefing newsletter — what matters in science, free to your inbox daily.

Get the most important science stories of the day, free in your inbox. Sign up for Nature Briefing