Specifically, Ref. [13] observed temperature oscillations at temperatures as high as 200 K in isotopically purified graphite, while the pioneering works [11, 12] observed temperature oscillations around 80–125 K in graphite at natural isotopic abundance.

Implemented in bte-barna following Sec. 4.2.2.1 of [54].

One can also directly see that $⟨{\theta }_{\nu }^0|v_{\nu }^i|{\theta }_{\nu }^0⟩\propto {\sum }_s{\int }_{\mathfrak{B}}{\omega }_{\mathit{q}s}^2{v}_{\mathit{q}s}^i{d}^{3}q=0$, because it is obtained by integrating an odd-parity integrand over the symmetric Brillouin zone $\mathfrak{B}$.

The relaxon [20] and variational [18] methods always provide results that are numerically compatible within 0.1%. For example, in natural graphite at 80 K the bulk thermal conductivity obtained from the variational method [18] at 80 K is $5735.08\mathrm{W}/\mathrm{mK}$ while with the relaxon formalism it is 5735.03 W/mK; in isotopically purified graphite (99.9% ${}^{12}\mathrm{C}$, 0.1% ${}^{13}\mathrm{C}$) the bulk thermal conductivity at 80 K obtained from the variational method is $27484.96\mathrm{W}/\mathrm{mK}$ while with the relaxon formalism it is 27485.09 W/mK; for $\mathrm{h}{}^{11}\mathrm{BN}$ the variational method at 60 K yields $3944.20\mathrm{W}/\mathrm{mK}$, while the result of the relaxon method is $3944.23\mathrm{W}/\mathrm{mK}$.

Specifically, Eq. (0.7.3.1) of Ref. [80] shows that for a partitioned nonsingular matrix $A$, $$A=\left[\begin{array}{ll}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right]$$(3)its inverse can be written in partitioned form as, $$A^{-1}=\left[\begin{array}{lc}{(A_{11}-A_{12}{A}_{22}^{-1}{A}_{21})}^{-1}& A_{11}^{-1}{A}_{12}{(A_{21}{A}_{11}^{-1}{A}_{12}-A_{22})}^{-1}\\ A_{22}^{-1}{A}_{21}{(A_{12}{A}_{22}^{-1}{A}_{21}-A_{11})}^{-1}& {(A_{22}-A_{21}{A}_{11}^{-1}{A}_{12})}^{-1}\end{array}\right].$$(4)

The volume viscosity is also called the bulk viscosity, here we use the term “bulk viscosity” for the value of the viscosity tensor in the bulk of the crystal, without accounting for finite size effects.

Strictly speaking, in proximity of a perfectly thermalized boundary (temperature fixed at a certain value and $\mathit{u}=\mathbf{0}$), the temperature gradient emerging from the VHE is weakly space dependent, due to the coupling between the temperature gradient and the drift velocity in Eq. (2). In practice, in a rectangular geometry of graphite with gradient applied around 70 K, the temperature gradient and drift velocity reach a practically space-independent value within about $1\mathrm{\mu }\mathrm{m}$ from the boundaries. Analytical details on the mechanisms determining this length scale can be found in Appendix G of Ref. [37], and numerical simulations in Fig. 4 of the same reference.

For a well-behaved vector field, the circulation computed on a closed streamline $\overrightarrow{\mathcal{S}}$ is nonzero: $\mathcal{C}={\oint }_{\mathcal{S}}{\mathit{Q}}^D·d\overrightarrow{\mathcal{S}}\ne 0$; therefore, one sees from Stokes’s theorem that this implies nonzero vorticity: $\mathcal{C}={\int }_{\mathcal{A}}(\nabla \times {\mathit{Q}}^D)·\mathit{m}d\mathcal{A}\ne 0$, where $\mathit{m}$ is the unit vector outward normal to the surface $\mathcal{A}$ enclosed by $\overrightarrow{\mathcal{S}}$.

Specifically, in the region where temperature inversion is most pronounced in Fig. 1, one has $\beta \{[\partial \mathit{T}(\mathit{r},t)]/(\partial r^y)\}\sim {\eta }^{yyyy}\{[{\partial }^2{u}^y(\mathit{r},t)]/(\partial r_y^2)\}$, a condition that permits temperature gradient and drift velocity to be aligned, thus temperature-gradient heat flux ${\mathit{Q}}^{\delta }=-{\kappa }_D\nabla \mathit{T}$ and drifting heat flux ${\mathit{Q}}^D=\alpha \mathit{u}$ to have opposite directions.

