1. Introduction

Traditionally, heat conduction is described as a diffusive process in terms of Fourier’s law, where the energy flux is proportional to the local thermal gradient. In nonmetallic crystalline solids, this description is valid only if phonons—the quantized vibrational modes that carry heat—locally and instantaneously perceive the inhomogeneity in the thermal field. Microscopically, this corresponds to situations where the characteristic size and timescale under consideration are significantly larger than the phonon mean free paths and scattering times, respectively^[1,2]. Therefore, Fourier’s law is generally inadequate to describe the thermal response in current electronic devices, where the inclusion of thin multilayered structures and miniaturized components introduce characteristic scales in the nanometer range^[3,4], two orders of magnitude smaller than the average room-temperature phonon mean free path in standard semiconductors such as silicon^[5].

To develop advanced heat transport models beyond diffusion and including boundary conditions, the usual starting point is the linearized Boltzmann transport equation (BTE) for phonons,

(1)$\frac{\partial f_{\lambda }}{\partial t}+{\overrightarrow{v}}_{\lambda }\cdot \overrightarrow{\nabla }{f}_{\lambda }=C\left(f_{\lambda }\right)$

where f_λ is the phonon distribution function of the λ-phonon mode (λ stands for the wave-vector $\overrightarrow{k}$ and polarization s), ${\overrightarrow{v}}_{\lambda }$ is the phonon group velocity, and C is the linearized collision operator, describing phonon interactions. While accessing solutions of this equation including arbitrary complex boundary conditions is not feasible, some approximations can be performed to predict the influence of nanoscale geometries or rapid excitations. The most common hypothesis to simplify operator C is the relaxation time approximation (RTA), which imposes that each phonon mode relaxes independently back to equilibrium with its own characteristic timescale, leading to completely uncorrelated evolution of the phonon mode populations. The RTA can successfully approximate the bulk thermal conductivity in crystals such as silicon, and has been extrapolated to describe heat transport at the nanoscale as a quasi-ballistic process^[6-9]. This picture provides a phenomenological interpretation of the reduction of the apparent thermal conductivity in the presence of nanoscale constraints, and has given rise to the notions of phonon suppression^[10], cumulative thermal conductivity^[5], mean free path spectroscopy^[9], phonon Lévy flights, and superdiffusion^[11]. Nevertheless, this approach neglects the influence of momentum-conserving phonon interactions^[12,13] and requires careful consideration of thermodynamic constraints^[14,15], which are generally not satisfied by the RTA. Moreover, the adequacy of using bulk phonon properties within the RTA-based quasi-ballistic picture to describe the influence of nanoscale confinement on heat flow remains under debate^[16-18].

Phonon hydrodynamics constitutes an alternative route to treat the BTE with special attention to the conservation laws and phonon mode correlations at the nanoscale and predicts the emergence of collective phonon evolution in a range of non-equilibrium conditions^[19]. Accordingly, phonon evolution is best understood not as a random walk of independent particles, but as a correlated phonon gas that can sustain transport phenomena akin to fluid-like viscosity or thermal waves. This description bridges microscopic (particle-based) and macroscopic (continuum) descriptions, providing new insights and a unified framework for modeling thermal transport. The hydrodynamic and RTA-based interpretations of nanoscale heat transport differ fundamentally in their description of phonon dynamics and lead to contrasting predictions, including the number of intrinsic phonon length scales expected to directly manifest in experiments^[18]. Determining which approach provides a deeper understanding of phonon evolution at the nanoscale remains an outstanding question and is essential for incorporating the key transport mechanisms beyond diffusion into predictive modeling tools. Ultimately, adequate theoretical treatment of the BTE for phonons will influence a broad range of engineering applications, including improved heat dissipation in nanoelectronics and enhanced performance in thermoelectric devices.

Extensive literature reviewing the topic of phonon hydrodynamics and its applications is available, including but not limited to previous studies^[19-23]. This manuscript is a complementary review aimed at unifying the different perspectives and flavors that this theoretical framework has adopted historically. We focus on illustrating the general applicability of the hydrodynamic approach as a tool to understand phonon evolution out of equilibrium, from both theoretical and experimental perspectives, and we critically revise restrictive applicability criteria that are commonly assumed. First, we review the diversity of routes to derive hydrodynamic-like equations from microscopic grounds and the BTE. Then, we discuss the key features and general thermodynamic robustness of the resulting phonon hydrodynamic models. Finally, we briefly but comprehensively review the experiments where hydrodynamic descriptions have successfully provided interpretations of non-diffusive effects. In this last section, we emphasize studies in which all the parameters of the hydrodynamic equations are quantified using ab initio calculations associated with specific averages over the phonon properties, as derived from the BTE.

2. Microscopic Foundation

The deviation from equilibrium of the phonon distribution can be expanded in terms of an arbitrary perturbation,

(2)$f_{\lambda }=f_{\lambda }^{\mathrm{eq}}-{\Phi }_{\lambda }\frac{\partial f_{\lambda }^{\mathrm{eq}}}{\partial {\epsilon }_{\lambda }}$

where ε_k is the phonon energy and Φ_λ generally depends on a set of macroscopic variables, such as the thermal gradient, weighted by mode-dependent coefficients. From the point of view of the variational principle, the solution of the linearized BTE, Eq. (1), for a given non-equilibrium situation corresponds to the phonon distribution function that optimizes the entropy production rate while accommodating all the boundary conditions and conservation laws. In other words, the BTE can be formulated as an optimization problem for a function directly related to entropy production and that only depends on Φ_λ and the collision operator. This interpretation was carefully discussed by J. M. Ziman in Chapters 7.7 and 7.8 of the previous study^[1]. Accordingly, it is crucial to select adequate macroscopic variables to expand Φ_λ depending on the non-equilibrium conditions under consideration. Otherwise, entropy production is optimized within a sub-space of solutions that does not accommodate sufficient complexity, and the resulting solution can be arbitrarily far from the true solution. This is particularly relevant in complex non-equilibrium situations, such as those in the presence of non-homogeneous heat flux profiles at the nanoscale. In light of this general variational principle, phonon hydrodynamic theory is a method to obtain refined heat transport models from specific macroscopic characterizations of Φ_λ, which are specifically designed to capture the main non-equilibrium features emerging at the nanoscale. As discussed in the remainder of this article, the analogy with fluid dynamics is pertinent in this context because an adequate choice of the expansion for Φ_λ in terms of slowly-evolving thermodynamic variables generally leads to hydrodynamic-like heat transport equations.

