Generalization of Fourier's law into viscous heat equations

Heat conduction in dielectric crystals originates from the propagation of atomic vibrations, whose microscopic dynamics is well described by the linearized phonon Boltzmann transport equation. Recently, it was shown that thermal conductivity can be resolved exactly and in a closed form as a sum over relaxons, $\mathit{i.e.}$ collective phonon excitations that are the eigenvectors of Boltzmann equation's scattering matrix [Cepellotti and Marzari, PRX $\mathbf{6}$ (2016)]. Relaxons have a well-defined parity, and only odd relaxons contribute to the thermal conductivity. Here, we show that the complementary set of even relaxons determines another quantity --- the thermal viscosity --- that enters into the description of heat transport, and is especially relevant in the hydrodynamic regime, where dissipation of crystal momentum by Umklapp scattering phases out. We also show how the thermal conductivity and viscosity parametrize two novel viscous heat equations --- two coupled equations for the temperature and drift-velocity fields --- which represent the thermal counterpart of the Navier-Stokes equations of hydrodynamics in the linear, laminar regime. These viscous heat equations are derived from a coarse-graining of the linearized Boltzmann transport equation for phonons, and encompass both limits of Fourier's law and of second sound, taking place, respectively, in the regimes of strong or weak momentum dissipation. Last, we introduce the Fourier deviation number as a descriptor that captures the deviations from Fourier's law due to hydrodynamic effects. We showcase these findings in a test case of a complex-shaped device made of graphite, obtaining a remarkable agreement with the very recent experimental demonstration of hydrodynamic transport in this material. The present findings also suggest that hydrodynamic behavior can appear at room temperature in micrometer-sized diamond crystals.

Heat conduction in dielectric crystals originates from the propagation of atomic vibrations, whose microscopic dynamics is well described by the linearized phonon Boltzmann transport equation. Recently, it was shown that thermal conductivity can be resolved exactly and in a closed form as a sum over relaxons, i.e. collective phonon excitations that are the eigenvectors of Boltzmann equation's scattering matrix [Cepellotti and Marzari, Phys. Rev. X 6, 041013 (2016)]. Relaxons have a well-defined parity, and only odd relaxons contribute to the thermal conductivity. Here, we show that the complementary set of even relaxons determines another quantity -the thermal viscosity -that enters into the description of heat transport, and is especially relevant in the hydrodynamic regime, where dissipation of crystal momentum by Umklapp scattering phases out. We also show how the thermal conductivity and viscosity parametrize two novel viscous heat equations -two coupled equations for the temperature and drift-velocity fields -which represent the thermal counterpart of the Navier-Stokes equations of hydrodynamics in the linear, laminar regime. These viscous heat equations are derived from a coarse-graining of the linearized Boltzmann transport equation for phonons, and encompass both limits of Fourier's law and of second sound, taking place, respectively, in the regimes of strong or weak momentum dissipation. Last, we introduce the Fourier deviation number as a descriptor that captures the deviations from Fourier's law due to hydrodynamic effects. We showcase these findings in a test case of a complex-shaped device made of graphite, obtaining a remarkable agreement with the very recent experimental demonstration of hydrodynamic transport in this material. The present findings also suggest that hydrodynamic behavior can appear at room temperature in micrometer-sized diamond crystals. The present formulation rigorously generalizes Fourier's heat equation, extending the reach of physical and computational models for heat conduction also to the hydrodynamic regime.

