Abstract
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 wavelike 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, wavelike 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.
Introduction
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 achieved^{1,2,3,4,5}. For instance, nonFourier heat transport regimes have been reported for hot spots dimensions inferior to the phonon mean freepath^{6,7,8}, in which energy is ballistically carried point to point, or have been engineered via nanopatterning of dielectric substrates^{9,10,11}. As a consequence of the existence of two nonthermal populations, wavelike thermal transport, often referred to as second sound^{12,13}, has been predicted in graphene, both in the frame of microscopic^{14,15,16,17} and macroscopic models^{18}. Temperature wavelike phenomena have been recently observed at high temperatures in graphene^{19} and 2D materials^{20} on subnanosecond timescales and scheme for their coherent control have been proposed^{21}. So far most of the effort has been devoted to phononic nonFourier heat transport^{17,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 proposed^{27}. On the contrary, electronic nonFourier heat transport remains relatively unexplored^{18,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 manybody properties, such as collective and decoupled diffusion of energy and charge^{30,31}. Tuning the interaction strength thus opens the possibility to investigate transport regimes with no counterpart in conventional weaklyinteracting materials^{32,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 intralayer 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 interlayer and intralayer hopping terms, may thus act as a tuning parameter to control the relative inter and intralayer 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.
We investigate the possibility for unconventional heat transport regimes by focusing on the impulsive thermal dynamics of the layered singleband 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 subpicosecond 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 SrVO_{3} (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 materials^{34}, to Mott transistors^{35} and transparent conductors^{36}. We argue that the degree of correlation of SVO, as measured by the interaction strength, is such that ballistic transport first, and wavelike thermal transport afterwords, are accessible on the subpicosecond 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 nonFourier heat transport.
The present work rationalizes the microscopic interactions underlying unconventional electronic heat transfer phenomena in LCM. Our findings enlarge the functionalities of quantum materials^{32,33} to the realm of nanoscale heat transport^{37}, beyond the case of radiative energy transfer^{38,39,40}.
Results and discussion
Model and observables
The Hamiltonian of the layered Hubbard model reads:
with
and
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 interplane hopping amplitudes. The sum of the inplane 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 inplane translational invariance so that we can introduce an inplane momentum k = (k_{x}, k_{y}) 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 (halffilling) corresponding to the perfect particlehole 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 interactiondriven 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 SrVO_{3}.
We study the nonequilibrium thermal dynamics in the frame of model (1) by means of a timedependent variational approach based on the generalized Gutzwiller approximation for layered systems^{41,42} (see Methods). This approach provides a versatile tool for describing in a nonperturbative 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 structure^{43,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 noninteracting 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 manybody configurations. The mutual feedback between the dynamics of the localized degrees of freedom and the lowenergy 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 noninteracting limit Z(U = 0) = 1, whereas at finite interaction Z < 1 and decreases as a function of U. Eventually, for a critical interaction strength, U_{c}, the system undergoes a metaltoinsulator Mott transition, corresponding to a vanishing quasiparticle weight, i.e. Z(U_{c}) → 0. In this regime, quasiparticle excitations are completely suppressed and the dynamics becomes dominated by highenergy incoherent excitations at energies ~ U^{45}. 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 pulse^{28}. We mimic the impulsive excitation by considering the twostep nonequilibrium protocol (see Methods): first, we initialize the system in a nonequilibrium state characterized by twoelectronic 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 noninteracting electrons at temperatures T_{bath}(n). All the layers are at a base temperature T_{bath}(n) = T_{c0} except the five topmost which are set at T_{bath}(n) = T_{h} ≫ T_{c0}, for n = 1, …, 5, a temperature characterizing the hot electrons. In the notation T_{c0}, the subscript “c” stands for cold and “0” indicates the instant right after the impulsive excitation, whereas the subscript “h” in T_{h} stands for hot. Throughout the paper, the temperatures used in the equilibration step are T_{c0} ≃ 35 K for the cold temperature in the bulk and T_{h} = 10 × T_{c0} 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 nonequilibrium 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 layerdependent quasiparticle nonequilibrium 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 T_{h}, 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 layerdependent temperature T(n, t) and of weight ρ_{cold}(n, t) ≡ 1 − ρ_{hot}(n, t)
with 0 ≤ ρ_{hot}(n, t) ≤ 1, f_{hot} = f(ϵ, T_{h}) and f_{cold}(n, t) = f(ϵ, T(n, t)). For any instant of time t and layer index n, we fit N^{neq}(ϵ, 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 T_{h}, 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 nonequilibrium distribution functions are described by different populations of hot electrons, and different cold temperatures T(n, t).
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 T_{c0}. 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)f_{hot}) 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 twotemperature 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 timedependent 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 spatiotemporal 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) = T_{n}, with T_{n} 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
where \(\hat{\rho }(t)\) represents the timeevolved 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 E_{n} 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 q_{n} at layer n, and along the zdirection perpendicular to the planes, by applying the continuity equation
where the discrete spatial derivative is defined with respect to the interlayer 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 nonconventional heat transport on the subpicosecond 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/T_{c0}(n, t) with ΔT = T(n, t) − T_{c0}, 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 = T_{c0}. At t ~ 150 fs the perturbation reaches the n = 15 layer and the dynamics that follows can be neatly divided in three steps.

(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 hotelectron 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 quasiequilibrium distribution.

(ii)
For t ≳ 400 fs the hotelectron 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 T_{c0} indicating that the transit of the ballistic front of hot electrons induced the heating of the population of cold electrons on the layer.

(iii)
Eventually the system equilibrates for t ≳ 0.9 ps with the residual damped temperature oscillations converging to T_{c0}.
We gain further insight into the thermal dynamics by comparing the dynamics of the local cold electronic temperature with the heat flux q_{n}(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 heatflux 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) = (T_{n+1}(t) − T_{n}(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. q_{n}(t) ∝−∇_{⊥}T(t), indicating that the heat transfer process is well described by a Fourierlike heat transfer law.
At intermediate times (0.4 < t < 0.9 ps), before Fourierlike 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 q_{n}(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 Fourierlike 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 subpopulations 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 heatflux ballistic front. As it will be further discussed, this feature can be considered as a temperature wavepacket propagating through the system.
Summarizing, the subpicosecond 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 wavelike propagation of the electronic temperature; (iii) a Fourierlike 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 T_{h}. 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 interlayer 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 interlayer motion and a smaller effective hopping, \({t}_{\perp }^{* }=Z{t}_{\perp }\).
We show these effects in Fig. 4 where we report the spatiotemporal dynamics of the hotelectron 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 v_{b} as a function of U for t_{⊥}/t_{∥} = 1. v_{b} is defined as the slope of the white dashed line in Fig. 4. The correlationinduced renormalization of \({t}_{\perp }^{* }\) strongly suppresses the energy propagation along the zdirection. For the sake of applications, we note that, in nanosystems with sizes of the order of the ballistic mean freepath, the thermal conductivity becomes a sizedependent property^{46,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.
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 counterpart^{13,14,19}. Similarly to the phononic case, this hydrodynamic regime manifests itself by a wavelike 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 wavelike 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 wavepacket group velocity from the simple relation \({X}_{\min }={v}_{g}t\), with v_{g} ~ 30 nm/ps of the same order of magnitude of the ballistic energy wavefront velocity. A similar result is obtained when tracking the timedependent 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} = v_{g}k/2π.
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 interlayer heatflux q sets in, t + τ_{q}
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 graphite^{18}.
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) − T_{c0}. We look for wavelike 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 planewaves of realvalued wave vectors k and complex frequencies ν. Underdamped planewave solutions for ΔT(t, z) are found if the condition τ_{q} > 2τ_{T} is met. These temperature waves are characterized by the complexvalued dispersion relation
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 C_{el} being the electronic thermal conductivity and specific heat, respectively. The analytical expressions for the realvalued ν_{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)∣/T_{c0} ≪ 1 (see Fig. 3a), the temperature dependence may be limited to the initial base temperature T_{c0}, i.e. R = R(T_{c0}) and α = α(T_{c0}).
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 SrVO_{3}. In particular, we extract the thermodynamic quantities from the literature, whereas we use the subpicosecond 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 angleresolved photoemission spectroscopy^{51} and optical conductivity^{36} data (see Supplementary Fig. 1). This timescale is compatible with the attribution in the microscopic model of an instantaneous local temperature on the subpicosecond timescale, as described in Fig. 3. We thus set τ_{T} = 5 fs. On the other hand, the heatflux 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 wavelike behaviour^{18}. 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 singleband Hubbard model at different base temperatures T_{c0}. As shown in Supplementary Fig. 2, τ_{q} is almost independent of T_{c0}, 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, C_{el} = γT with γ = 2.4 × 10^{2} Jm^{−3} K^{−2} ^{52}. As for κ_{T,el}(T) we retrieve it from the temperaturedependent electrical conductivity, σ(T), of SVO single crystals^{52} upon application of the WiedemannFranzLorentz relation: κ_{T,el} = LσT, L = 2.44 ⋅ 10^{−8} WΩK^{−2} being the Lorentz number. The temperaturedependent κ_{T,el} ranges from ≃10 Wm^{−1}K^{−1} at 300 K to ≃20 Wm^{−1}K^{−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 Qfactor (bottom panel) as a function of wavelength λ = 2π/k and base temperature T_{c0}. The temperature wave frequency ν_{∗} ~ 10.5 THz, obtained from the microscopic model at the base temperature T_{c0} = 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 Qfactors ~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 singleband Hubbard model.
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 subpicoseconds dynamics of the singleband 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 stateoftheart materials growth techniques and timeresolved spectroscopies. More in general, the frequencies and Qfactor 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 subpicosecond electronic hydrodynamic regime is thus expected to be very robust against temperature, giving rise to the emergence of temperature wavelike oscillations in real materials at ambient conditions.
Similarly to the ballistic transport regime, the electronic correlations are key to control the wavelike 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 hotelectron 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 correlationdriven 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.
Recovery of Fourierlike 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 nonequilibrium hotelectron 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 wavepacket propagating with speed v_{g} (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 interlayer 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 q_{n}(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/U_{c} ~ 1. In doing so we estimate κ_{T,el} ≃ 440 Wm^{−1} K^{−1}, which is in the order of typical zerofrequency 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,el}^{53}. For instance, when U/U_{c} = 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 zerofrequency conductivity reported for SVO^{52}. Thus, despite its simplicity, the thermal dynamics of the layered Hubbard model predicts the correct order of magnitude of Fourierlike 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 nonequilibrium 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 wavelike oscillations in the hydrodynamic regime. Eventually, the Fourierlike 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 wavelike oscillations in a parameter space fully compatible with stateoftheart timeresolved calorimetry techniques^{54}. 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 layers^{55} and with the possibility of engineering artificial periodicities to select highQ modes of temperature waves. The recent introduction of the temperonic crystal^{21}, 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 interlayer twist angle, have recently been reported in graphene superlattices^{56,57,58}. The present work paves the way to the control of electronic ballistic propagation and to the engineering of nanodevices exploiting the wavelike nature of the electronic heat transfer on the subpicosecond timescale.
Methods
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 indepth review about the method^{41,42,44}. We start from the definition of a variational density matrix of the form
where \({\hat{\rho }}_{* }(t)\) represents the density matrix of an effective system of noninteracting 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 inplane index i and layer index n. The projectors \({{{{{{{{\mathcal{P}}}}}}}}}_{in}(t)\) control the relative weights of the local manybody configurations. We obtain dynamical equations of motions by using the DiracFrenkel timevariational principle applied to the Gutzwiller ansatz (9). This variational dynamics has been originally introduced for the description of the dynamics of pure quantum states^{41}, 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
In the above equations \({H}_{* }[\hat{{{\Phi }}}(t)]\) is a singleparticle Hamiltoinan obtained by replacing the fermionic operators in the hopping Hamiltonian H_{0}, \({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}\). H_{in} 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 singleband case at halffilling, 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 nontrivial 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 inplane translation invariance, so that the matrices Φ_{in} depend only on the layer index n. The nonequilibrium dynamics is completely determined by the solution of the set dynamical Eqs. ((10)–(11)) with initial conditions set by the nonequilibration protocol described below.
Nonequilibrium protocol
The nonequilibrium 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 twoelectronic populations by running an equilibration dynamics in which each layer is coupled to an external bath at temperature T_{n}. To this extent we supplement the unitary dynamics for the quasiparticles (10) with a dissipative term \({{{{{{{{\mathcal{L}}}}}}}}}_{{{{{{{{\rm{bath}}}}}}}}}\)
which is obtained by considering, on each layer, a noninteracting 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 singleparticle hopping processes \({H}_{{{{{{{{\rm{bathsys}}}}}}}}}=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 socalled BornMarkov 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 singleparticle losses, \({{{{{{{{\mathcal{L}}}}}}}}}^{}\), and gains, \({{{{{{{{\mathcal{L}}}}}}}}}^{+}\),
with
The couplings \({{{\Gamma }}}_{{{{{{{{\bf{k}}}}}}}}}^{\pm }(t)\) are defined by the distribution functions of the baths as
where Γ is a constant proportionals to the couplings between layer and baths Γ ∝ v^{2} and r(t) = θ(t − t_{−∞})θ(− t) is a box function function highlighting that the coupling is switchedon at a time t_{−∞} < 0 and switchedoff 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 Fermidistribution of the nth bath at temperature T_{bath}(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 quasimomentum k to energy using the bath dispersions ϵ_{bath}(k), the equilibration condition, in terms of the quasiparticle energy distribution function, reads
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.
Nonequilibrium distribution functions
Here, we describe the definition of the nonequilibrium distribution functions \({N}_{n}^{{{{{{{{\rm{neq}}}}}}}}}(\epsilon ,t)\). We consider the time evolution of the layerdependent quasiparticle occupation numbers
In order to track the evolution of the occupation numbers with respect to the initial condition (17), we compare N_{n}(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 nonequilibrium distribution functions from the layerdependent occupation numbers, Eq. (18), in the same way as done in Eq. (17) for the equilibrated distribution functions, namely
We mention that such a definition of a nonequilibrium 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 database^{59}.
References
 1.
Chen, G. Nanoscale Energy Transport and Conversion (Oxford University Press, 2005).
 2.
Volz, S. et al. Nanophononics: state of the art and perspectives. Eur. Phys. J. B 89, 15 (2016).
 3.
Li, N. et al. Colloquium: Phononics: manipulating heat flow with electronic analogs and beyond. Rev. Mod. Phys. 84, 1045–1066 (2012).
 4.
Luo, T. & Chen, G. Nanoscale heat transfer—from computation to experiment. Phys. Chem. Chem. Phys. 15, 3389–3412 (2013).
 5.
Cahill, D. G. et al. Nanoscale thermal transport. ii. 2003–2012. Appl. Phys. Rev. 1, 011305 (2014).
 6.
Siemens, M. E. et al. Quasiballistic thermal transport from nanoscale interfaces observed using ultrafast coherent soft xray beams. Nat. Mater. 9, 26 (2010).
 7.
Minnich, J. et al. Thermal conductivity spectroscopy technique to measure phonon mean free paths. Phys. Rev. Lett. 107, 095901 (2011).
 8.
Johnson, J. A. et al. Direct measurement of roomtemperature nondiffusive thermal transport over micron distances in a silicon membrane. Phys. Rev. Lett. 110, 025901 (2013).
 9.
HoogeboomPot, K. M. et al. A new regime of nanoscale thermal transport: collective diffusion increases dissipation efficiency. Proc. Natl Acad. Sci. USA 112, 4846 (2015).
 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).
 11.
Frazer, T. D. et al. Engineering nanoscale thermal transport: Size and spacingdependent cooling of nanostructures. Phys. Rev. Appl. 11, 024042 (2019).
 12.
Guyer, R. A. & Krumhansl, J. A. Thermal conductivity, second sound, and phonon hydrodynamic phenomena in nonmetallic crystals. Phys. Rev. 148, 778–788 (1966).
 13.
Beck, H., Meier, P. & Thellun, A. Phonon hydrodynamics in solids. Phys. Stat. Sol. a 24, 11–63 (1974).
 14.
Lee, S., Broido, D., Esfarjani, K. & Chen, G. Hydrodynamic phonon transport in suspended graphene. Nat. Commun. 6, 6290 (2015).
 15.
Ding, Z. et al. Phonon hydrodynamic heat conduction and knudsen minimum in graphite. Nano Lett. 18, 638–649 (2018).
 16.
Cepellotti, A. et al. Phonon hydrodynamics in twodimensional materials. Nat. Commun. 6, 1–7 (2015).
 17.
Li, X. & Lee, S. Crossover of ballistic, hydrodynamic, and diffusive phonon transport in suspended graphene. Phys. Rev. B 99, 085202 (2019).
 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).
 19.
Huberman, S. et al. Observation of second sound in graphite at temperatures above 100 K. Science 364, 375 (2019).
 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).
 21.
Gandolfi, M., Giannetti, C. & Banfi, F. Temperonic crystal: a superlattice for temperature waves in graphene. Phys. Rev. Lett. 125, 265901 (2020).
 22.
Cepellotti, A. & Marzari, N. Transport waves as crystal excitations. Phys. Rev. Mater. 1, 045406 (2017).
 23.
Torres, P., Alvarez, F. X., Cartoixà, X. & Rurali, R. Thermal conductivity and phonon hydrodynamics in transition metal dichalcogenides from firstprinciples. 2D Mater. 6, 035002 (2019).
 24.
Machida, Y., Matsumoto, N., Isono, T. & Behnia, K. Phonon hydrodynamics and ultrahigh–roomtemperature thermal conductivity in thin graphite. Science 367, 309–312 (2020).
 25.
Beardo, A. et al. Observation of second sound in a rapidly varying temperature field in Ge. Sci. Adv. 7, eabg4677 (2021).
 26.
Melis, C., Fugallo, G. & Colombo, L. Room temperature second sound in cumulene. Phys. Chem. Chem. Phys.23, 15275–15281 (2021).
 27.
Simoncelli, M., Marzari, N. & Cepellotti, A. Generalization of Fourier’s law into viscous heat equations. Phys. Rev. X 10, 011019 (2020).
 28.
Gandolfi, M. et al. Emergent ultrafast phenomena in correlated oxides and heterostructures. Phys. Scr. 92, 034004 (2017).
 29.
Sobolev, S. Local nonequilibrium electron transport in metals after femtosecond laser pulses: a multitemperature hyperbolic model. Nanoscale Microscale Thermophys. Eng. 1–13 (2021).
 30.
Hartnoll, S. A. Theory of universal incoherent metallic transport. Nat. Phys. 11, 54–61 (2015).
 31.
Lee, S. et al. Anomalously low electronic thermal conductivity in metallic vanadium dioxide. Science 355, 371–374 (2017).
 32.
Tokura, Y., Kawasaki, M. & Nagaosa, N. Emergent functions of quantum materials. Nat. Phys. 13, 1056–1068 (2017).
 33.
Basov, D. N., Averitt, N. D. & Hsieh, D. Towards properties on demand in quantum materials. Nat. Mat. 16, 1077–1088 (2017).
 34.
Moyer, J. A., Eaton, C. & EngelHerbert, R. Highly conductive srvo3 as a bottom electrode for functional perovskite oxides. Adv. Mater. 25, 3578–3582 (2013).
 35.
Zhong, Z. et al. Electronics with correlated oxides: SrVO_{3}/SrTiO_{3} as a mott transistor. Phys. Rev. Lett. 114, 246401 (2015).
 36.
Zhang, L. et al. Correlated metals as transparent conductors. Nat. Mater. 15, 204–210 (2015).
 37.
OrdonezMiranda, J., Ezzahri, Y., Joulain, K., Drevillon, J. & AlvaradoGil, J. J. Modeling of the electrical conductivity, thermal conductivity, and specific heat capacity of VO_{2}. Phys. Rev. B 98, 075144 (2018).
 38.
Cesarini, G. et al. Quantitative evaluation of emission properties and thermal hysteresis in the midinfrared for a single thin film of vanadium dioxide on a silicon substrate. Int. J. Therm. Sci. 146, 106061 (2019).
 39.
BenAbdallah, P. & Biehs, S.A. Nearfield thermal transistor. Phys. Rev. Lett. 112, 044301 (2014).
 40.
OrdonezMiranda, J., Ezzahri, Y., TiburcioMoreno, J. A., Joulain, K. & Drevillon, J. Radiative thermal memristor. Phys. Rev. Lett. 123, 025901 (2019).
 41.
Fabrizio M. The OutofEquilibrium TimeDependent 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) https://doi.org/10.1007/9789400749849_16.
 42.
Mazza, G., Amaricci, A., Capone, M. & Fabrizio, M. Electronic transport and dynamics in correlated heterostructures. Phys. Rev. B 91, 195124 (2015).
 43.
Sandri, M., Capone, M. & Fabrizio, M. Finitetemperature gutzwiller approximation and the phase diagram of a toy model for v_{2}o_{3}. Phys. Rev. B 87, 205108 (2013).
 44.
Lanatà, N., Deng, X. & Kotliar, G. Finitetemperature gutzwiller approximation from the timedependent variational principle. Phys. Rev. B 92, 081108 (2015).
 45.
Georges, A., Kotliar, G., Krauth, W. & Rozenberg, M. J. Dynamical meanfield theory of strongly correlated fermion systems and the limit of infinite dimensions. Rev. Mod. Phys. 68, 13–125 (1996).
 46.
Mingo, N. & Broido, D. A. Carbon nanotube ballistic thermal conductance and its limits. Phys. Rev. Lett. 95, 096105 (2005).
 47.
Muñoz, E., Lu, J. & Yakobson, B. I. Ballistic thermal conductance of graphene ribbons. Nano Lett. 10, 1652–1656 (2010).
 48.
Bae, M.H. et al. Ballistic to diffusive crossover of heat flow in graphene ribbons. Nat. Commun. 4, 1734 (2013).
 49.
Caddeo, C. et al. Thermal boundary resistance from transient nanocalorimetry: a multiscale modeling approach. Phys. Rev. B 95, 085306 (2017).
 50.
Tzou, D. Y. Macroto Microscale Heat Transfer: The Lagging Behavior (John Wiley & Sons, 2014).
 51.
Aizaki, S. et al. Selfenergy on the low to highenergy electronic structure of correlated metal SrVO_{3}. Phys. Rev. Lett. 109, 056401 (2012).
 52.
Inoue, I. H., Goto, O., Makino, H., Hussey, N. E. & Ishikawa, M. Bandwidth control in a perovskitetype 3d^{1}correlated metal Ca_{1−x}Sr_{x}VO_{3}. I. Evolution of the electronic properties and effective mass. Phys. Rev. B 58, 4372–4383 (1998).
 53.
Karrasch, C., Kennes, D. M. & HeidrichMeisner, F. Thermal conductivity of the onedimensional FermiHubbard model. Phys. Rev. Lett. 117, 116401 (2016).
 54.
Giannetti, C. et al. Ultrafast optical spectroscopy of strongly correlated materials and hightemperature superconductors: a nonequilibrium approach. Adv. Phys. 65, 58 (2016).
 55.
Hwang, H. Y. et al. Emergent phenomena at oxide interfaces. Nat. Mater. 11, 103–113 (2012).
 56.
Cao, Y. et al. Correlated insulator behaviour at halffilling in magicangle graphene superlattices. Nature 556, 80?84 (2018).
 57.
Seyler, K. L. et al. Signatures of moirétrapped valley excitons in MoSe_{2}/WSe_{2} heterobilayers. Nature 567, 66–70 (2019).
 58.
Lu, X. et al. Superconductors, orbital magnets and correlated states in magicangle bilayer graphene. Nature 574, 653–657 (2019).
 59.
Mazza, G., Gandolfi, M., Capone, M., Banfi, F. & Giannetti, C. Dataset “Ultrafast thermal dynamics in strongly correlated systems". YARETA (2021) https://doi.org/10/gnhcw9.
Acknowledgements
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 (ERC319286QMAC). 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 (ANR16IDEX0005) 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
Affiliations
Contributions
G.M. developed the theoretical framework and performed the numerical study of the nonequilibrium 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
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 http://creativecommons.org/licenses/by/4.0/.
About this article
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). https://doi.org/10.1038/s41467021270812
Received:
Accepted:
Published:
DOI: https://doi.org/10.1038/s41467021270812
Comments
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.