The development and implementation of the hydrodynamic interpretation of phonon transport has a long and rich history. Originally, hydrodynamic behavior was predicted to manifest at ultra-low temperatures in specific materials where momentum-conserving (Normal) phonon collisions are much more frequent than resistive ones^[24-28]. Not only are Normal phonon-phonon collisions analogous to molecular collisions in fluids, but the first derivation of a hydrodynamic heat transport equation resembling the Navier-Stokes equation was derived from the BTE using the dominance of Normal collisions as a central hypothesis^[26]. More recently, a substantial body of work has shown that the abundance of Normal collisions is a sufficient but not necessary condition for the emergence of phonon hydrodynamics^[29-38]. A more fundamental necessary condition for the emergence of hydrodynamic transport effects at the nanoscale is the comparatively slow relaxation of the phonon distribution perturbation associated with the heat flux, relative to higher-order nonequilibrium perturbations^[39-41]. This condition is not restricted to situations where most phonon interactions conserve the phonon momentum, but is generally satisfied at length or time scales comparable to the average mean free path or scattering time associated with resistive phonon interactions. Phenomenologically, this implies that the heat flux profile cannot locally or instantaneously align with the thermal gradient, thus leading to non-local and memory transport effects^[42]. This view, schematically illustrated in Figure 1, extends the significance of hydrodynamic descriptions to a wide range of semiconductor materials and temperatures. The applicability of current hydrodynamic models, however, has been shown to be limited to moderately small Knudsen numbers, which prevents predictive modeling at arbitrary small length and timescales^[30,43].

Display Full Size

Figure 1. Illustration of non-equilibrium phonon evolution. We consider a nanoscale constriction connecting a heat source terminal subjected to a pulsed injection of thermal energy and a heat sink terminal at constant temperature. (a,b) Schematic representation of a possible non-equilibrium phonon distribution function in terms of the amplitude of its statistical moments at the position denoted by the red star. The insets show the phonon population as a function of the wave-vector direction. At time scales much smaller than the heat flux relaxation time τ after pulsed excitation, the deviation from equilibrium is arbitrarily complex and the distribution function accommodates high-order statistical moments. At time scales comparable to τ, all the non-equilibrium statistical moments have already relaxed due to phonon scattering except the ones directly related to energy and heat flux. This intermediate non-equilibrium state, which is generally established before thermal equilibrium is recovered, displays hydrodynamic-like transport features that can be modeled using the GKE; (c) Representative heat flux profiles along the cross-section denoted by the blue line at the different time scales. The hydrodynamic signature manifests here as a viscous non-local correlation of the heat flux profile restricted to the heat flux boundary layers (denoted in gray color). GKE: Guyer-Krumhansl equation.

In this section, we briefly review some of the key derivations of hydrodynamic-like heat transport models from the BTE. We focus on the different expansions for the non-equilibrium phonon distribution assumed in each case, which limit the applicability of the resulting models to different scenarios, and highlight the generality of this methodology to interpret nondiffusive thermal transport effects.

2.1 Guyer and Krumhansl derivation and historic perspective

In the collective limit, where Normal phonon interactions are much more frequent than resistive ones, an arbitrary non-equilibrium phonon state rapidly decays toward the displaced phonon distribution—the maximum entropy state that accommodates the total crystalline momentum in the absence of boundaries—before ultimately decaying to equilibrium via momentum-relaxing resistive interactions^[24]. This displaced distribution reads,

(3)$f_{\lambda }^{\text{disp}}=\frac{1}{\exp \left[({\omega }_{\lambda }-{\overrightarrow{k}}_{\lambda }\cdot \overrightarrow{u})ħ/\left(k_{B}T\right)\right]-1}$

where ω_λ is the phonon frequency, T is the temperature, k_B is the Boltzmann constant, and $\overrightarrow{u}$ is the mode-independent drift velocity. This distribution collapses to the equilibrium one for $\overrightarrow{u}$ = 0. Thus, the presence of a drift velocity is a mesoscopic signature of a strong out-of-equilibrium correlation between all the phonon modes, which enables establishing a microscopic analogy with standard fluid mechanics^[25,27]. This insight anticipated the emergence of strong hydrodynamic effects^[44-47], such as the Knudsen minimum, in materials and temperatures that can accommodate the collective limit, such as graphene at 100 K^[48-49].

To investigate the implications of collective phonon evolution, a hydrodynamic heat transport equation was originally derived from the BTE in the seminal work of Guyer and Krumhansl^[26]. They first decomposed the collision operator into terms of Normal and resistive interactions, C = N + R. Instead of the phonon-mode basis, they expanded the distribution function in terms of the eigenvectors of the Normal collision operator N. Some of these eigenvectors, which are specific linear combinations of phonon modes, can be related to key thermodynamic quantities such as the energy and the crystalline momentum. When formulated in this basis, the BTE can be substantially simplified by enforcing thermodynamic constraints. For instance, energy conservation implies that the first N-eigenvector does not relax under the action of either R or N, while momentum conservation imposes that the three N-eigenvectors associated with the crystal momentum remain invariant under N. In this picture, the bulk thermal conductivity can be decomposed in terms of collective and kinetic contributions, thus providing physical intuition about the phonon transport mechanisms^[26,50-52]. Most importantly, the BTE can be projected into the form of a hydrodynamic-like heat transport equation by assuming the dominance of Normal collisions $N\gg R$ (collective limit)^[26]. The resulting equation, known as the Guyer-Krumhansl equation (GKE), is structurally similar to the Navier-Stokes equation of classical fluids in the laminar regime, and relates the heat flux $\overrightarrow{q}$ and the temperature T,

(4)$\overrightarrow{q}+\tau \frac{\partial \overrightarrow{q}}{\partial t}=-\kappa \overrightarrow{\nabla }T+{\ell }^2({\nabla }^2\overrightarrow{q}+\zeta \overrightarrow{\nabla }\overrightarrow{\nabla }\cdot \overrightarrow{q})$