I. INTRODUCTION
Thermal transport in insulating crystals takes place through the evolution and dynamics of the vibrations of atoms around their equilibrium positions. The first predictive theoretical framework to describe thermal transport was developed by Peierls in 1929 [1][2][3], who envisioned a microscopic theory in terms of a Boltzmann transport equation (BTE) for the propagation of vibrational excitations (phonon wavepackets). In the 1960s significant progress took place in this field, propelled by newly discovered hydrodynamic phenomena in crystals, with striking signatures such as Poiseuille-like heat flow [4] and second sound [5]. The former manifests itself with a heat flux that is akin to the flow of a fluid in a pipe (i.e. showing a parabolic-like profile with a maximum in the center and minimum at the boundaries, due to viscous effects); the latter instead results in heat propagation in the form of a coherent temperature wave, rather than a diffusing heat front. Second sound has been observed experimentally in a handful of solids; first, in solid helium [5], followed by sodium fluoride [6,7], bismuth [8], sapphire [9], and strontium titanate [10,11] -all at cryogenic conditions. Importantly, * michele.simoncelli@epfl.ch neither Poiseuille flow nor second sound can be described by the macroscopic Fourier's equation, which is limited to a diffusive description of heat propagation.
These experimental observations have been accompanied by several pioneering efforts aimed at providing a quantitative description of heat hydrodynamics at the mesoscopic level, i.e. in terms of partial differential equations (PDEs) that are simpler than the microscopic integro-differential BTE (we use "mesoscopic model" to denote any description requiring more fields or PDEs than the Fourier's PDE for the temperature field [12]). Sussmann and Thellung [13], starting from the linearized BTE (LBTE) in the absence of momentumdissipating (Umklapp) phonon-phonon scattering events, derived mesoscopic equations in terms of the temperature and of the phonon drift velocity, i.e. the thermal counterparts of pressure and fluid velocity in liquids. Further advances came from Gurzhi [14,15] and Guyer & Krumhansl [16,17] who, including the effect of weak crystal momentum dissipation, obtained equations for damped second sound and for Poiseuille heat flow. Among early works, we also mention the discussions of phonon hydrodynamics in the framework of many-body theory, as in Refs. [18,19]. While correctly capturing the qualitative features of phonon hydrodynamics, all the theoretical investigations mentioned above assume simplified phonon dispersion relations (either power-law [14,15] or linear isotropic [13,16,17]), or neglect momentum dissipation [13]. A more rigorous and general formulation -albeit valid only in the hydrodynamic regime of weak Umklapp scattering -was introduced by Hardy, who extended the study of second sound [20] and, together with Albers, of Poiseuille flow in terms of mesoscopic transport equations [21].
Further applications of the LBTE, combined with state-of-the-art first-principles simulations, have also recently predicted the existence of hydrodynamic phenomena at non-cryogenic temperatures in low-dimensional or layered materials such as graphene [31,39,40], other 2D materials [31], carbon nanotubes [41] and graphite [42]. Fittingly, and remarkably, these theoretical suggestions have now been confirmed by the experimental finding of second sound in graphite [43] at ∼ 100 K.
In the hydrodynamic regime, where Poiseuille flow or second sound occur, Fourier's law fails [43][44][45][46], depriving us of the most common tool used to predict the temperature profile in a device. The LBTE, in principle, allows to predict accurately thermal transport under these conditions, but its complexity prevents a straightforward application to materials with complex geometries (used in experiments and relevant for applications), thus posing limitations to the study of how a material's shape alters transport [35]. Recent research efforts have been directed at developing mesoscopic models that correct the shortcomings of Fourier's law and extend it at a lower computational cost than the solution of the full LBTE. Different strategies have been suggested; some approaches reduce the complexity of the LBTE by neglecting the effects of the phonon modes repopulation due to scattering events (the so-called single-mode relaxation-time approximation (SMA)), thus allowing for analytical [47][48][49][50] or asymptotic [51] solutions. From the LBTE in the SMA, mesoscopic models that generalize Fourier's law accounting for ultrafast thermal processes or ballistic effects have been derived [52][53][54][55][56]. Other works have derived mesoscopic models without relying on the LBTE [57], or have generalized the Guyer-Krumhansl equation to account for the effect of the boundaries on the heat flow [58][59][60][61][62]. A hydrodynamic transport model has been derived from the LBTE in the Callaway approximation, defining a phonon viscosity which can be computed from atomistic data [63].
Here, we provide a general and universal solution to the challenge of extending Fourier's law all the way to the hydrodynamic regime, deriving from the LBTE two novel coupled mesoscopic heat transport equations that cover exactly and on equal footing Fourier diffusion, hydrodynamic propagation, and all regimes in between.
To this aim, we first show that one can define the thermal viscosity of a crystal starting from an exact solution of the LBTE in terms of the eigenvectors of the scattering matrix (i.e. the relaxons introduced in Ref. [27] to determine thermal conductivity), and evaluate from it the crystal-momentum flux generated in response to a drift-velocity gradient. The relaxons' parity [27] highlights the complementary character of thermal conductivity and viscosity, with the former being determined by odd relaxons, and the latter by even relaxons.
Next, we use a coarse-graining procedure to derive two novel "viscous" heat equations: these are two coupled equations for the local temperature and drift-velocity fields, and are parametrized in terms of the thermal conductivity and viscosity. The viscous heat equations represent the thermal counterpart of the Navier-Stokes equations for fluids in the laminar regime, and, as mentioned, include Fourier's law and second sound in the limits of strong and weak crystal momentum dissipation, respectively.
Last, we introduce the Fourier deviation number (FDN), a dimensionless parameter that quantifies the deviation from Fourier's law due to hydrodynamic effects. We test this formalism on graphite, diamond and silicon, showing that the FDN predicts a temperature and size for the window of hydrodynamic thermal transport that replicates the explicit solution of the viscous heat equations, but at a negligible computational cost. Most importantly, the prediction of the present formulation for the hydrodynamic window of graphite is found to be in excellent agreement with recent, pioneering experiments [43]. In passing, we also predict that hydrodynamic behavior can appear in diamond at room temperature for micrometer-sized crystals.

II. THERMAL VISCOSITY
A microscopic description of thermal transport is given by the LBTE: where ν labels a phonon state (i.e. an index running on all the phonon wavevectors q and phonon branches s), v ν is the phonon group velocity, V is the crystal volume [64], and Ω νν is the phonon-phonon scattering matrix [27]. Eq. (1) governs the evolution of the deviation n ν (r, t) of the phonon populations from equilibrium: where N ν (r, t) are the out-of-equilibrium phonon populations at position r and time t,N ν = (e ων /(k BT ) − 1) −1 is the equilibrium Bose-Einstein distribution at temper-atureT , and ω ν are the phonon frequencies. From the solution of the LBTE one can derive the local lattice energy E(r, t)= 1 V ν ω ν N ν (r, t) and the total crystal momentum P (r, t)= 1 V ν qN ν (r, t) [20]. The former is often studied in connection with the thermal conductivity [27], while the latter becomes relevant in the hydrodynamic regime of thermal transport [13,15,65]. We note in passing that this latter emerges only in "simple" crystals, i.e. those where phonon interbranch spacings are much larger than their linewidths [38].
The energy flux generated in response to a temperature gradient determines the thermal conductivity; correspondingly, the crystal-momentum flux generated in response to a perturbation of the drift velocity determines the thermal viscosity (for the electronic analogous, see Ref. [66]). Therefore, we start by considering a crystal in the hydrodynamic regime of thermal transport (i.e. carrying a finite amount of crystal momentum); under the constraint of fixed total energy and crystal momentum, the local equilibrium is given by the phonon drifting The drifting distribution differs from the Bose-Einstein distribution due to the presence of a drift velocity u (a parameter expressing the amount of local momentum, just as the temperature does for the local energy); it depends implicitly on r, t through T (r, t) and u(r, t). Next, we study the effect of small perturbations in the temperature and drift velocity. To this aim, we expand the out-of-equilibrium drifting distribution (2) in proximity of the local thermal equilibrium [20], finding is referred to as the specific momentum, τ α is the relaxation time of relaxon α (i.e. the inverse eigenvalue associated to the eigenvector θ α ν of the symmetrized scattering matrixΩ νν [27]), and w j iα is the velocity tensor w j iα = 1 V ν φ i ν v j ν θ α ν for relaxon θ α ν and eigenvector φ i ν (see Appendix B for a full demonstration). φ i ν are three special eigenvectors linked to the crystal momentum of the system; to see this, we first decompose the scattering matrix asΩ νν =Ω N νν +Ω U νν , whereΩ N νν andΩ U νν contain only momentum-conserving (normal) and momentum-dissipating (Umklapp) processes, respectively. Since the normal part of the scattering matrix conserves crystal momentum, there exists a set of 3 eigenvectors φ i ν (i = 1, . . . , 3 where 3 is the dimensionality of the system) with zero eigenvalue forΩ N νν , which are associated to the conservation of crystal momentum in the 3 Cartesian directions. Because the viscosity describes the response of the crystal momentum flux to a change of drift velocity, it is not surprising that the eigenvectors φ i ν appear in its definition. In fact, the local equilibrium distribution (Eq. (3)) is linear in the drift velocity, with the proportionality coefficients being these special eigenvectors (see Appendix A for a proof), and therefore they appear in the viscosity as well, to describe a perturbation to this local equilibrium. We note in passing that the thermal viscosity defined in Eq. (9) has the dimensions of a dynamic viscosity (Pa·s) and satisfies the property µ ijkl = µ ilkj = µ kjil .
We report in Fig. 1 the first-principles estimates of the thermal viscosity for graphite, diamond, and silicon. We choose these three materials as prototypes for two different behaviors, the former two displaying hydrodynamic thermal transport [25,42,43,45,69,70], as opposed to the more conventional case of silicon [25,27]. We account for finite-size effects by combining the bulk viscosity in Eq. (9) with its ballistic limit via Matthiessen's rule for a sample having size (or grain size) of 10 µm; the same renormalization is applied also to the thermal conductivity (see Appendix C for details). This is an approximate treatment of surface scattering effects, which provides an estimate of size effects at a much lower computational complexity compared to more refined models [69,71]. In general, when samples' sizes become comparable or smaller than the carriers' diffusion lengths [61,[72][73][74][75][76], transport coefficients cannot be rigorously defined and a more accurate treatment is required, solving the much more complex and computationally expensive space-dependent LBTE [24,35,47,77]. Matthiessen's approach is a good approximation when the values of thermal conductivity and viscosity do not differ significantly from their bulk counterparts, as is the , and silicon (c), as a function of temperature (the off-plane tensor components for graphite are reported in Fig. 9). Insets: total thermal conductivities (in-plane components for graphite and diagonal components of the isotropic tensor for diamond and silicon) as a function of temperature. The dashed lines refer to the bulk materials (Eq. (9) for viscosity and Eq. (B12) for conductivity); the solid lines show finite-size predictions computed using Eq. (C1) and Eq. (C2) for a sample having size (or grain size) 10 µm (see main text and Appendix C; for diamond and silicon the component µ ijji is negligible).
case for the materials, sizes and temperatures considered in this work (e.g., in Fig. 1 the size renormalization is negligibly small in the temperature range considered for graphite and silicon, and becomes so above 300 K for diamond).
For graphite in the low-temperature ( < ∼ 100 K) regime, the bulk in-plane viscosity components are constant or slowly increasing with temperature; at higher temperatures, the viscosity components increase with temperature up to reaching their saturation values. For a graphite polycrystalline sample with a grain size of 10µm, finite-size effects on viscosity (and on conductivity) are small. For diamond in the temperature regime below 90 K, the bulk values of the thermal viscosity increase with increasing temperature, then in the intermediate temperature regime (100 K < ∼ T < ∼ 300 K) they decrease with temperature, and finally for increasingly higher temperatures they gently increase up to saturating to the high-temperature limiting values. However, size effects renormalize this behavior for temperatures below 300 K, leaving viscosity components that increase with temperature. In silicon, increasing temperature yields a decrease, up to a high-temperature saturation value, of the bulk thermal viscosity; size effect renormalize the viscosity below ∼ 100 K yielding viscosity components that increase with temperature (i.e. a behavior analogous to that of diamond). We note that in the high-temperature limit, i.e. when the thermal conductivity decay as T −1 [3,78], the viscosity components tend to constant values for all these three materials. The total bulk thermal conductivities are shown in the insets of Fig. 1 for comparison (we only show one component; for diamond and silicon the conductivity tensor is isotropic, for graphite we show the in-plane component). The off-plane transport coefficients for graphite are discussed in Fig. 9a. In graphite, the bulk thermal conductivity below 100 K does not increase monotonically as the temperature is decreased due to the presence of natural-abundance isotopic scattering. Natural-abundance isotopic scattering has been considered also for diamond and silicon, and the thermal conductivity at low temperatures is much lower than for isotopically-pure diamond [25,[79][80][81][82][83][84] or silicon [85]. We have also verified that the effects of phonon coherences are negligible in these crystals [38]. We note in passing that, even if the thermal conductivities of diamond and silicon differ by more than one order of magnitude, their largest thermal-viscosity components differ only by a factor of 3. These results may be compared with the case of water, whose dynamic (shear) viscosity is 8.5·10 −3 Pa·s at room temperature, indicating that the thermal viscosity found here is comparable or larger. To a good approximation, water is an incompressible fluid, and thus its largest viscosity component is ijij, also called "first viscosity" or "shear viscosity". For compressible fluids, the iiii components of the viscosity tensor -also called "second viscosity" or "volume viscosity" [86] -are nonnegligible. Here, in contrast with water, the iiii component of the thermal viscosity tensor is the largest, which Relaxons' contributions to the bulk thermal conductivity and to the thermal viscosity at 300 K for graphite (a, in-plane components), diamond (b) and silicon (c). Each dot represents a relaxon, with its color labeling its relaxation time, and its area being proportional to the sum of its percentage contributions to the thermal conductivity and viscosity. The dashed lines are plotted as a guide to the eye, to underscore how even and odd relaxons are fully decoupled. Odd relaxons, which determine the thermal conductivity, yield negligible (zero) contributions to the thermal viscosity; conversely, even relaxons determine the thermal viscosity and yield negligible (zero) contributions to thermal conductivity.
elicits an analogy to a compressible fluid. It is important to mention that the present formulation may need further extension to deal with unstrained 2D materials, since their flexural phonon modes have a quadratic dispersion ω ν ∝ q 2 at T = 0K [87]. Under these conditions, long-wavelength phonons lead to q · u being larger than ω ν , causing negative values of the phonon drifting distribution which are not compatible with a semiclassical description of transport. Further work is needed to address this issue, for example considering the phonon renormalization due to coupling between bending and stretching degrees of freedom of the monolayer [88] or the presence of a substrate [89], or introducing the Wigner-function formalism [38,90].
Finally, it is worth checking the complementary character of the thermal conductivity and viscosity that arises from their decoupled relaxons' contributions. As commented above, decomposing the thermal conductivity [27] and viscosity (9) in terms of single relaxons it is possible to show that the thermal viscosity is uniquely determined by the even part of the relaxon spectrum, while the thermal conductivity is determined uniquely by the odd part of the relaxon spectrum [27]. In Fig. 2 we highlight the contributions of each relaxon to the total thermal conductivity and viscosity, confirming numerically this picture.

III. VISCOUS HEAT EQUATIONS
We show here that heat conduction can be described by two novel viscous heat equations that cover on the same footing both the Fourier and hydrodynamic limits, and all intermediate regimes. These are two coupled equations in the temperature T (r, t) and drift-velocity u(r, t) fields, which are parametrized by the thermal viscosity and conductivity. These equations represent the thermal counterpart of the Stokes equations of fluid dy-namics -i.e. the Navier-Stokes equations in the linear regime, whose solution yields the laminar flow -where temperature takes the role of pressure and the phonon drift velocity that of the fluid velocity. In the kinetic regime, when momentum-dissipating (Umklapp) scattering processes dominate [31], these viscous heat equations become equivalent to Fourier's heat equation.
As underscored before, hydrodynamic thermal transport is characterized by energy conservation and crystal momentum quasi-conservation (the latter being exactly conserved only in absence of Umklapp processes [13]). Conserved quantities in the LBTE dynamics can be related to the eigenvectors of the full or normal scattering matrix with zero eigenvalues [20] (see also Appendix A). Four of these eigenvectors (i.e. phonon distribution functions) can be identified. The first one is the Bose-Einstein eigenvector φ 0 ν ∝ ω ν ∝ ∂Nν ∂T ∝ñ T ν , which is the eigenvector of zero eigenvalue for the symmetrized full scattering matrixΩ νν ; its zero eigenvalue is associated to energy conservation in scattering events (both normal and Umklapp). The other three eigenvectors are the drift eigen- is the dimensionality of the system) already introduced for the evaluation of the viscosity; these φ i ν are eigenvectors with zero eigenvalue for the normal part of the scattering ma-trixΩ N νν and are associated to the conservation of crystal momentum by normal scattering events [13,16,20]. We note that these four eigenvectors constitute the first two terms of the phonon distribution expansion in Eq. (4). We can thus derive the mesoscopic equations that govern the evolution of the temperature (T (r, t)) and drift velocity (u(r, t)) fields projecting the microscopic LBTE in the subspaces spanned by φ 0 ν and by φ i ν (i = 1, . . . , 3). In order to derive a closed-form equation for the drift velocity, when projecting in the subspace spanned by φ i ν we consider the effects of momentum dissipation only within that subspace. The result is the following set of equations (see Appendix D for a detailed derivation): ν is a velocity tensor that arises from the non-diagonal form of the diffusion operator in the basis of the eigenvectors of the normal part of the scattering matrix (see Appendix D),T is the reference (equilibrium) temperature on which a perturbation is ap-plied (see Sec. V and Appendix A 1), A i is the specific momentum in direction i (defined in Sec. II), µ ijkl is the thermal viscosity tensor, κ ij is the thermal conductivity tensor [27], and D ij U = 1 V νν φ i νΩ U νν φ j ν is the momentum dissipation rate; the latter is caused both by the presence of Umklapp processes as well as boundary scattering (i.e. D ij U is sensitive to size effects like the thermal conductivity and viscosity). Treating boundary scattering as in Ref. [69], we compute D ij U as D ij νν φ j ν is the momentum dissipation caused by Umklapp scattering in the bulk and D ij U,boundary ( is the momentum dissipation caused by scattering with boundaries at a distance L S . The scalar equation (10) and the vectorial (3-components) equation (11) rule the coupled evolution of the scalar temperature field and of the vector drift-velocity field; they constitute the main result of this work and we name these "viscous heat equations". These transport equations are reminiscent of the linearized Stokes equations for fluids: to see this more clearly, we note that local energy E(r, t) and crystal momentum P (r, t) are proportional to temperature and drift velocity respectively (E(r, t) = C T (r, t) and where C is the specific heat and A is the specific momentum). Exploiting these relationships, it is possible to rewrite the viscous heat equations in a more familiar form, namely as energy E and momentum P i balance equations: where, on the basis of the phonon population expansion in Eq. (4), we distinguish the drifting heat flux into the contributions from the local temperature gradient Q δ,i (r, t) = − j κ ij ∇ j T (r, t) and from the drift veloc- Similarly, the momentum flux receives separate contributions from the local temperature Π ij T (r, t) = CA i /T W j i0 T (r, t) and from variations of the drift velocity through Eq. (8). Finally, ∂P i ∂t Umkl accounts for the dissipation of crystal momentum by Umklapp processes or scattering with boundaries; further details can be found in Appendix D. The distinction between temperature-driven and drifting components of the heat flux (Q δ (r, t) and Q D (r, t), respectively) is essential for hydrodynamic transport. We will show in Sec. VI that Fourier's law is recovered in the limiting case where crystal momentum dissipation dominates over viscous effects and is the fastest timescale of the phonon dynamics.
It is worth mentioning that the viscous heat equations introduced here differ from the Stokes equations for fluids in two major ways. First, there is no analogous to the mass conservation satisfied by Stokes equations, since the total phonon number is not a constant of motion (e.g. a phonon coalescence event decreases the number of phonons in the system). Second, while collisions between molecules in the fluid conserve momentum, scattering among phonons does not conserve crystal momentum in presence of Umklapp processes.
The most relevant feature of the viscous heat equations is their capability to describe hydrodynamic thermal transport in terms of mesoscopic quantities, i.e. temperature and drift velocity, resulting in a much simpler and computationally less expensive approach than the microscopic LBTE. The parameters entering Eqs. (10,11) can be determined from first-principles calculations or, possibly less accurately, from classical potentials, and are tabulated in Appendix E for graphite, diamond and silicon.
In order to be solved, the viscous heat equations require appropriate boundary conditions on the temperature and drift velocity. Boundary conditions on temperature have been widely studied in conjunction with Fourier's heat equation [91]: typically, one makes assumptions on the system's capability to exchange heat at the boundaries, and on the temperature at those boundaries (Neumann and Dirichlet boundary conditions, respectively [91]). In the next section V we consider a system in which the temperature is fixed at some boundaries, while the others are assumed to be adiabatic (that is, the heat flux across these boundaries is zero). In contrast, boundary conditions on the drift velocity, i.e. on crystal momentum at the sample's borders, have not been studied as extensively. Since crystal momentum is not conserved at boundaries [13], we impose a no-slip condition of zero drift velocity u(r, t) on all boundaries, ensuring thus zero drifting heat Q D (r, t) ∝ u(r, t). As discussed in past work [3,35], more comprehensive boundary conditions require quantifying phonon reflection at surfaces, and are beyond the scope of this work.
The viscous heat equations (10,11) improve upon past work on different levels. First, they are valid for general dispersion relations for phonons, and thus overcome the limitations -pointed out by Hardy and Albers [21] already in 1974 -of the pioneering mesoscopic models developed in the 1960s by Guyer-Krumhansl [16,17] or Gurzhi [15], which assumed linear-isotropic or powerlaw phonon dispersion relations, respectively (assumptions which arose from the hypothesis that hydrodynamic phenomena would occur at cryogenic temperatures only). The limitations of the Guyer-Krumhansl [16,17] and Gurzhi [15] models were overcome in 1974 by Hardy and Albers [21], who derived from the LBTE a set of mesoscopic equations for energy and crystal momentum valid for a general phonon dispersion relation (i.e. not necessarily linear-isotropic). Hardy and Albers' formulation relies on the hydrodynamic approximation -i.e. that the fastest timescale of phonon dynamics is that of normal processes -that is valid only within the hydrodynamic regime, where Umklapp collisions are rare events. In fact, Hardy and Albers' equations are limited just to the hydrodynamic regime and do not incorporate Fourier's law as a limit (and thus any intermediate regimes). The viscous heat equations address exactly this issue, and Fourier's law emerges when crystal momentum dissipation (due to Umklapp processes or scattering with boundaries) dominates over viscous effects and is the fastest timescale of the phonon dynamics (see Sec. VI). Finally, the viscous heat equations take into account the entire collision matrix, at variance with recent mesoscopic models derived from the LBTE either in the SMA [60] or in the Callaway approximation [63].

IV. SECOND SOUND
Second sound is the coherent propagation of a temperature wave [13,17,19,20,40,43,[92][93][94][95], and it is an effect properly described by the viscous heat equations. From a phenomenological point of view, second sound appears when the temperature field satisfies the following damped wave equation [20] (we define x as the second sound propagation direction): where τ ss and v ss are the second-sound relaxation time and undamped propagation velocity, yet to be determined. In contrast to Fourier's law, which has the form of a classical (non-relativistic) diffusion equation and states that a temperature-gradient variation causes an instantaneous variation of the heat flux, Eq. (14) has the physical property that a sudden localized change of temperature propagates in space with a finite speed (i.e. it is not felt instantly everywhere in space) [12,[96][97][98]. The temperature profile that solves Eq. (14) has the form of a damped wave: where the second sound frequencyω(k) depends on the second-sound wavevector k. In Appendix F we derive the second-sound equation from the viscous heat equations (10,11) following two different approaches (bottomup and top-down). In the first bottom-up approach we find the conditions for which the damped wave equation (14) emerges form the viscous heat equations Eq. (10,11). When this happens, the solution of Eq. (14) is the damped wave equation for temperature (15) shown above, with the second-sound dispersion relation given byω(k) = v 2 ss k 2 − (2τ ss ) −2 (this can be easily verified substituting Eq. (15) into Eq. (14)). This allows to express τ ss and v ss in terms of the parameters appearing in the viscous heat equations; in particular, . The propagation velocity of second sound is affected by damping and depends on the wavevector k: it is given by the group velocity v g (k) = ∂ω(k) ∂k = kv ss [k 2 + (2τ ss v ss ) −2 ] − 1 2 , and we note that it reduces to the undamped propagation velocity v ss in the undamped limit τ ss → 0.
These results are consistent with empirical expectations on second sound: in the limit of weak crystal momentum dissipation the second-sound relaxation time increases, while the velocity becomes smaller, making second sound more likely to be observed in the hydrodynamic regime [17,20]. In fact, when D xx U → 0 we find that τ ss → (D xx U ) −1 and v ss → W x x0 . We note that the viscous heat equations describe not only the propagation of the temperature field, but also that of the drift velocity. In Appendix F we show that when second sound emerges the drift-velocity field propagates as a damped wave as well (i.e. similar to Eq. (15)), with the same relaxation time and velocity of temperature, but with a phase shift of π/2.
Phase shifts between hydrodynamic and resistive components of the frequency-dependent heat flux have also been discussed recently in the context of an improved Callaway approximation for the LBTE [99].
As an alternative, we took inspiration from Ref. [40], which derived the second-sound dispersion relations by taking advantage of the Laplace transform of the LBTE, to identify solutions in the form of a damped wave following a top-down approach. In particular, we take a damped-wave ansatz for temperature, Eq. (15), and a similar one for the drift-velocity field (with the same frequency and decay time of temperature); we substitute these into the viscous heat equations (10,11) and determine the conditions under which these are acceptable solutions. As detailed in Appendix F, we find that the dispersion relations of second sound in the long-wavelength limit reduce to ω In this long wavelength limit and in the hydrodynamic regime and v g (k) ≈ W , which is consistent with the first bottom-up approach presented to derive the equations for second sound. We further stress that the LBTE can only rigorously describe second sound in the long-wavelength limit k → 0 in order for the temperature to be slowly varying (for which the two approaches shown in this section provide the same result); for smaller wavelengths, the definition of temperature itself becomes questionable [100].
We also recall that Enz [19] and Hardy [20] distinguished between "drifting" and "driftless" second sound. The former emerges when normal processes dominate and is described in terms of balance equations for energy and momentum; the latter is determined by a uniform energy flux that decays exponentially. The second sound discussed here is of the drifting kind, as it emerges from a set of balance equations for energy and crystal momentum derived from the LBTE.

V. CASE STUDY
We showcase the solution for the viscous heat equations for graphite around the equilibrium temperaturē T = 70 K in the geometry shown in Fig. 3, often used as an illustrative example in textbooks on fluid dynamics. The equations are solved numerically using a finiteelement solver implemented in Mathematica [101], imposing a temperature of 80 K on the left side (x = 0 µm) and of 60 K on the right side (x = 15 µm), assuming all ∂r i ∂r j = 0. Boundary conditions are as follows: the temperature is set at 80K/60K on the left/right boundaries, a zero total heat flux is imposed at the other boundaries, and zero drift velocity is imposed on all boundaries (no-slip boundary condition). With this choice of boundary conditions, the phonon distribution reduces to the Bose-Einstein distribution on the left/right boundaries, where the temperature is fixed. Despite having local equilibrium at the left/right boundaries, a non-zero drift velocity is present in the material as a consequence of the coupling in Eq. (10) between drift velocity and temperature. Panel c) shows the x-component of the heat flux along sections taken at x = 1.5 µm (orange, gray, red) and x = 9.0 µm (black, blue, green); solid lines are the results obtained from the viscous heat equations, dashed lines are the results from Fourier's law, and dotted lines are the heat fluxes due to the temperature gradient within the viscous heat equations (Q δ ).
boundaries at x = 0 and x = 15 µm to be adiabatic, and imposing a no-slip condition on u at all boundaries.
We show in Fig. 3 the viscous heat flux Q δ + Q D , and the flux predicted by Fourier's law Q Fourier (left panels) and contrast the Fourier and viscous heat flux components across two different sections (right panel).
We stress that Fourier's law lacks a description of the contribution to the heat flux derived from the local drift velocity [13,102,103]. As a result, Fourier's law misses qualitative and quantitative properties of the heat flux profile. The largest differences are observed in proximity of spatial inhomogeneities, such as boundaries or corners. For example, Q D quickly increases (decreases) in proximity of the thermal reservoir on the hot-left (cold-right) side of the sample; these changes in Q D determine opposite changes in Q δ . Microscopically, these variations are caused by the rapid transition of the phonon distribution from the Bose-Einstein equilibrium distribution imposed at the boundaries to an out-of-equilibrium distribution carrying non-zero total crystal momentum (i.e. a drift velocity) inside the sample.
We report in Fig. 3c the total heat flux profiles along two transversal sections of the sample, contrasting the prediction from the viscous heat equations (solid lines) with that of Fourier's law (dashed lines). Along these sections, Fourier's law predicts a flat heat flux profile, while the viscous heat equations yield a Poiseuille-like profile. The results from the viscous heat equations are thus substantially different from Fourier's predictions, and the behavior of the heat flux can be understood from a simple analytical 1D solution of the viscous heat equations in the absence of Umklapp processes [13]: as discussed in Appendix G, the flux is described by hyperbolic functions with a characteristic length scale b = µκ AC(W ) 2 (an estimate of the friction lengths, see Ref. [35]). At distances from the surface larger than b we recover the flat behavior typical of the bulk. We also note that these results mimic those from the (computationally very ex- pensive) space-dependent solution of the LBTE (either in the frequency-dependent SMA approximation [77] or considering the full scattering operator [35]), which generates a minimum flux on surfaces and maximum at the sample's center. We further note that, at variance with classical fluid dynamics and as pointed out in Ref. [35], the total heat flux does not drop to zero at the boundaries: the no-slip condition sets the drifting heat flux Q D to zero, but the temperature-driven component Q δ is still allowed to be nonzero.
In Fig. 4 we plot the difference between the temperature profile predicted by Fourier's law and the viscous heat equations along longitudinal (top panel) and transversal (bottom panel) directions. The insets in Fig. 4b show the results of the viscous equations (solid blue line) and Fourier's law (dashed red line) along the section y = 4.5 µm. Along the transversal direction (Fig. 4a), Fourier's law and Eqs. (10,11) predict temperature profiles which are substantially different in the presence of variations of the sample's shape (green line corresponding to x = 3.3 µm), while they are merely shifted by a positive or negative offset away from these; the precise amount depends on the distance from the fixed-temperature boundaries. These differences become more clear by inspecting the longitudinal direction (Fig. 4b), where the discrepancy between the temperature predicted by Fourier's law and Eqs. (10,11) is largest at x 3 µm, i.e. where the sample of Fig. 3 changes geometry. We show in the inset of Fig. 4b that also the longitudinal temperature profile (for y = 4.5 µm) changes when going from Fourier's law (dashed red line) to the viscous heat equations (solid blue line). The difference is maximized close to variations of the sample's shape at x = 3 µm and at the boundaries x = 0 and x = 15 µm; in this latter case, the viscous heat equations predict steeper-than-Fourier's law or non-linear temperature gradients that are reminiscent of those obtained in molecular dynamics simulations [104][105][106][107][108][109][110][111][112][113] and in explicit solutions of the LBTE [35].

VI. FOURIER DEVIATION NUMBER
In this section, we introduce a descriptor that parametrizes the conditions under which hydrodynamic heat conduction is observed; we will refer to this as the "Fourier deviation number" (FDN). In particular, we aim at distinguishing the diffusive regime from the hydrodynamic regime: in the former case the viscous heat equations become equivalent to Fourier's law, while in the latter case Fourier's law no longer holds and the viscous heat equations are required.
In order to investigate how the hydrodynamic deviations from Fourier's law depend on the sample's sizes and reference temperatureT , we solve the viscous heat equations and Fourier's equation for several samples having different dimensions, and at various reference tempera-turesT . We perform this analysis for graphite (Fig. 5a), diamond (Fig. 5b) and silicon (Fig. 5c). In particular, we consider samples that are similar -in the geometric sense -to the reference sample shown in Fig. 3, i.e. we generate samples of different dimensions starting from the reference sample of Fig. 3 and varying the sizes via a uniform scaling. Following this protocol, we vary the sample's length l TOT , reported on the y−axis of Fig. 5ac, from 0.1 to 100 µm (e.g. the reference sample of Fig. 3 corresponds to y = 15 µm, y = 100 µm corresponds to a sample obtained magnifying uniformly by a factor of 100/15 6.66 the sizes of the reference sample). The color in Fig. 5a-c represents the normalized L 2 distance between the temperature profile predicted by Fourier's law and the viscous heat equations for a given sample length l TOT and reference temperatureT : and in order to ease the qualitative interpretation of these results later, we evaluate L 2 in the spatiallyhomogeneous region G defined by x > 1 5 l TOT . We also inspected the effects of shape on the deviation L 2 , finding that the magnitude of L 2 is larger for geometries with spatial non-homogeneities that imply larger values of the drift velocity's second derivative; nevertheless, the dependence of L 2 on l TOT andT is qualitatively unchanged for different shapes. Results in Figs. 5a-c have been computed accounting for finite-size effects as discussed in Sec. II and Appendix C. Additionally, to better compare with the experimental results of Ref. [43], we make the hypothesis of working with polycrystalline graphite with an average crystal grain size of 10 µm as in Ref. [43]: as a result, the transport coefficients of graphite are renormalized by the system size for l TOT < 10 µm, and by the grain size for l TOT > 10 µm (therefore, for l TOT > 10 µm size effects are given only by the boundary conditions effects on the temperature and drift-velocity fields). In graphite, the difference between Fourier's law and the viscous heat equations is largest in the temperature range 60-80 K and for sizes 10µm < ∼ l TOT < ∼ 20µm. Turning our attention to diamond, we first recall that its thermal properties are characterized by a very large thermal conductivity, which originates from having large group velocities and weak Umklapp scattering [79][80][81][82][83][84]; the latter condition being favorable to the emergence of hydrodynamic effects. This is confirmed by the results shown in Fig. 5b, with the hydrodynamic deviation being largest around room temperature and dimension l TOT > ∼ 1 µm (we note in passing that, at this temperature and size, the surface-renormalization of transport coefficients is small and well compatible with the mesoscopic approach described in this work). These results suggest that a hydrodynamic window exists also for diamond, and that therefore hydrodynamic behavior (e.g. second sound) might be measurable in this material at temperatures even larger than graphite.
For silicon, instead, the deviation from Fourier's law is smaller and takes place at lower temperatures. Silicon is an example of a material for which the thermal conductivity computed within the SMA approximation is very close to the conductivity computed from the exact solution of the LBTE [23,25,27,114], and it is known that the SMA works reasonably well in systems where momentum-dissipating (Umklapp) scattering events dominate over those conserving crystal momentum (normal) [34,114]. Therefore, the negligible magnitude of hydrodynamic effects predicted by the viscous heat equations for silicon is consistent with the predominance of Umklapp over normal scattering events known to occur in this material. We last remark that, in contrast to graphite, results for diamond and silicon have been computed for single crystals, i.e. transport coefficients are not limited by grains' size.
In order to capture intuitively and inexpensively all these trends we introduce an approach inspired by the definition of the Reynolds number, and we rewrite the viscous heat equations (10,11) in adimensional form (we follow the standard procedure used e.g. in fluid dynamics, which is also called "Buckingham Pi theorem" [115,116]). First, we extract the magnitudes of the tensors appearing in the viscous heat equations, factorizing the largest component: we define a set of dimensionless variables r * = r/L, u * = u/u 0 , and T * = T /δT , where L, u 0 and δT are a characteristic size, drift velocity and temperature perturbation (more on these later). Substituting these variables in Eqs. (10,11), and limiting ourselves to the steady-state regime, we obtain: where we have introduced the quantities The role of these parameters is to account for the correct order of magnitude of the five summations appearing in Eqs. (17,18); they carry a Cartesian superscript i whenever they may depend on direction, and their subscripts indicate the number of indexes summed in their definition. For example, for a 3×3 isotropic tensor χ ij = χ m δ ij (where χ m is a constant and δ ij is the Kronecker tensor), g(χ) 2 = (χ m ) −1 ij χ ij = 3 and g i (χ) 1 = (χ m ) −1 j χ ij = 1 ∀ i. Because we are computing the FDN for the planar geometry discussed in Fig. 3, we can discard the z component (i.e. the Cartesian direction indexed by 3) from Eqs. (17,18). Then, for diamond, silicon and in-plane graphite, transport properties are isotropic (i.e. all the 2 nd -rank tensors appearing in Eq. (17) and Eq. (18) are diagonal). From this it follows that g(w· √ a) 2 = 2; g(k) 2 = 2; g i (w) 1 = 1 ∀ i = 1, 2; and g i ( √ a⊗a·d) 1 = 1 ∀ i = 1, 2. Finally, the tensor m ijkl is not isotropic  (10,11), shown as a function of temperature and total length of a sample similar (in the geometric sense) to the shape of Fig. 3 for the cases of graphite (a), diamond (b) and silicon (c). The color quantifies the L 2 distance between the temperature profiles predicted by Fourier's law (TFourier(x, y)) and the viscous heat equations (10, 11) (TViscous(x, y)), computed over the region G corresponding to points with x > 1 5 lTOT where lTOT is the total length of the sample geometrically similar to that in Fig. 3 (see main text for details). Bottom row, Fourier Deviation Number (FDN, as defined in Eq. (22)) for the same materials as a function of temperature and total length of the sample (for the geometry of Fig. 3, the characteristic size L appearing in FDN is the shortest length-scale, i.e. 1/5 of the total length lTOT). It is apparent that the FDN correctly identifies the deviations between the solutions of the Fourier and of the viscous heat equations solutions.
The final expressions for the dimensionless parameters π 1 , π 2 , π 3 are thus: From the derivation of these parameters, it is clear that they can be interpreted in terms of average values (we indicate with . . . the average over space) of physical quantities: π 1 ∼ Q D / Q δ , π 2 ∼ ∂Π T ∂r / ∂P ∂t Umkl and π 3 ∼ ∂Π δE ∂r / ∂P ∂t Umkl . Still, to evaluate these parameters we need to know the characteristic size L and temperature perturbation δT , and estimate the characteristic values of the drift velocity u 0 . Focusing on the setup discussed in the previous section (Fig. 3), we have L = 1 5 l TOT (corresponding to the shortest length scale of the geometry considered, which has total length l TOT along x) and δT =10K. As shown in Appendix H, the characteristic drift velocity u 0 is found by interpolating the asymptotic behavior at low (u L ) and high (u H ) temperatures u −1 0 = u −1 L +u −1 H . In the low-temperature limit, where Umklapp scattering is frozen, viscous effects determine the drift velocity, and one can show that u L = . In the high temperature limit, the drift velocity is mainly determined by the crystal momentum dissipation rate and u H = δT L . We can therefore estimate the values of all the π 1 , π 2 , and π 3 factors.
Looking at the definition of π 1 , π 2 and π 3 , it is straightforward to identify the conditions for which the viscous heat equations reduce to Fourier's law. When π 3 1, i.e.
for crystal momentum dissipation dominating over viscous effects and/or very large characteristic size L, the temperature gradient in Eq. (18) is proportional to the drift velocity: . While it is intuitive to understand the emergence of Fourier's law when crystal momentum dissipation (e.g. due to Umklapp scattering) dominates over viscous effects, the Fourier-like behavior of large-size samples can be rationalized qualitatively recalling the simple analytical 1D solution of the viscous heat equations in the absence of Umklapp processes discussed in Appendix G: the heat flux is described by hyperbolic functions with a characteristic length scale b = µmκm AmC(Wm) 2 , and for a characteristic size L b, a flat Fourier-like heat flux profile is recovered. It is worth mentioning that π 3 1 -and thus Fourier-like behavior -can also emerge for very small (sub-micrometer) length scales, where ballistic effects are relevant and strongly renormalize the transport coefficients, as discussed in Appendix C. The emergence of Fourier-like behavior in the ballistic regime (with a size-dependent thermal conductivity alike to that employed here) has been explained in a recent work [117], confirming that the results for sub-micrometer sizes obtained with the aforementioned approximated treatment of boundary scattering are qualitatively correct. If we insert the condition π 3 1 in Eq. (17), we obtain a Fourier-like equation for the temperature: ∂r * i ∂r * j =0. So, the thermal conductivity is corrected by a term proportional to π 1 π 2 (since in all the setups considered here, g(k)2g i ( a)2g i (w)1 = 1), but the qualitative behavior is that of Fourier's law. Numerically, we find that π 1 π 2 1 for graphite in the high-temperature limit (see Fig. 6a), so that the viscous heat flux is well approximated by Fourier's flux; the same limiting behavior is observed also in diamond and silicon. On the other hand, deviations from the temperature profile predicted by Fourier's law appear when both π 1 and π 3 are large. In fact, π 3 > ∼ 1 implies that the temperature gradient in Eq. (18) is coupled to the second derivative of the drift velocity, and one can not obtain a Fourier-like equation as before; large values of π 1 imply that the drift velocity affects the evolution of temperature in Eq. (17). It follows that Fourier's law is mathematically recovered for crystal momentum dissipation dominating over viscous effects and/or large characteristic size (π 3 1), but deviations between the viscous and Fourier temperature profile can be small even when viscous effects dominate over crystal momentum dissipation (π 3 1) if the coupling between drift velocity and temperature in Eq. (17) is small (π 1 1). We note in passing that while this reasoning has been made here at steady state, it can be straightforwardly extended to the time-dependent case, considering time variations of the drifting velocity slow compared to the crystal momentum dissipation timescale.
On the basis of this reasoning, it is convenient to in-troduce a Fourier deviation number (FDN) as: which provides a simple estimate of the deviations from Fourier's law: the larger the FDN, the larger the deviations from Fourier's law. In Fig. 5d-f, we plot the FDN for graphite, diamond and silicon as a function of sample's total length l TOT and reference temperaturē T . Remarkably, the FDN captures the trends of the exact solution of the viscous heat equations (Fig. 5a-c), thus identifying accurately the regime where hydrodynamic behavior emerges. The magnitude of FDN is directly proportional to the magnitude of hydrodynamic effects: a larger FDN implies a larger difference between  [43]. The terms π2 and π1π2 are also reported in panel a) because in the high-temperature limit Fourier's law is recovered when π1π2 1 (see main text for details). Notably, FDN predicts hydrodynamic effects to be largest at temperatures and sizes that are in excellent agreement with these observed in the experiments.
the Fourier and viscous temperature profiles.
A detailed analysis of the π terms as a function of temperature for a graphite sample geometrically similar to that in Fig. 3 and with a fixed total length l TOT = 10µm is shown in Fig. 6a, where l TOT = 10µm has been chosen to match the length fixed in the experiments by Huberman et al. [43] to investigate the magnitude of hydrodynamic effects as a function of temperature. Importantly, the FDN at fixed l TOT = 10µm predicts hydrodynamic effects to be largest at temperatures that match very well those measured in recent experiments for second sound in graphite [43]. Then, we rationalize the trend of FDN as a function of temperature in terms of the π terms entering in Eq. (22): π 1 increases with temperature for T < ∼ 100 K and saturates to a constant value at high temperature and π 3 decreases asymptotically like T −2 . Fig. 6b shows the analysis of FDN as a function of the sample's size l TOT (varied according to the aforementioned protocol) at fixed temperature 85 K, i.e. the temperature fixed in the experiments by Huberman et al. [43] to investigate the magnitude of hydrodynamic effects as a function of size. Also in this case, FDN predicts hydrodynamic effects to be largest in a size range that matches very well the sizes at which hydrodynamic effects have been measured. We thus conclude that the present model is capable of describing quantitatively the hydrodynamic window that has been recently measured in graphite. In Appendix I, we also show that for diamond extending the sizes in Fig. 5e up to 10 mm yields an increase of FDN for temperatures around 50-60 K and dimensions ∼ 1 mm. This result is in qualitative agreement with the predictions for diamond's second sound window performed using the reduced isotropic crystal model or the Callaway model [45]. Finally, the parameters entering in the π i (i = 1, 2, 3) are exactly the same used for the numerical solutions of the viscous heat equations in Fig. 5ac.

VII. CONCLUSIONS
We have introduced a framework to describe heat conduction beyond Fourier's law and to capture in the process the hydrodynamic transport regime where the phonon gas assumes a drift velocity and heat propagation resembles fluid dynamics. Just as a perturbation of the temperature generates an energy flux that is proportional to thermal conductivity, a perturbation of the drift velocity generates a crystal momentum flux, with the proportionality tensor coefficient between the two being a thermal viscosity. We have shown that thermal conductivity and viscosity can be determined exactly and in a closed form as a sum over relaxons (i.e. the eigenvectors of the phonon scattering matrix). Relaxons have a well-defined parity and even relaxons contribute exclusively to the thermal viscosity (odd relaxons contribute exclusively to the thermal conductivity [27]).
Most importantly, the microscopic LBTE has been coarse-grained into two novel mesoscopic viscous heat equations, which are coupled equations parametrized by the thermal conductivity and viscosity that govern the evolution of the temperature and drift-velocity fields. The viscous heat equations provide a general description of heat transfer encompassing both the Fourier and the hydrodynamic regimes, and all intermediate cases, allowing for the emergence of second sound, and of Poiseuillelike heat flow associated to a temperature profile deviating from Fourier's law. In addition, these viscous heat equations show that the thermal conductivity is not sufficient to describe thermal transport in general terms, but also its complementary quantity, the thermal viscosity, must be taken into account.
We have characterized the hydrodynamic behavior in terms of the Fourier deviation number (FDN), a dimensionless parameter that quantifies hydrodynamic deviations from Fourier's law. At a negligible computational cost, the FDN enables to study the temperature-and-size window at which hydrodynamic phenomena such as second sound emerge.
The full solution of the viscous heat equations allows to predict measurable temperature and heat flux profiles in complex shaped devices in the mesoscopic heat transport regime at a much reduced cost, and more transparently, than solving the full LBTE, therefore paving the way towards understanding shape and size effects in complex and microscopic devices, and especially when novel physics arises, such is the case in phononic devices [118][119][120][121][122][123][124][125][126][127][128][129][130]. These results are particularly relevant for the emerging field of materials that display hydrodynamic thermal transport, often characterized by large thermal conductivities [31,[39][40][41][42][43][44][45][46]. Such behavior has been investigated in devices made either of graphite, diamond or silicon. For graphite, the present formulation predicts the hydrodynamic windows (as a function of temperature for a fixed size, or as a function of size for a fixed temperature) in excellent agreement with recent experimental findings [43], while suggesting that hydrodynamic behavior can also appear in diamond at room temperature for micrometer-sized crystals.
Finally, we remark that the present methodology can be adapted to describe viscous phenomena for electronic conduction [131][132][133][134], or any other transport phenomena described by a linearized Boltzmann equation.

VIII. ACKNOWLEDGEMENTS
We acknowledge support from the Swiss National Science Foundation under project ID P2ELP2 168546 and the MARVEL NCCR. Computational resources have been provided by the Swiss National Supercomputing Center (CSCS) and PRACE.
The scattering matrix Ω ν,ν appearing in Eq. (1) is not symmetric, but it can be recast in a symmetric (and thus diagonalizable) form by means of the the following transformation [20,26,27]: where also the distribution n ν (r, t) appearing in Eqs. (6,7) has to be transformed for consistency. The symmetrized scattering operatorΩ ν,ν is real and symmetric, and can thus be diagonalized [20,26,27]: where θ α ν denotes a relaxon (i.e. an eigenvector), α is the relaxon index and the inverse eigenvalue τ α is the relaxon lifetime. We also define a scalar product [27]: used to orthonormalize eigenvectors. In order to show that eigenvectors with zero eigenvalues are related to conserved quantities in the LBTE dynamics, we rewrite the scattering operator distinguishing scattering events that conserve crystal momentumnormal (N) -from those that do not -Ukmlapp (U): As stated in the main text, there are four eigenvectors with zero eigenvalues, that we will discuss in the next sections.

The Bose-Einstein eigenvector: local temperature
Applying the transformation (A2) to Eq. (5) and considering the steady-state conditions, one obtains the following equation for the even part (see below for the odd part): The distributionñ T ν (r, t) is obtained applying the symmetrization (A2) to the distribution n T ν (r, t) appearing in equation (4) and it is thus evident thatñ T ν (r, t) ∝ ω ν (T (r, t) −T ). From the energy conservation in scattering events (both normal and Umklapp), it follows that [20,135] As a result, we identify ∂Nν ∂T as an eigenvector with zero eigenvalue θ 0 ν , that, after normalization, is where the specific heat C is From equation (A7), it follows thatñ T ν (r, t) disappears from Eq. (A6). Therefore, by removing the symmetrization (A2), Eq. (A6) gives Eq. (7) in the main text. In the context of the decomposition (A5), the Bose-Einstein eigenvector (A8) is an eigenvector to both the normal and Umklapp scattering operator, and will be denoted as φ 0 ν when we will later consider the basis of eigenvectors of the normal scattering operator in Appendix D.

The drift eigenvectors: local drift velocity
Starting from Eq. (5) at the steady-state, one obtains the following equation for the odd part: The distributionñ D ν (r, t) is obtained applying the symmetrization (A2) to the distribution n D ν (r, t) appearing in equation (4). We note in particular thatñ D ν (r, t) The drifting distributionñ D ν (r, t) is the stationary distribution for a system conserving crystal momentum. Therefore, recalling the decomposition (A5), we have that sinceΩ N ν,ν accounts only for Normal scattering events that conserve crystal momentum [13,16,20]. From Eq. (A15) it is possible to identify three eigenvectors of Ω N ν,ν with zero eigenvalues [20]: where i = 1, 2, 3 and A i is a normalization constant. The drift eigenvectors (A12) are, in general, not orthogonal [21]. Nevertheless, we work in a Cartesian coordinate system for q and u(r, t), so that these 3 eigenvectors are orthogonal and can also be normalized, choosing A i from the condition φ i |φ i = 1. In computing the normalization constants A i , we note that they can be expressed in terms of physically meaningful quantities. In particular, we note that the crystal momentum density associated to the drifting distribution is: and its derivative with respect to the drift velocity is: Comparing Eq. (A12) with Eq. (A14) we note that A i = ∂P i ∂u i eq and ∂P i ∂u j eq =0 ⇐⇒ i = j. Therefore, we will refer to A i as the specific momentum, due to its formal similarity with specific heat. It can be shown that in the high temperature limit, A i (T → ∞) ∝ T (see also Appendix E).
Finally, it is worth mentioning that in going from equation (5) to equation (6), the term n D ν (r, t) has disappeared due to the following approximation: This is correct both at high and low temperatures, because at high temperatures the strong crystal momentum dissipation ensuresñ D ν (r, t) ∝ u(r, t) ≈ 0 (see also Fig. 10), and at low temperatures Umklapp processes are frozen (Ω U νν ≈ 0).

Local equilibrium
From Eq. (A7) and Eq. (A15) it follows that, in the hydrodynamic regime, the distributions n T ν (r, t) and n D ν (r, t) are left unchanged by the dynamics described by the LBTE; therefore, these are local equilibrium distributions. It follows that n T ν (r, t) and n D ν (r, t) do not appear in Eq. (6) and Eq. (7) and thus do not contribute to the thermal conductivity and viscosity, which respectively describe the response to a perturbation of the local temperature and drift velocity. It is worth mentioning that, in the kinetic regime, n D ν (r, t) vanishes and n T ν (r, t) is still a stationary distribution for the LBTE.

Appendix B: Thermal viscosity
The total crystal momentum flux tensor Π ij tot [20] is defined as Due to the odd parity of q i and v j ν , only the even part of the phonon distribution contributes to the crystal momentum flux. Using the decomposition (4) introduced in the main text, we identify three contributions to the crystal momentum flux: whereΠ ij is the equilibrium (constant) crystal momentum flux, which is not affected by the LBTE's dynamics; Π ij T (r, t) is the momentum flux related to the local equilibrium temperature and Π ij δE (r, t) is the out-ofequilibrium momentum flux generated in response to deviations from local equilibrium conditions and is further discussed below.
As stated by Eq. (8), the thermal viscosity tensor η ijkl relates a drift velocity perturbation to the momentum flux generated as a response to that perturbation. Eq. (8) is valid in the mesoscopic regime mentioned in Sec. II and V, where the thermal viscosity is space-independent and non-homogeneities arise mainly from variations of the local temperature and local drift velocity. In particular, η ijkl is determined by deviations from local equilibrium and thus depends only on n δE . To determine the distribution n δE , we must solve the LBTE linearized in the drift velocity gradient Eq. (7). To this aim, we first symmetrize the LBTE using the transformations (A2), finding Next, using the relaxon approach discussed in Ref. [27] for thermal conductivity, we write the response to the perturbation ∇u as a linear combination of even eigenvectors:ñ Substituting this relation in the LBTE, and noting that the left term is related to the eigenvector φ i of the normal scattering matrix (Eq. (A12)), we obtain Taking the scalar product with a generic eigenvector θ α ν , we find where w j iα is a velocity tensor given by Thanks to the odd parity of φ i ν and v j ν , the velocity w j iα is different from zero only for even eigenvectors α; Eq. (B6) can thus be trivially solved for f ij α . With the knowledge of the LBTE solution f ij α at hand, the crystal momentum flux tensor is readily computed. We thus express Π ij δE in the relaxon basis, finding Substituting Eq. (B6) in Eq. (B8), we obtain the expression for the asymmetric thermal viscosity tensor: The tensor η ijkl is, in general, not symmetric under exchange of indexes j ↔ l. We will show later that, in a way analogous to fluids, only the sum over the spatial derivatives of the momentum flux tensor is relevant for the mesoscopic description of hydrodynamic thermal transport: where in the last term of Eq. (B10) we have rewritten the divergence of the momentum flux tensor using the symmetrized viscosity tensor µ ijkl : It is worth drawing a parallel between thermal viscosity and conductivity, where the latter can be written as [27] Notably, w i 0α is different from zero only for odd eigenvectors. As a result, thermal conductivity and viscosity are two quantities describing the transport due to the odd and even part of the spectrum respectively, i.e. energy and crystal momentum.

Single-mode relaxation-time approximation
In this section we derive the expression for thermal viscosity within the single-mode relaxation-time approximation (SMA). Using the SMA, the LBTE at Eq.
The deviation from equilibriumñ δE ν is readily found as We insert this result in the definition of momentum flux, obtaining From the definition of the thermal viscosity (8), the SMA thermal viscosity is therefore (B16) where we highlight the fact that in the SMA approximation η ijkl SMA = µ ijkl SMA , since in the SMA the viscosity tensor η ijkl SMA has the j↔l symmetry. A comparison between the exact bulk thermal viscosity (9) and the SMA bulk thermal viscosity (B16) is shown in Fig. 7. We highlight how the SMA approximation -which neglects the offdiagonal elements of the scattering operator and works well when the Umklapp processes dominate over normal processes [16,25,34] -overestimates the largest component of the thermal viscosity, especially at low temperatures. This overestimation is more pronounced in diamond compared to silicon, in agreement with results from previous works where the SMA approximation was reported to yield quite accurate results for the thermal conductivity of silicon [27] but not for diamond [25].
Appendix C: Size effects on the thermal conductivity and viscosity In order to obtain a simple estimate of size effects, we compute the total effective thermal conductivity and viscosity using a Matthiessen sum of the diffusive and ballistic limit: The ballistic conductivity and viscosity are computed for a sample of size L S as κ ij ballistic = K ij S · L S and µ ijkl ballistic = M ijkl · L S , where the prefactors K S and M ijkl are These prefactors are obtained after setting τ ν = L S |vν | in the SMA expressions of k and µ. L S is the length that determines boundary effects, i.e. the sample's size for a single crystal or the grain size for a polycrystalline sample. The numerical values of K ij S and M ijkl can be computed from first-principles and are tabulated in Tabs. I, II, III.

Appendix D: Viscous heat equations
In this section we derive an extension to Fourier's law from the LBTE, which describes mesoscopic hydrodynamic thermal transport in terms of the temperature T (r, t) and drift velocity u(r, t) fields. We focus on the mesoscopic regime where transport coefficients are well-defined and non-homogeneities arise from variations of the temperature and drift velocity. In contrast to Sec. II, here the temperature and drift velocity gradients are space-dependent, and consequently the LBTE is not linearized in the temperature and drift velocity gradients. We start recalling that hydrodynamic thermal transport emerges when most collisions between phonon wavepackets conserve the crystal momentum. This can happen, for example, when the mean free path for normal collisions Λ N is much smaller than the boundary scattering's length L S (for a single crystal L S is the sample's size, in the case of a polycrystalline sample L S is the size of the grains) or the mean free path for Umklapp collisions Λ U : Λ N L S , Λ U [15,31]. Under these conditions, the local equilibrium is expressed in terms of the four special eigenvectors φ 0 ν (also denoted θ 0 ν , since this is a common eigenvector for the full and normal scattering operator, see Appendix A), φ i ν (i = 1, 2, 3) described in Sec. (A 1,A 2), and of the local temperature T (r, t) and drift velocity u(r, t) fields. As explained in Sec. (A 1,A 2), these four special eigenvectors do not contribute to thermal conductivity and viscosity (i.e. they do not appear in Eqs. (6) and (7)), since these latter coefficients only describe the response to a perturbation of the local equilibrium.
In order to exploit the relationship between the drift velocity u(r, t) and the drift eigenvectors φ i ν , we choose to work with the basis of eigenvectors of the normal scattering matrixΩ N νν . To this aim, we diagonalizeΩ N νν as where β is an eigenvalue index, φ β ν an eigenvector, and 1 τ N β is an eigenvalue. Among these eigenvectors, we know the analytic expression for the 4 of them associated with energy and momentum conservation, which we label with β = 0 for the energy eigenvector (in this section we will use φ 0 ν to label the Bose-Einstein eigenvector (A8)), and β = 1, . . . , 3 for the momentum eigenvectors, Eq. (A12).
Noting that the set of "normal" eigenvectors {φ β ν } ("normal" in the sense that they diagonalizeΩ N νν , the part of the scattering matrix that accounts for normal processes only) is a complete basis set [16], we write the deviation from equilibriumñ ν (r, t) as a linear combination of these eigenvectors: After inserting Eq. (D2) in (1), we write the LBTE in the basis of the eigenvectors of the normal scattering operator This equation is formally equivalent to the LBTE, but allows us to take advantage of the knowledge of the first 4 eigenvectors to derive mesoscopic equations.

The projection of the LBTE on the 1 st (Bose-Einstein) eigenvector: energy moment
Here we show how to obtain an energy balance equation. First, we note from Eq. (4) that the phonon population expansion of Eq. (D2) can be recast as The coefficients z β in front of the 4 special eigenvectors ofΩ N νν are associated to temperature and drift velocity, which fully determine local equilibrium; in detail We now project the LBTE (D3) in the subspace spanned by the Bose-Einstein eigenvector φ 0 ν , i.e. we take the scalar product of Eq. (D3) with φ 0 ν , finding where we used the fact that φ 0 ν is an eigenvector of zero eigenvalue toΩ U νν (see Sec. (A 1)) and we defined the velocity tensor Note that the velocity W j αβ differs from the velocity w j αβ introduced in Sec. (II) for thermal viscosity; the difference arises from the use in Eq. (D8) of the "normal" eigenvectors (ofΩ N νν ) rather than the general eigenvectors ofΩ νν (W j αβ = w j αβ in presence of Umklapp processes).
Substituting Eqs. (D5, D6) in (D7) we obtain To elucidate the meaning of this equation, we note that the harmonic heat flux can be written as [136]: where we used the fact that only odd components of the phonon distribution contribute to the heat flux. Therefore, the heat flux receives contributions from both the drift velocity [103] and the temperature gradient [27]. In the basis of "normal" eigenvectors, the drifting contribu-tion can be written as while for the contribution from the deviation from local equilibrium n δO ν (r, t), we find where only odd eigenvectors contribute (W i 0β = 0 for even β eigenvectors). As explained in Sec. (A 3), the thermal conductivity is determined only from odd eigenvectors that are not related to local equilibrium (that is, all the odd eigenvectors minus the three drift eigenvectors). At this point it is crucial to recall that the relaxons have a well defined parity (even or odd), deriving from the symmetries of the full scattering operator Ω νν =Ω −ν,−ν [20,137]. Because the same symmetries apply toΩ N ν,ν , also the eigenvectors of the normal scattering operator have a well defined parity [20,21,137]. The odd component of the heat flux (and hence the thermal conductivity) is determined only from the odd part of the LBTE's solutionñ δO ν (r, t) [27,137], and this can be written either as a linear combination of odd relaxons (as done in Ref. [27]) or, equivalently, as a linear combination of odd eigenvectors of the normal scattering matrix, as done here.
From the completeness of these basis sets it follows that the heat flux (D12) arising from the odd out-of-equilibrium phonon distribution determined from equation (6), and described as a linear combination of odd eigenvectors of the normal scattering matrix, is equivalent to the heat flux of Ref. [27] that is written as a linear combination of odd relaxons (eigenvectors of the full scattering matrix). In an equivalent way, the equation for the odd part of the LBTE (6) is equivalent to the equation used to determine the thermal conductivity in Ref. [27], the only difference being the basis chosen for the decomposition of the solution (the relaxons in Ref. [27] and in the computation of the thermal viscosity reported in Fig. 1, and the eigenvectors of the normal scattering matrix here). The choice of the basis set does not affect the resulting local heat flux, which can therefore be related to the local temperature gradient via the thermal conductivity [27]: In Eq. (D13) the thermal conductivity is spaceindependent, i.e. Eq. (D13) is valid in the mesoscopic regime where the thermal conductivity and viscosity are space-independent and non-homogeneities arise solely from variations of the local temperature and local drift velocity. Eq. (D13) can be used to rewrite Eq. (D9) in terms of the local temperature and drift velocity fields, obtaining Eq. (10) of the main text: Eq. (D14) is clearly not sufficient to fully describe the hydrodynamic heat conduction problem in which both T (r, t) and u(r, t) are nonzero. In the next section we will derive a complementary set of equations that completes the formulation.
2. The projection of the LBTE on the 2 nd /3 rd /4 th eigenvectors (the momentum eigenspace) In this section, we derive a set of balance equations for crystal momentum. We start by recalling from Eq. (B2) that the momentum flux receives contributions from three different terms. Of these three, the first term is a constant related to the equilibrium temperature, and thus is not changed by the LBTE. Therefore, we focus only on the momentum flux related to the local equilibrium temperature Π ij T (r, t) and the out-of-equilibrium momentum flux generated in response to a drift velocity gradient Π ij δE (r, t). Using the expression of the four special eigenvectors discussed in Sec. (A 1, A 2), we rewrite these two momentum fluxes in the basis of eigenvectors of the normal scattering matrix, finding: and where we used the velocity tensor defined in Eq. (D8). Next, as in the previous section, we take the scalar product of Eq. (D3) with φ i ν (i = 1, 2, 3), obtaining the following three equations indexed by i = 1, 2, 3 where we used the fact that φ i ν are eigenvectors ofΩ N νν with zero eigenvalues and we defined From the propertyΩ U ν,ν =Ω U −ν,−ν it can be shown that D iβ U vanishes when β indexes an even eigenvector. Since the coefficients z β (r, t) for the first four eigenvectors (β = 0, 1, 2, 3) are known, it is convenient to rewrite Eq. (D17) as: In the hydrodynamic regime, Umklapp momentum dissipation is weak and thus D iβ U → 0 (φ i ν is approximately an eigenvector with a vanishing eigenvalue forΩ U νν ). Therefore, we simplify Eq. (D19) noting that Then, we use the expression of the coefficients z β (r, t) (β = 0, 1, 2, 3), and substitute Eqs. (D5,D6) in the simplified Eq. (D19), obtaining: Next, we note that only even eigenvectors different from the Bose-Einstein eigenvector determine the even distribution n δE ν (r, t) appearing in the expression for the outof-equilibrium momentum flux tensor (D16). As shown in Eq. (B10), in the mesoscopic regime the sum over the spatial derivatives of Π ij δE (r, t) can be expressed in terms of the viscosity and second derivative of the drift velocity. We thus find Combining this with Eq. (D14) we obtain 4 equations to be solved in terms of temperature T (r, t) and drift velocity u(r, t); these are the viscous heat equations at the core of this work and that are further discussed in the main text. As a final remark it is worth mentioning that in the kinetic regime, characterized by strong crystal momentum dissipation, the inequality (D20) may not be valid. Nevertheless, in such regime the strong crystal momentum dissipation (max ij [D ij U ] → ∞) yields the following (stronger) inequality implying that the viscous heat equations reduce to Fourier's law, as discussed in Sec. VI. From a mathematical point of view, the projection of the LBTE in the Bose-Einstein subspace, performed computing the scalar product between the LBTE in the normal eigenvectors basis (D3) and the Bose-Einstein eigenvector (A8) ∝ ω ν , is equivalent to calculating the energy moment of the LBTE. Analogously, the projection in the momentum subspace, performed calculating the scalar product between the LBTE (D3) and the momentum eigenvectors (A12) ∝ q i , is equivalent to computing the momentum moment of the LBTE.
Appendix E: Parameters entering the viscous heat equations.
We report in Fig. 8 the trends as a function of temperature of all the parameters needed to solve the viscous heat equations on the geometry of Fig. 3. Due to the crystal symmetries of the materials considered, several components of vectors and tensors are equivalent. 2 ndrank tensors such as thermal conductivity (κ ij ), momentum dissipation (D ij U ) and velocity (W j 0i ) are isotropic and diagonal for diamond and silicon; for graphite they are isotropic and diagonal for in-plane (x, y) directions (as used in the geometry of Fig. 3). Similarly, the spe- cific momentum A i is isotropic for diamond and silicon, and independent from the in-plane direction in graphite (A x = A y ). The values of D ij U , A i , C and W j 0i as a function of temperature are plotted in Fig. 8 (we only report the in-plane or isotropic components, indexes are omitted). The numerical values of these parameters can be found in Tabs. I, II, III; in addition, these tables contain also the prefactors discussed in Sec. II, III and Appendix C needed to compute the thermal conductivity in the ballistic limit κ ij ballistic . In particular, the parameter K ij S (defined for silicon, diamond and in-plane graphite), is defined as K ij S = κ ij ballistic /L S , where L S is the length that determines surface effects, i.e. the sample's size for a single crystal or the grain size for a polycrystalline sample. In the tables, we also report the non-zero irreducible components of the bulk heat flux viscosity tensor, limiting graphite to in-plane components. The parameters M ijkl are the prefactors (defined in Appendix C) used to compute the ballistic thermal viscosity as µ ijkl ballistic = M ijkl ·L S . Similarly, the parameter F ij U are defined as F ij U = D ij U,boundary (L S ) · L S . Additionally, Fig. 9 shows the off-plane transport coefficients (panel a) and the off-plane non-negligible values of A i , D ij U , and W j 0i (panel b) for graphite; these parameters are not relevant for the geometry studied in Sec. V and are reported for completeness.
wards, we will analyze second sound within the viscous heat equations following a top-down approach, i.e. requiring the temperature field to have the mathematical form of a damped wave and analyzing if and under which conditions this can emerge from the viscous heat equations.

Second sound from the viscous heat equations (bottom-up approach)
For simplicity, we consider a system such that the tensors W i 0j and κ ij appearing in the viscous heat equations (10,11) are isotropic; the generalization to an anisotropic case is analogous to what is reported here.
Without loss of generality, we considerx as the direction of second sound propagation. For simplicity we consider an isotropic system; the general derivation can be obtained straightforwardly generalizing the procedure reported here. In an isotropic system, the drifting heat flux Q D (r, t) is collinear with the drift velocity and the heat flux due to local temperature changes Q δ (r, t) is collinear with the temperature gradient. Thus, it follows that the only nonzero component of the drift velocity must be along the second sound propagation direction u x = u (for simplicity we omit all tensor indexes in the rest of this section, since the only the component having all the indexes equal to x is needed for this discussion). With these conditions, the viscous heat equations (10,11) become: In order to observe second sound, temperature changes need to propagate following a damped wave equation. To this aim, we require that drift velocity and temperature are related as (F3) where τ ss is the second sound relaxation time and 0<|f |<1 is a constant, both to be determined. To better understand this requirement, we insert Eq. (F3) in Eq. (F1), finding the desired temperature damped-wave equation: (F4) Next, we show that condition (F3) implies that also the drift velocity field follows a damped-wave equation.
To this aim, we take the derivative with respect to x of Eq. (F3), finding whereÔ is a differential operator. Next, by applying the operatorÔ to both sides of equation (F2) and using condition (F5), we obtain ∂t .

(F6)
If we consider only the lowest-order derivatives in Eq. (F6), we obtain a simplified equation that holds in the close-to-equilibrium regime where variations in space and time are small. In particular, neglecting higher-thansecond order derivatives gives: Therefore, if then both constants c 1 and c 2 are positive and the evolution of the drift velocity is that of a damped wave equation like the one for temperature.
The coefficients τ ss and f are determined solving Eqs. (F4,F7) and imposing the second sound condition (F3). The solutions are of the form and The condition (F13) is derived from the requirement c 1 = 1 τss , and is consistent with the damped wave requirement (F8); condition (F14) is derived from the requirement that c 2 = κ Cτss(1−f ) and the second sound velocity (F15) has been obtained substituting Eq. (F13) and Eq. (F14) into Eq. (F4). We note that for D U → 0 (that is, negligible crystal momentum dissipation) v g (k) → v ss → W , i.e. the second sound propagation velocity approaches the drifting second sound velocity defined by Hardy [20]) and τ ss → D U −1 . Finally, the second sound condition (F3) imposes the following relation between the coefficients: C T (k) must be set according to initial conditions and the form of C u (k) follows from equation (F16). We note from Eq. (F16) that, when second sound occur, temperature and drift velocity are both damped waves with a phase shift of π/2.

Second sound from the viscous heat equations (top-down approach)
In the previous section, we obtained second sound properties by finding the conditions under which the viscous equation for temperature (Eq. (10)) becomes a damped wave equation. However, we can also obtain a second sound equation taking inspiration from the approach outlined in Ref. [40], i.e. by looking for the conditions upon which the microscopic degrees of freedom of the transport equation evolve as a damped wave. In particular, we want to find the conditions such that the solution of Eqs. (10,11) are u(x, t) = u 0 e i(kx−ω(k)t) e −t/τss , where δT and u 0 are in general complex numbers to allow for a phase difference between the two waves. We note in particular that this guess for solution requires that both temperature and drift velocity oscillate/decay at the same frequency/rate, which is consistent with the conditions (F9,F10) obtained in the previous section. Using this guess for the solution, the derivation of the dispersion relation and the decay time easily follows. To this aim, we substitute Eqs. (F17,F18) in the viscous heat equations (10,11) and find: where we introduced a complex frequencyω(k) =ω(k) − i τss to simplify the calculation. The real part of this complex frequency is the oscillation frequency of second sound, whereas the imaginary part describes its decay time. Next, we rewrite Eq. (F20) as: and substitute this expression into Eq. (F19), finding that gives: (−iCω(k) + κk 2 )(AD U + µk 2 − iAω(k)) + CAW 2 k 2 = 0 . (F23) This is a quadratic equation that determines the dispersion relations forω k , given by: This equation can be solved to obtain the complex frequencyω(k) and thus the oscillation frequency and decay time of second sound as a function of the wavevector k. Solving for this quadratic equation, we obtain: In order to compare this result with the expression for second sound derived in the previous section, it is worth recalling that the semiclassical description of thermal transport used throughout this work holds for longwavelength perturbations. Therefore, we simplify the previous expression retaining terms to smallest order in k, finding: We can readily see that in the limit of small wavevectors (k → 0) the non-trivial solution isω k ≈ −iD U , that is, the second sound oscillation frequency goes to zero, and has a decay time set by the crystal momentum dissipation rate: τ ss ≈ 1 D U . To estimate the behavior of the oscillation frequency, is instead convenient to recall the hypothesis of small Umklapp rates. In fact, if we set D U = 0, we findω(k) ≈ ±W k, that is, second sound disperses linearly with k, and has a velocity v g (k) = ∂ω k ∂k ≈ W . These limiting approximations (v g ≈ W and τ ss ≈ 1 D U ) are consistent with the results found in the previous section.

Appendix G: Analytical 1D example
The viscous heat equations (10,11) can be solved analytically in a handful of toy models. Here, for simplicity, we neglect dissipation of momentum by Umklapp processes (D ij U = 0) and consider steady-state heat diffusion along the transversal direction of a thin film, so that the problem becomes effectively a 1D problem, with x labeling the orthogonal direction position. In addition, since we are considering a 1D geometry, we label A = A x , W = W x 0x , µ = µ xxxx and assume that temperature and velocity fields depend only on the position x. Under these considerations the viscous heat equations (10,11) reduce to: To solve the problem, we specify the following no-slip boundary conditions on a 1D geometry having length 2l: and that is, we assume boundaries at thermal equilibrium.
We look for solutions of the form: T (x) =T + c sinh bx .
After some algebra, one finds the solution This analytical solution shares several qualitative similarities with the numerical example discussed in the main text, and more clearly highlights how the factor 1/b represents a length scale at which surface scattering affects thermal transport, which is in turn dependent on both conductivity and viscosity. Moreover, we note that the mathematical form of the solution has the same qualitative behavior of the problem studied by Sussmann and Thellung [13], which serves as a verification of the present model. At variance with their work however, the prefactors introduced here allow us to go beyond the Debye approximation.
Appendix H: Estimate of the characteristic drift velocity In this section, we estimate the characteristic value of the drift velocity (u 0 ) in the high (u H ) and low (u L ) temperature regimes. These characteristic values are determined substituting in the viscous heat equations (10,11) the characteristic values of the temperature (and related derivatives) and solving them approximatively for the velocity. With this aim, we start estimating the characteristic temperature gradient in the setup of Fig. 3 when a temperature differenceT ± δT is imposed on the two opposite sides (at x = 0 and x = 15 µm).
In this setup we clearly distinguish two regions. Choosing l = lTOT 10 as a length unit, where l TOT is the total length of the sample (e.g. l TOT = 15µm for the sample in Fig. 3), the left-hand-side region has length l L = 2 l and width w L = 2 l, while the right-hand-side region has length l R = 8 l and width w R = 6 l. Energy conservation requires that the current in the left-hand side must be equal to the right-hand side. Therefore, the heat flux on the left Q L must be three times the heat flux on the right side Q L = 3Q R . Using Fourier's law Q = −k∇T , and supposing that the thermal conductivity is constant throughout the sample, it follows that the temperature gradients in the two regions are related as ∇ x T L = 3∇ x T R . Requiring the total temperature drop to be equal to the temperature difference imposed by the boundary conditions, we can write Focusing from now on on the larger region on the right, it follows that the temperature drop taking place is approximately given by ∆T R = − 8 7 δT . To determine a characteristic value of drift velocity u 0 , we substitute ∇T x R ∆T R l R in Eq. (11) and for simplicity consider isotropic symmetry, valid e.g. for silicon, diamond and for graphite in the in-plane directions. u 0 can be determined focusing on the steady state limit of Eq. (11). We simplify its estimate considering separately the limits of low and high temperatures.
In the high temperature limit, the term related to momentum dissipation (∝ D ij U ) is much larger than the viscous term (∝ µ ijkl ) (see Fig. 8), therefore Eq. (11) can be approximated as: Using the estimated temperature gradient ∇T x R ∆T R l R , the high-temperature characteristic value of drift velocity u H is found to be At low temperatures, viscosity dominates over the momentum dissipation term (see Fig. 8), so that Eq. (11) is approximated as ∂T (x, y) ∂x µ xxxx ∂ 2 u x (x, y) ∂x 2 +µ xyxy ∂ 2 u x (x, y) ∂y 2 , (H4) where we considered only the two largest components of the viscosity tensor. To estimate the average value of these second derivatives, we note that, as shown in Fig. 3, u(x, y) has a bell-like profile in the sample interior, which vanishes at the boundaries. We thus proceed with a few assumptions that allow us to make an estimate of u L . First, we suppose that the average values of the second derivatives of the drift velocity in Eq. H4 are of the same order of magnitude, that is ∂ 2 u x ∂x 2 ∼ ∂ 2 u x ∂y 2 , which can be checked numerically. It follows that Eq. (H4) can be simplified as: Next, we note that the variation of u is stronger along the y coordinate. To mimic the Poiseuille-like shape, we assume the velocity profile to be constant along the x direction, and parabolic along the y direction with vanishing velocity at the boundaries (y = 0 µm and y = w R ), so that u(x, y) (−a · y(y − w R ), 0, 0). With these approximations, we can estimate the average value of the parabolic velocity profile, i.e. the characteristic value of the drift velocity at low temperatures u L , as: where we used w R = 6 l and l R = 8 l. We thus recover the expression for u L given in section VI in the main text. The average value of u is interpolated in between the high and low temperature limit using Matthiessen's rule, as: u −1 0 = u −1 H + u −1 L . The estimates of u H , u L and u 0 for graphite are reported in Fig. 10. We also compare this rough estimate with the average value of u x (x, y) computed from the numerical solution of the viscous heat equations on the region x > 1 5 l TOT of the geometry discussed in the main text and denoted denoted with u x . Despite the qualitative arguments used to derive u 0 , the estimate is able to capture qualitative trends and approximately reproduce average results from the numerical solution of the viscous heat equations.

Appendix I: Hydrodynamic behavior in macroscopic diamond samples
In this section we show how extending the sizes in Fig. 5e up to 10 mm yields another peak for FDN for dimensions around 1 mm and temperatures around 50-60 K. This result has been confirmed by full solutions of the viscous heat equations similar to these reported in Fig. 5b. The increase of FDN at low temperature and large size reported in Fig. 11 is in qualitative agreement with the predictions for the second-sound window performed using the reduced isotropic crystal model or the Callaway model [45]. However, in contrast with this work, Ref. [45] reports that an isotopic concentration much lower than the natural isotope abundance is necessary for the observation of second sound in diamond. The analysis here and in Fig. 5 are limited to a minimum reference temperature ofT = 50 K because, for lower temperatures, anharmonicity becomes decreasingly smaller while size effects become increasingly larger, resulting in convergence issues for the numerical calculations (which would also require a refined, computationally-expensive, treatment of surface effects [35]). Most importantly, the temperature/size range reported in Fig. 5 highlights the most accessible experimental conditions under which hydrodynamic behav-ior can appear in diamond, namely micrometer-sized diamond samples (thus much less expensive than millimetersized samples) and non-cryogenic temperatures.

Appendix J: Computational details
First-principles calculations are performed with the Quantum ESPRESSO distribution [138,139]. For all the materials analyzed, the LDA functional is used due to its capability to accurately describe the structural and vibrational [140] properties of graphite [69], diamond [25] and silicon [27], and its compatibility with the D3Q code [32,33] for first-principles calculations of anharmonic (third-order) interatomic force constants. Details on the computation of the second-and third-order interatomic force constants are reported in Ref. [27] for silicon, in Ref. [25] for diamond and in Ref. [69] for graphite. The LBTE's scattering matrix Ω νν is computed as in Ref. [25] and accounts for third-order anharmonicity [32] and isotopic disorder [141,142] at natural abundance. Thermal conductivity and viscosity calculations for silicon and diamond are performed using 27×27×27 q-point grids and a Gaussian smearing of 2 cm −1 and 8 cm −1 , respectively. Thermal conductivity and viscosity calculations for graphite are performed using a 49×49×3 qpoint grid and a Gaussian smearing of 8 cm −1 . The use of an odd, Gamma-centered q-points mesh is crucial to correctly account for the parity symmetries of the scattering operator. I. Parameters entering the viscous heat equations for graphite. We report here only the in-plane components of the tensors needed to perform the calculation of Fig. 3: κ ij P = κPδ ij , κ ij C = κCδ ij , K ij S = KSδ ij , D ij U = DU δ ij , F ij U = FU δ ij , A = A i ∀ i, W j 0i = W j i0 = W δ ij , where the indexes i, j represent the in-plane directions x, y only (i, j = 1, 2).  TABLE II. Parameters entering the viscous heat equations for diamond. As discussed above, due to the symmetries of diamond's crystal, κ ij P = κPδ ij , κ ij C = κCδ ij , K ij S = KSδ ij , D ij U = DU δ ij , F ij U = FU δ ij , A = A i ∀ i, W j 0i = W j i0 = W δ ij , where i, j = 1, · · ·, 3.  III. Parameters entering the viscous heat equations for silicon. As discussed above, due to the symmetries of silicon's crystal, κ ij P = κPδ ij , κ ij C = κCδ ij , K ij S = KSδ ij , D ij U = DU δ ij , F ij U = FU δ ij , A = A i ∀ i, W j 0i = W j i0 = W δ ij , where i, j = 1, · · ·, 3.