The VHE solution in Fig. 1 was computed in less than two minutes on a laptop, while the corresponding LBTE solution in Fig. 5 took about 72 hours on one node (128 cores) of the Kelvin2 supercomputer from the Northern Ireland High Performance Computing center.

We note that the VHE have a mathematical form analogous to the damped and diffusion-extended Navier-Stokes equations, used to describe rarefied fluids flowing in porous media in isothermal and laminar conditions, with VHE’s temperature mapping to fluid’s pressure and VHE’s drift velocity mapping to fluid velocity. To see this, we highlight how the linear damping term $-{\gamma }^{ij}{u}^j(\mathit{r},t)$ appearing in the second VHE (2) is analogous to the dissipative term appearing in the linearized “damped” Navier-Stokes equations [120, 121, 122, 123, 124] for fluids flowing through a porous medium. Then, the diffusive term appearing in the first VHE (1) $-{\kappa }_D^{ij}\{[{\partial }^2\mathit{T}(\mathit{r},t)]/[\partial r^i\partial r^j]\}$ is analogous to the self-diffusion term appearing for pressure in the extended Navier-Stokes equations [85, 125, 126] used to describe the flow of an isothermal, compressible, and rarefied gas [127].

The DPLE encompasses Cattaneo’s second-sound equation [128] as a special case, as discussed in SM II and [42].

We investigate monoisotopic $\mathrm{h}{}^{11}\mathrm{BN}$ because Refs. [6, 129] suggest that heat hydrodynamics in $\mathrm{h}{}^{11}\mathrm{BN}$ is stronger than in hexagonal boron nitride with natural isotopic-mass disorder (19.9% ${}^{10}\mathrm{B}$ and 80.1% ${}^{11}\mathrm{B}$).

We recall that the curl of a gradient is zero, so Fourier’s heat flux $Q^{\delta ,i}=-{\kappa }_D^{ij}{\nabla }^j\mathit{T}$ is irrotational.

In this time-dependent simulation we choose an equilibrium temperature of 80 K to match the temperature at which transient hydrodynamic heat propagation has been observed in recent experiments [11, 12].

$\dot{q}(\mathit{r},t)=\mathcal{H}·\theta (t_{\mathrm{heat}}-t)\exp \{-[(x+x_c{)}^2/(2{\sigma }_x^2)]-[y^2/(2{\sigma }_y^2)]\}$. To ensure that the perturbation created causes variations within 10% of the equilibrium temperature, we used the following parameters: $\mathcal{H}=0.013[\mathrm{W}/(\mathrm{\mu }{\mathrm{m}}^3)]$, $t_{\mathrm{heat}}=0.4ns$, $x_c=5\mathrm{\mu }\mathrm{m}$, ${\sigma }_x=2\mathrm{\mu }\mathrm{m}$, ${\sigma }_y=2.8\mathrm{\mu }\mathrm{m}$.

$\dot{q}(\mathit{r},t)=\mathcal{H}·[\sin (2\pi ft)+1]\exp \{-[(x+x_c{)}^2/(2{\sigma }_x^2)]-[y^2/(2{\sigma }_y^2)]\}+\mathcal{H}·[\sin (2\pi ft+\pi )+1]\exp \{-[(x-x_c{)}^2/(2{\sigma }_x^2)]-[y^2/(2{\sigma }_y^2)]\}$.

The simulations in Fig. 3 were performed accounting for grain-boundary scattering as in Ref. [37], and considering a grain size of $20\mathrm{\mu }\mathrm{m}$. This value was estimated considering the largest grains in Fig. S3 of Ref. [11], Fig. S1 of Ref. [12], Fig. S11 of Ref. [13], and Fig. 1(e) of Ref. [6], which are all in broad agreement with $20\mathrm{\mu }\mathrm{m}$.

Experiments [13] for isotopically purified graphite are available only at temperatures higher than 100 K, preventing us from comparing VHE and DPLE in the low-temperature limit, where Figs. 4 and 4 show that viscous effects are largest.

Refs. [11,13] quantified the hydrodynamic strength as the dip depth of the transient thermal grating signal, see, e.g., Fig. 1(a) in Ref. [13].