where κ is the bulk thermal conductivity, τ is the heat flux relaxation time, ℓ is the non-local length, and ζ is a dimensionless viscosity coefficient, all resulting from specific averages of the phonon group velocities, specific heats, and relaxation times.

The GKE, in combination with appropriate boundary conditions^[30,53] and the energy conservation equation ($-C_v\partial T/\partial t+Q=\nabla \cdot \overline{q}$, with C_v being the volumetric heat capacity and Q an external heat source), generalizes the heat diffusion equation by capturing collective phonon phenomena in the form of memory and non-local effects^{[19,42,46,54,55]}. These include a reduction of the apparent conductivity relative to the bulk diffusive limit under nanoscale boundary constraints or inhomogeneous excitations, and the emergence of second sound—a wave-like evolution of the temperature field at small timescales.

The Guyer and Krumhansl framework represents the starting point for the implementation of hydrodynamic heat transport modeling. According to the original derivation, this treatment is restricted to the collective limit, where in principle the GKE can be extended to very small length scales due to the small Knudsen number associated with abundant Normal collisions. Early experimental efforts identified signatures of strong hydrodynamic effects in this limit at ultra-low temperatures^[56,57], and more recent ones in graphite at higher temperatures, including fully developed Poiseuille profiles^[58,59] and damped thermal waves^[60,61].

2.2 Perturbation expansion

More recently, the regime of applicability of Eq. (4) has been extended into weak hydrodynamic regimes beyond the collective limit, analogous to fluid dynamics in rarefied gases^[62]. The first modern derivation of phonon hydrodynamics was introduced in the previous work^[30]. The derivation considers the BTE with a single scattering time (gray-model) and a single group velocity (Debye approximation), and employs a Perturbation Expansion around a nonequilibrium solution characterized by four phonon moments,

(5)$f_{\lambda }=f_{\lambda }^{\mathrm{eq}}+\frac{3}{C_v{v}_g^2}\frac{\partial f^{eq}}{\partial T}\overrightarrow{q}\cdot {\overrightarrow{v}}_g+\frac{\overline{\tau }}{C_v}\frac{\partial f^{eq}}{\partial T}\overrightarrow{\nabla }\cdot \overrightarrow{q}-\frac{3\overline{\tau }}{C_v{v}_g^2}{\overrightarrow{v}}_g{\overrightarrow{v}}_g:\overrightarrow{\nabla }\overrightarrow{q}$

where $f_{\lambda }^{\mathrm{eq}}$ is the equilibrium distribution at the local temperature T, $\overline{\tau }$ is an average phonon scattering time, ${\overrightarrow{v}}_g$ is the average group velocity, and C_v is the volumetric heat capacity.

The heat transport equation resulting from the BTE and the distribution in Eq. (5) is formally equivalent to the GKE, Eq. (4). Thus, this modern derivation expands the applicability of phonon hydrodynamics beyond the traditional low-temperature regime, and highlights the significance of phonon hydrodynamics in nanoscale crystalline semiconductors even at ambient conditions in semiconductors such as silicon. In addition, the nonequilibrium distribution function in Eq. (5) can be used to derive consistent boundary conditions by imposing microscopic energy and heat flux balance constraints^[30]. The resulting boundary conditions are formulated at the level of integrated variables—temperature and heat flux—rather than individual phonon modes, thereby enabling implementation in complex geometries and nanostructures, including interfaces^{[53,54,63,64]}.

It is worth stressing that the Perturbation Expansion assumes an average relaxation time for all phonon modes, as reflected in Eq. (5), and does not provide a microscopic expression to quantify the non-local length ℓ and the flux relaxation time τ using first principles calculations of the spectral phonon transport properties^[5,65,66].

2.3 Flux derivatives formalism

An alternative derivation of Eq. (4) from the BTE was later presented in the studies^[32,67], which introduced a formalism that systematically solves the BTE for arbitrary materials and temperatures at moderately small Knudsen numbers. The phonon distribution is assumed to be well approximated by an expansion in terms of the heat flux and its spatial and time derivatives, which are treated as independent variables,

(6)$f_{\lambda }=f_{\lambda }^{\mathrm{eq}}+{\overrightarrow{\beta }}_{\lambda }\cdot \overrightarrow{q}+{\overrightarrow{\gamma }}_{\lambda }\cdot \frac{\partial \overrightarrow{q}}{\partial t}+{\bar{\overline{G}}}_{\lambda }:\overrightarrow{\nabla }\overrightarrow{q}$

and the BTE is used to determine the mode-dependent prefactors weighting each term in the expansion. Once this distribution function is fully characterized using ab initio inputs, the energy and momentum projections of the BTE can be used to recover the GKE. Crucially, this approach provides explicit expressions for the transport parameters κ, ℓ, τ, and ζ as averages over the spectral properties of the complete phonon distribution, thereby going beyond the gray-model approximation, and accounting for realistic phonon dispersion relations beyond the Debye approximation^[32]. The transport parameters values from ab initio calculations for graphene, silicon, germanium, and diamond across a range of temperatures are provided in the studies^[32,52].

The Flux Derivatives Formalism generalizes the approach of the previous work^[68], in which the perturbation is expanded in terms of a constant temperature gradient, to the nanoscale, where inhomogeneous and rapidly varying perturbations occur. This framework implicitly assumes that, at the nanoscale, arbitrary non-equilibrium distribution functions rapidly decay to the intermediate non-equilibrium state described by Eq. (6), rather than directly relaxing back to equilibrium. As illustrated in Figure 1, this hierarchical relaxation of the statistical moments of the phonon distribution is the key microscopic reason behind the manifestation of hydrodynamic effects at small scales. From a thermodynamic perspective, the choice of the heat flux as the key magnitude to characterize the non-equilibrium state is motivated by its slower relaxation compared to other non-equilibrium features under spontaneous phonon dynamics and interactions at moderately small Knudsen numbers. Microscopically, this approximation captures the limited phase space for phonon evolution restricted by the boundary conditions and conservation laws.

Importantly, this approach does not implicitly assume the gray-model approximation for the phonon scattering landscape. Instead, the microscopic expressions for the transport parameters depend on the phonon properties of all modes, and can be calculated via iterative methods by accounting for the full collision operator^[32]. Therefore, the inclusion of a single characteristic non-local length ℓ and relaxation time τ in the GKE to model the influence of boundaries and heat sources on thermal evolution is not an approximation, but rather a defining feature of phonon hydrodynamic behavior^[37,69]. This characteristic distinguishes the hydrodynamic framework from RTA-based descriptions, in which each phonon mode interacts with boundaries independently and the multiscale mean free path spectrum explicitly manifests^[18].

Finally, we note that non-local and memory effects at the nanoscale compatible with the GKE have recently been predicted using alternative methodologies^[38,41,70], wherein the GKE can be obtained as a limiting case of more general transport equations. Future work should aim to reconcile the different derivations from first principles, and identify the significance of higher-order and non-linear transport effects beyond the ones captured by the GKE^[38,71-73].

2.4 Viscous heat equation

As illustrated by the analysis of Guyer and Krumhansl, the eigenvectors of the collision operator constitute a useful basis for studying the BTE for phonons. In contrast to the phonon mode basis, the C-eigenvectors decay exponentially toward equilibrium with well-defined relaxation times. The only exception is the eigenvector associated with the energy density, which does not relax under collisions. Since different C-eigenvectors remain uncorrelated under phonon interactions, this change of basis implicitly captures collective phonon behavior while reducing the collision operator in the BTE to a diagonal form. This advantage was exploited in a series of seminal works by Robert Hardy and coauthors^[25,39], which distinguished between odd eigenvectors with a non-null contribution to the heat flux and even eigenvectors, which represent symmetric non-equilibrium perturbations and therefore do not contribute to the net heat flux. Under spatially homogeneous perturbations, the relaxation time of the odd C-eigenvectors is therefore directly related to the timescale for the decay of energy currents in the system^[41].

This perspective was later adopted in recent work^[74], where the decay times of the C-eigenvectors, denoted as relaxons, were quantified from first-principles calculations for different materials. In addition, this approach enabled solving the BTE in the presence of a homogeneous thermal gradient and computing the bulk thermal conductivity without requiring iterative methods. However, despite the usefulness of the relaxon basis to model simple non-equilibrium situations, relaxons are not eigenvectors of the drift operator, as they evolve into other relaxons due to phonon propagation. Therefore, describing complex non-equilibrium situations beyond homogeneous thermal gradients using the relaxon basis has remained elusive, requiring further approximations to account for boundary effects^[75].

To address these limitations, a mesoscopic hydrodynamic-like heat equation was later derived from the BTE using the relaxon basis^[76]. The key assumption here is to expand the non-equilibrium state in terms of linear perturbations of the displaced distribution, Eq. (3):

(7)$f_{\lambda }=f_{\lambda }^{\mathrm{eq}}+\frac{\partial f_{\lambda }^{\mathrm{disp}}}{\partial T}∆T+\frac{\partial f_{\lambda }^{\mathrm{disp}}}{\partial \overrightarrow{u}}\cdot \overrightarrow{u}+∆f_{\lambda }$

where ∆T refers to the local temperature deviation from global equilibrium, $\overrightarrow{u}$ is the drift velocity, and ∆f_λ is any additional deviation from equilibrium. Accordingly, the key mesoscopic variable is the drift velocity $\overrightarrow{u}$, which is a meaningful magnitude in situations where Normal collisions dominate and the phonon momentum is mechanically conserved over long timescales^[24]. Hence, this approach is appropriate in the collective limit, where $\overrightarrow{u}$ is a slowly-evolving variable and Eq. (7) can approximate a long-lived intermediate state prior to the final relaxation toward equilibrium. This is useful in specific materials like graphene at low temperatures^[48,49], but inadequate to describe standard semiconductors such as silicon, where the drift velocity is not a meaningful mesoscopic variable. Out of the collective regime, this approach essentially reduces to the standard local equilibrium approximation associated with diffusive transport^[38].

The mesoscopic transport equation consistent with Eq. (7) and the BTE contains terms analogous to viscous effects in fluids, but exhibits significant qualitative differences compared to the GKE. The Viscous Heat Equation predicts the existence of a diffusive heat flux background that obeys Fourier’s law with the bulk conductivity, combined with an additional heat flux contribution that displays hydrodynamic behavior^[76]. In contrast to the GKE, this equation thus predicts that hydrodynamic effects induce an increase in apparent thermal conductivity relative to the diffusive background, rather than a reduction. Consequently, this formalism does not attempt to capture the non-local reduction of the heat current caused by the presence of boundaries, but instead aims to model strong hydrodynamic effects in the collective limit, such as the Knudsen minimum^[45]. To model the reduction in conductivity with respect to the diffusive limit generally observed at the nanoscale, the transport parameters of the viscous heat equation must be corrected by using reduced phonon mode scattering times^[76].

2.5 Direct BTE solvers

The construction of mesoscopic hydrodynamic-like heat transport equations is not the only possibility to describe heat flux correlations and collective phonon effects. Phonon hydrodynamics can also be studied by directly solving the linearized BTE, Eq. (1), with a collision operator fully characterized from ab initio inputs. Since this approach seeks to capture the full microscopic complexity of phonon dynamics, it is worth emphasizing a few general considerations to ensure a physically meaningful description. First and foremost, it is necessary to avoid the RTA, which explicitly suppresses collective behavior and correlations between different phonon modes in direct microscopic solvers. Second, the macroscopic variables to expand the phonon distribution in Eq. (2) need to be carefully selected depending on the non-equilibrium conditions under consideration. This is particularly important beyond local equilibrium and in the presence of inhomogeneous heat flux profiles, where the local thermal gradient is not sufficient to characterize the non-equilibrium constraints. Third, low-frequency phonon modes play an important role in determining the hydrodynamic length and timescales. Hence, direct BTE solutions must be meticulously converged in terms of the phonon wave-vector space discretization.

In recent years, a variety of methods have been developed to solve the linearized BTE beyond the RTA. Under homogeneous thermal gradients, the distribution function and the associated bulk thermal conductivity can be obtained by iterative methods^[65,77]. These iterative solutions demonstrate that the displaced distribution described by Eq. (3) can be approached in some materials such as graphene even at relatively high temperatures^[48,49,52]. This is an indication that some materials and temperatures can host the collective limit and the associated strong hydrodynamic regime, leading to pronounced non-diffusive responses in the presence of boundaries^[59] or rapid excitations^[60]. Nevertheless, under complex non-equilibrium conditions and nanoscale boundary constraints, iterative solutions of the BTE are generally not feasible. To address this, direct solvers of the BTE usually simplify the collision operator. Beyond the RTA, the Callaway approximation is commonly employed^[12],

(8)$C\left(f_{\lambda }\right)=-\frac{f_{\lambda }-f_{\lambda }^{\mathrm{eq}}}{{\tau }_{\lambda }^R}-\frac{f_{\lambda }-f_{\lambda }^{\mathrm{disp}}}{{\tau }_{\lambda }^N}$

where ${\tau }_{\lambda }^R$ is the average scattering time of a specific phonon mode due to resistive interactions, such as Umklapp or isotopic scattering, and ${\tau }_{\lambda }^N$ is the scattering time considering only Normal interactions. Accordingly, momentum-conserving collisions can only relax the nonequilibrium state toward the displaced distribution, Eq. (3). This approximation is commonly employed in numerical BTE solvers to account for boundaries and predict non-diffusive transport effects^{[13,31,78,79]}, and has recently been used to recover the GKE derivation in the collective limit^[73].

However, Eq. (8) does not inherently capture the complexity of intermediate non-equilibrium states established at the nanoscale, such as those described by Eq. (5) and Eq. (6), or Eq. (7), which can cause discrepancies with more sophisticated solvers using the full form of the collision operator^[80]. Beyond Callaway’s approximation, the BTE can be solved using a deviational Monte Carlo approach^[16,61] or the Green’s function approach^[81]. These techniques have been successfully implemented to describe the emergence of second sound under ring-shaped^[61] and transient grating^[60] optical excitations in the absence of boundaries. Hence, these direct BTE solvers constitute alternative routes for predicting phonon hydrodynamic effects relative to mesoscopic heat transport equations such as the GKE.

3. Thermodynamic Foundation

The fundamental condition underlying the emergence of collective phonon behavior in the form of hydrodynamic effects is the slow relaxation of the non-equilibrium part of the distribution accommodating the heat flux^[39]. The implicit existence of this slowly-relaxing variable constitutes the common necessary condition across all approaches to derive phonon hydrodynamic models.

In formalisms where the heat transport equation is derived from the BTE by expanding the phonon distribution function in terms of the heat flux and related variables, such as the Perturbation Expansion or the Flux Derivatives Formalism, this condition is explicitly imposed by assuming that the non-equilibrium state mainly depends on the shape of the heat flux vector field. This directly implies that other forms of non-equilibrium perturbations relax much faster and can be neglected when approximating the phonon mode populations. As illustrated in Figure 1, this is a reasonable hypothesis at relatively short length or timescales, where the limited phase space for phonon interactions reveals a hierarchical decay of the non-equilibrium state^[1].

In formalisms assuming the collective limit, such as the original derivation of the GKE, the abundance of Normal collisions slows down the relaxation of the total crystal momentum compared to other potentially excited mesoscopic variables. Microscopically, this corresponds to a slow decay of the asymmetric features of the distribution function, which in turn entails a slow relaxation of the total heat flux. Consequently, the collective limit assumption constitutes a sufficient condition for the emergence of phonon hydrodynamics. It is important to emphasize, however, that this is not a necessary condition, particularly at the nanoscale where the heat flux does not evolve locally or instantaneously.

In summary, the slow relaxation of the heat flux implies non-local and memory effects on the thermal evolution, reminiscent of fluid-like behavior. These effects are captured in the GKE by introducing a single characteristic decay length ℓ and timescale τ, which directly quantify the relaxation dynamics of the heat flux. This is in stark contrast with RTA-based descriptions of nanoscale heat transport, where the full spectrum of decay timescales associated with each phonon mode manifests explicitly^[6,9]. Thus, the existence of a single interaction scale characterizing nanoscale dynamics represents the primary and most general signature of phonon hydrodynamics, rather than any specific microscopic criterion. This crucial aspect of hydrodynamic theory underlies the similarities between the GKE and alternative models based on the existence of a single interaction scale in time and space, such as the Dual Phase Lag model^[82,83].

Finally, we note that hydrodynamic equations derived from the BTE, such as the GKE, are fully consistent with the thermodynamic principles^[23,84]. The energy conservation constraint is explicitly satisfied in the formalisms described above. Moreover, the second law can also be reconciled with hydrodynamic-like equations by considering advanced non-equilibrium thermodynamic formulations. For instance, Extended Irreversible Thermodynamics^[85] proposes the generalization of the Gibbs equation by including dissipative fluxes as independent variables characterizing the local non-equilibrium state. This generalized thermodynamic relation can be combined with conservation laws to demonstrate that the GKE predicts a positive entropy production under arbitrary non-equilibrium conditions^[84]. Notably, the idea of using the heat flux as a key variable to describe the non-equilibrium thermodynamic state motivates the variational treatment of the phonon distribution in the Perturbation Expansion method or the Flux Derivatives Formalism described above.

4. Experiments and Phenomenology

Hydrodynamic phonon transport manifests in the form of viscous-like non-local correlations of the heat flux, or memory and thermal waves^[27,36,42]. The emergence of these effects requires that the relevant length and timescales under consideration approach the nonlocal length ℓ and the heat flux relaxation time τ, respectively. Hydrodynamic phenomena are therefore amplified at low temperatures, where scattering rates decrease, or in materials where resistive processes such as Umklapp, phonon-defect, isotope, or alloy scattering are weak, allowing the distribution function to deviate from local equilibrium and increasing ℓ and τ. The interplay of heat viscosity, memory, and diffusion at different scales, as predicted by the GKE under a simple non-equilibrium situation, is illustrated in Figure 2.

Display Full Size

Figure 2. Emergence and interplay of weak hydrodynamic transport effects at different length and time scales. We consider a pulsed excitation with a Gaussian-shape injecting heat during a time window T within a thin film with thickness L. Three different situations are shown both at early and late times t after excitation. (b, e) At length scales L and time scales T larger than the non-local length ℓ and the relaxation time τ, respectively, non-local and memory effects are negligible and the heat flux and temperature fields evolve following standard diffusion; (a, d) For L~ℓ and long timescales with respect to τ, the heat flux is reduced close to the film boundaries, which reduces the total energy current and delays thermal relaxation. This viscous effect tends to prevent the propagation of thermal waves; (c, f) Finally, if the excitation time T and the observation time t are comparable to τ, a fraction of the thermal energy propagates as a damped wave.

Experimental identification of hydrodynamic effects and quantitative comparison between theory and experiments requires characterizing the specific heat capacities and group velocities of all the phonon modes along with the phonon-phonon transition rates from first principles calculations^[66]. These inputs can be used either to directly model the experiment by solving the linearized BTE^[60,61,86], or to evaluate the transport parameters of mesoscopic hydrodynamic models, such as κ, ℓ and τ in the GKE^[30,35,87]. A systematic comparison of measured versus calculated GKE parameters in a variety of experiments is provided in the previous work^[32].

In this section, we review the experimental methods that have been used to detect phonon hydrodynamic behavior, and we identify the range of predictive applicability of hydrodynamic models with ab initio calculated inputs across different experiments. Overall, this section shows that the hydrodynamic perspective provides a unifying framework for predicting a variety of experimental configurations where Fourier’s law breaks down.

4.1 Stationary heat flow in nanostructures

4.1.1 Heat flux boundary layers

The energy current established in response to a thermal gradient across a semiconductor nanostructure is generally smaller than the diffusive prediction based on Fourier’s law^[88-90]. This reduction in apparent thermal conductivity is associated with the momentum-destroying phonon-boundary scattering events. For a free surface, the microscopic balance of the flux of heat flux according to the non-equilibrium distribution functions in Eq. (5) or Eq. (6) results in a slip boundary condition^[1,30,43]. This condition, combined with the GKE, predicts a reduced energy transfer rate and a non-uniform heat flux profile near the structure edges that accommodates a slip boundary flux, in analogy to fluid flow in rarefied gases^[42,55,91]. In general, resistive phonon-phonon interactions limit the heat flux correlations in space and time. Consequently, the curvature of the flux profile is confined to a region next to the boundaries, known as the heat flux boundary layer, with a characteristic size comparable to the non-local length (Figure 1 and Figure 2).

Although the flux profile is not directly accessible in experiments, predictions of the apparent thermal conductivity in steady-state can be directly compared with measurements. In contrast to standard fluid dynamics, the GKE with slip boundary flux predicts a linear dependence of the apparent conductivity on system size^[92]. Experiments show good agreement with the hydrodynamic predictions for moderate Knudsen numbers in standard semiconductors such as silicon^[30,43,93]. In particular, the reduced conductivity observed experimentally is accurately predicted by the GKE with parameter values obtained from ab initio calculations^[32] across different geometries and a wide range of temperatures for L down to 2ℓ (i.e., Kn = ℓ/L ≤ 1/2), where L is the smallest characteristic size of the nanostructure and ℓ is the non-local length in the GKE^[43]. Accordingly, the applicability of current hydrodynamic models in standard semiconductors is limited to relatively large inter-boundary distances (a few hundreds of nanometers at 300 K), which precludes the use of this description with ab initio inputs to model ultrathin layers present in nanoelectronic devices. Moreover, thermal conductivity measurements in alloy nanostructures reveal a gradual reduction of the conductivity as the device length scale is decreased over multiple orders of magnitude^[94]. This multiscale onset of non-diffusive behavior observed in alloys can be described using various implementations of the RTA^[9,15,95], and is incompatible with the hydrodynamic description based on a single interaction scale .

4.1.2 Poiseuille flow

In certain materials and temperatures, such as monolayer graphene around 100 K, the collective limit is approached because a large fraction of phonon-phonon interactions conserve momentum^[48,49,96]. In these conditions, phonon momentum is almost exclusively destroyed in the boundaries, and the heat flux profile becomes strongly correlated across the full nanostructure. Due to the rapid redistribution of phonon momentum via Normal collisions, the slip heat flux in the boundary decays to zero, and the GKE predicts a fully correlated quadratic heat flux profile, analogous to Poiseuille flow in standard fluid dynamics. This is a signature of strong collective phonon evolution, which is phenomenologically different from the presence of slip flux and boundary layers confining the viscous effect close to the system edges.

The manifestation of Poiseuille flow has traditionally been associated with specific trends in thermal conductivity as a function of temperature T^[27]. At sufficiently high temperatures, the mean free path of both resistive (Λ_R) and Normal (Λ_N) phonon-phonon interactions can become significantly smaller than the sample characteristic size L. In this limit of bulk diffusion, the thermal conductivity follows a T^-1 trend, reflecting the temperature-dependence of the resistive scattering time. In contrast, at ultra-low temperatures, both Λ_N and Λ_R can easily become much larger than the device size L. In this ballistic limit, the mean free paths are geometrically confined, and the thermal conductivity follows a T³ trend, reflecting the temperature dependence of the specific heat capacity. Deviations from these trends at intermediate temperatures can be associated with strong phonon hydrodynamic behavior. By increasing the temperature above the ballistic limit, a temperature dependence of the conductivity with an exponent larger than 3 can be established if Λ_N < L < Λ_R This is a manifestation of Poiseuille heat flow, which increases the maximum flux in the internal regions of the nanostructure relative to the ballistic limit. Finally, in between the Poiseuille regime and the diffusion limit, an exponential trend of conductivity versus temperature can also be observed. This effect is a signature of the Ziman limit^[27] and reflects the dominance of Normal collisions with respect to resistive ones in situations where the mean free paths are not geometrically confined. Indications of Poiseuille flow and the Ziman limit regimes have been observed in some materials, including isotopically purified graphite^[59], black phosphorus^[58], and even sapphire^[97]. Additionally, in the collective limit, the complex interplay between phonon-boundary and phonon-phonon interactions can induce a minimum in the thermal conductivity as a function of the device size, which is known as the Knudsen minimum^[45,98]. However, direct experimental evidence for this effect is scarce^[22]. The particular temperature and size dependence of the conductivity in the collective limit has been successfully described using direct solutions of the BTE under Callaway’s approximation^[13]. To precisely capture collective phonon behavior within this microscopic framework, careful consideration of isotope scattering, surface roughness and grain-boundary effects is required^[86,99].

At this point, it is worth stressing that the manifestation of hydrodynamic behavior is not limited to materials and temperatures that accommodate the Poiseuille flow and the Ziman limit. Viscous heat flux correlations generally emerge across a broader range of materials, temperatures, and nanoscale conditions, at length and time scales comparable to the resistive mean free path and scattering times. Accordingly, we emphasize that the relevance of including hydrodynamic corrections to the diffusive transport equation is not limited to the collective limit.

4.2 Nanostructured heat sources

A crucial case of study is the process of heat release from nanoscale structures toward a semiconductor substrate^[100-102]. The heat flux profile underneath nanoscale interfaces cannot immediately and locally align with the thermal gradient due to the lack of resistive phonon collisions, leading to inefficient thermal energy evacuation compared to the diffusive prediction. The resulting thermal relaxation of the nanoscale heat sources cannot be described using an effective Fourier’s law with a fitted conductivity or a fitted interfacial thermal resistance^[87,101]. From the hydrodynamic perspective, this is a consequence of the emergence of heat vorticity near the corners of the interfaces, where $\nabla \times \overrightarrow{q}\ne 0$ and $\overrightarrow{q}$ is not parallel to $-\nabla T$^[31,37,103].

According to the GKE, the regions that accommodate vorticity display a well-defined characteristic size ℓ (cf. Figure 3). This provides a unified explanation for the transition from diffusive to non-diffusive behavior by reducing the characteristic size of the heat source, as has been observed using different experimental techniques. First, the GKE with ab initio parameters^[32] accurately predicts the stationary thermoreflectance profile measured around electrically-excited gold nanostructures with different sizes and shapes on top of a silicon substrate^[40,103]. In contrast, the full temperature profile cannot be reproduced using Fourier’s law even with effective parameters^[101,103]. Second, the same hydrodynamic modeling and parameter values also capture the thermoelastic relaxation of nanoscale heat sources on silicon following pulsed optical excitation^[87]. Remarkably, the temperature decay in these experiments displays a double-exponential function form, with the first time scale associated to the interfacial thermal resistance and the second to the non-local length scale. Hence, these results provide direct evidence of a single characteristic phonon interaction scale, strongly supporting the occurrence of hydrodynamic phonon evolution in silicon at room temperature^[87].

Display Full Size

Figure 3. Illustration of the temperature and heat flux profiles in a metallic nanostructure releasing heat toward a semiconductor substrate. We consider a stationary situation where a constant and homogeneous power density is injected in the nanostructures. The left side of the figures show a hydrodynamic shape as predicted by the GKE, which is representative of situations where the interface width is comparable to the non-local length ℓ. The right side shows a diffusive shape corresponding to an interface width much larger than ℓ. GKE: Guyer-Krumhansl equation.

Although hydrodynamic modeling shows general agreement with experiments for nanoscale heat sources with characteristic sizes down to tens of nanometers, it has been shown to fail where the spacing between adjacent heat sources is smaller than ~2ℓ. In these cases, non-local heat flux correlations and vorticity effects appear to be geometrically constrained, leading to faster thermal relaxation compared to isolated heat sources^[104]. This counter-intuitive effect is not captured by the GKE with ab initio inputs, but can be reproduced by geometrically constraining the non-local length value^[87].

4.3 Transient grating experiments

The Transient Grating (TG) technique is a popular experimental method to investigate the onset of non-diffusive transport effects at the nanoscale. Here, two optical pulses interfere to generate a sinusoidal power density pattern in the sample, and the subsequent relaxation is measured to identify the transport mechanisms during the temperature field homogenization^[105,106]. This technique enables identifying a reduction of the apparent thermal conductivity at high temperatures by reducing the transient grating period in a variety of materials, including silicon^[106] and graphite^[60].

Within the hydrodynamic framework, this apparent reduction in conductivity is interpreted as a viscous effect arising from a strongly non-uniform heat flux profile. In situations where the grating period P approaches the ab initio value of the non-local length ℓ, the GKE predicts the emergence of viscosity and a highly impeded heat flow, in good agreement with experiments in silicon^[107]. However, similar to the experiments in nanostructures, the predictive power of the GKE with ab initio inputs is limited to relatively small ratios ℓ/P, preventing accurate predictions of thermal relaxation under nanometer scale optical excitations^[108]. Moreover, the viscous effect based on a characteristic non-local length scale, as predicted by the GKE, is inadequate to explain TG experiments in alloys, where the apparent thermal conductivity decreases over multiple orders of magnitude with decreasing grating period^[109,110].

Transient grating experiments have also been employed to investigate the other form of hydrodynamic behavior, namely memory effects or second sound. Specifically, an oscillating, wave-like evolution of the thermal field has been observed at temperatures ranging from 100 to 200 K in graphite under small grating period excitations^[60,111]. This phenomenon, and other signatures of second sound propagation in similar conditions^[61], have been directly predicted using numerical solutions of the BTE with the full collision operator^[60,81]. According to the GKE, second sound propagation is only possible in the absence of strong viscous effects^[83]. In this regime, where the viscous contribution is small compared to the diffusive and memory terms, the GKE solutions resemble those of the hyperbolic heat equation^[112,113], thereby capturing second sound propagation in the TG configuration^[52].

4.4 Time- and frequency-domain thermoreflectance

Time- and Frequency-Domain thermoreflectance experiments^[114] have been established as one of the main experimental techniques not only to measure intrinsic thermal properties of a wide range of materials, but also to uncover non-diffusive thermal transport mechanisms emerging at the nanoscale^[115,116]. Here, the sample is subjected to a time-modulated optical excitation, and the thermal response is monitored by measuring the optical reflectivity of the sample^[117]. The ability of the sample to dissipate heat from the transducer or surface region where the energy is injected affects both the amplitude of the thermoreflectance oscillations and their phase relative to the laser power, thus informing about the underlying transport mechanisms within the semiconductor. Once again, experimental measurements at small scales cannot be described using effective Fourier descriptions with a fitted conductivity or interfacial resistance^[53,116]. These non-diffusive signatures are a consequence of the same mechanisms inducing the slow relaxation of the nanoheaters discussed above. Indeed, the GKE with the same ab initio calculated inputs can describe non-Fourier signatures in both time-domain^[35] and frequency-domain^[53] thermoreflectance experiments in silicon and germanium over a wide range of temperatures.

Furthermore, this experimental configuration in the absence of metallic transducers can be used to investigate thermal evolution in response to highly confined energy density sources injecting heat directly into the semiconductor. Similar to nanoscale heat source relaxation experiments, in these conditions the non-local length can become geometrically limited by the penetration depth of the optical excitation. Consequently, the viscous term in the GKE is attenuated, which unlocks the propagation of thermal waves under high-frequency excitations^[33,118]. The emergence of second sound in transducerless FDTR experiments can be distinguished from ballistic transport effects because the non-diffusive behavior can be predicted using a single heat flux relaxation time-scale τ, rather than a multiscale phonon scattering spectrum. The observation of this single-scale feature in experiments unambiguously indicates that a fraction of the thermal energy can propagate as a wave even at high temperatures in standard semiconductors, such as germanium^[33].

4.5 Molecular dynamics

As evidenced by multiple experiments, the applicability of hydrodynamic modeling with ab initio parameter values is generally limited to length scales comparable or larger than the intrinsic non-local length ℓ, characteristic of a given material and temperature. In principle, this prevents interpretation of molecular dynamics simulations of standard semiconductors such as silicon, which typically consider very small nanostructures due to computational limitations. Nevertheless, a growing body of research shows that some non-equilibrium features observed in molecular dynamics simulations are reminiscent of hydrodynamic behavior, including non-local^[119-122] and memory^{[33,119,123,124]} effects. In particular, the heat flux profile established in structures down to a few tens of nanometers size displays a hydrodynamic-like shape that can be accurately fitted using the GKE with effective parameters^{[120,125,126]}, suggesting that the phonon population remains strongly correlated even at extremely high Knudsen numbers. This, in turn, indicates that current formulations of hydrodynamic theory with ab initio inputs could be refined and extended to extraordinarily small length scales, down to tens of nanometers.

5. Concluding Remarks and Outstanding Questions

Hydrodynamic models derived directly from the BTE, such as the GKE, represent a powerful framework to connect the ab initio characterization of the phonon transport properties with macroscopic, continuum-level descriptions. As demonstrated by a variety of experiments, this approach is useful to identify the breakdown of Fourier’s law of heat diffusion and predict the thermal response of devices at moderately small length and time scales. Nevertheless, the predictive capabilities of hydrodynamic models with ab initio inputs are generally limited to length or time scales comparable to or larger than the intrinsic value of the non-local length ℓ and the flux relaxation time τ. Although recent theoretical and experimental studies suggest that the applicability of this framework could be extended to smaller scales, a robust connection with ab initio calculations and molecular dynamics simulations is still missing at high Knudsen numbers. In particular, understanding the apparent reduction of the non-local length observed in experiments and molecular dynamics simulations at high Knudsen numbers^[87,120] stands as a key objective to push the applicability limits of the GKE with ab initio inputs to smaller scales.

From a phenomenological or experimental point of view, the most-defining feature of hydrodynamic phonon transport is the manifestation of single intrinsic phonon length and time scales, ℓ and τ, characterizing the interaction with boundaries and external heat sources, rather than uncorrelated phonon mode evolution associated with the multiscale phonon mean free path spectrum^[18]. This single-scale characteristic is compatible with a wide range of experimental observations in semiconductors such as graphite and silicon. In contrast, a variety of heat transport experiments in alloy semiconductors, such as SiGe, reveal a gradual onset of non-diffusive transport effects spanning multiple orders of magnitude in length scale that is not generally observed in the absence of alloy scattering. Therefore, elucidating and incorporating the role of alloy scattering within the phonon hydrodynamic framework remains as an important outstanding question.

Another crucial open question concerns the fundamental non-equilibrium conditions underlying the emergence of second sound. Experiments reporting this effect involve extremely confined heat sources within the material sample, which effectively attenuate viscosity and enable the manifestation of memory effects. A deeper understanding of the interplay between non-local and memory effects in the presence of extremely confined volumetric heat sources is therefore essential to establish the general conditions necessary for second sound propagation. To address this, further development of direct numerical solvers of the BTE in complex non-equilibrium situations represents a promising option to benchmark hydrodynamic modeling.

More generally, the efforts discussed in this review to identify phonon hydrodynamic behavior and to validate generalized hydrodynamic descriptions of phonon transport at the nanoscale parallel a growing body of work on electron hydrodynamics^[44,127-130]. Establishing clear analogies between these two domains, and uncovering possible coupled phonon-electron hydrodynamic regimes^[131,132], stands among the most promising directions for advancing our fundamental understanding of nanoscale transport phenomena.

In conclusion, the generalized hydrodynamic perspective presented in this review suggests that nanoscale thermal evolution can be interpreted using fluid-dynamics concepts such as viscosity and vorticity across a broad range of conditions beyond the restrictive applicability limitations that have been traditionally assumed. Since hydrodynamic modeling tools such as the GKE are formulated at the continuum level of description, they can be easily solved in complex conditions using finite elements^[43,113] or alternative numerical and analytical methods^[29,63]. This framework thus provides a robust theoretical basis for developing advanced thermal management strategies that leverage fluid-dynamics analogies, with potential applications to increase the efficiency and lifetime of energy-conversion and data-processing devices^[3,4,113].

Authors contribution

Beardo A: Conceptualization, writing-original draft, writing review & editing.

Tur-Prats J: Conceptualization, visualization, writing review & editing.

Camacho J, Alvarez FX: Conceptualization, writing review & editing.

Conflicts of interest

The authors declare no conflicts of interest.

Ethical approval

Not applicable.

Consent to participate

Not applicable.

Consent for publication

Not applicable.

Availability of data and materials

Not applicable.

Funding

This work is supported by the Spanish Ministerio de Ciencia, Innovacion y Universidades under Grant No.PID2021-122322NB-I00 and TED2021-129612B-C22 (MCIU/AEI/10.13039/501100011033/FEDER UE) and the AGAUR-Generalitat de Catalunya under Grant No.2021-SGR-00644.

Copyright

© The Author(s) 2026.