Thermal transport from 1D- and 2D-confined nanostructures on silicon probed using coherent extreme UV light: General and predictive model yields new understanding

Heat management is crucial in the design of nanoscale devices as the operating temperature determines their efficiency and lifetime. Past experimental and theoretical works exploring nanoscale heat transport in semiconductors addressed known deviations from Fourier's law modeling by including effective parameters, such as a size-dependent thermal conductivity. However, recent experiments have qualitatively shown new behavior that cannot be modeled in this way. Here, we combine advanced experiment and theory to show that the cooling of 1D- and 2D-confined nanoscale hot spots on silicon can be described using a general hydrodynamic heat transport model, contrary to previous understanding of heat flow in bulk silicon. We use a comprehensive set of extreme ultraviolet scatterometry measurements of non-diffusive transport from transiently heated nanolines and nanodots to validate and generalize our ab initio model, that does not need any geometry-dependent fitting parameters. This allows us to uncover the existence of two distinct timescales and heat transport mechanisms: an interface resistance regime that dominates on short time scales, and a hydrodynamic-like phonon transport regime that dominates on longer timescales. Moreover, our model can predict the full thermomechanical response on nanometer length scales and picosecond time scales for arbitrary geometries, providing an advanced practical tool for thermal management of nanoscale technologies. Furthermore, we derive analytical expressions for the transport timescales, valid for a subset of geometries, supplying a route for optimizing heat dissipation.


I. INTRODUCTION
Advances in fabrication have scaled the characteristic dimensions of complex systems to the few nanometer range and even thinner. At these length scales, conventional macroscopic (bulk) models can fail to accurately describe nanoscale behavior, due to the dominance of interfaces and surfaces. Specifically, thermal transport from nanoscale heat sources on semiconductor substrates strongly deviates from bulk diffusive transport predictions.
Experiments show that as the heat source size is reduced, the heat transport efficiency falls well below what is predicted by bulk diffusion, both for structured optical excitation [1][2][3][4] or for optically-excited nanostructured transducers [5][6][7][8][9][10][11][12] . Moreover, recent experiments have uncovered that both the size and spacing of periodic nanoheater arrays strongly influence thermal transport, resulting in new counter-intuitive behaviors 8,9 . However, there is still no consensus on the underlying physics-in large part because there is no comprehensive model to describe these new nanoscale thermal transport regimes. This precludes smart design for good thermal management in next-generation nanodevices.
New theoretical proposals based on truncated Levy flights 13 , suppression of phonons 8,14 or relaxons 15 have explained certain aspects of non-diffusive thermal transport for specific geometries. Phonon hydrodynamics [16][17][18][19][20][21][22][23][24] has been also successfully used to explain thermal transport behavior on 2D materials 25 , such as graphene 26 , and even in bulk materials at very low temperatures 27 . As this behavior is known to occur when "normal" phonon scattering events (i.e. processes that conserve quasi-momentum) dominate over "resistive" ones, the existence of hydrodynamic transport in bulk semiconductors, like silicon, at room temperature has been explicitly discarded 28 .
It is widely accepted that solving the Boltzmann transport equation with ab initio calculated parameters is the most precise way to describe the transport of phonons, which are the dominant heat carriers for semiconductor and dielectric materials 29-31 . However, several difficulties associated with this approach limit its use at a practical level. First, this equation is challenging to solve in general due to the complexity of phonon collisions. To overcome this challenge, the relaxation time approximation is often used to simplify the collision expression 13,32 . However, this approximation does not guarantee energy conservation, which can lead to invalid results 13,33 . Second, complex geometries are challenging because one must model how each phonon mode interacts with every boundary present. Finally, coupling this equation to other phenomena, such as thermoelectricity or thermoelasticity, exponentially increases the computational requirements. These challenges can prevent microscopic models, like the Boltzmann transport equation, from being directly compared to experimental data.
Because microscopic models cannot be applied in many complex geometries, experiments often use an intermediate layer, or mesoscopic models, to compare results and theory. Most mesoscopic models to date are based on Fourier's law of heat diffusion with the addition of phenomenological effective parameters. This approach fits effective parameters to experiments, and then formulates theoretical models to connect the fitted values to ab initio calculations. Recent works have used this effective Fourier model to analyze heat dissipation away from metallic nanostructures of varying size and spacing [5][6][7][8][9]11 . This can quantify the deviation from the diffusive prediction by fitting either an effective thermal boundary resistance between the transducer and substrate [7][8][9] or an effective thermal conductivity of the substrate 5,6,12,34-36 . These techniques have significantly advanced our understanding, making it possible to develop new experimental mean free path spectroscopy techniques 1 , as well as uncovering new transport regimes dominated by the heat source spacing 7,8 . However, using Fourier's law as a mesoscopic model, even with effective parameters, can obscure the underlying physics and fails to predict thermal transport observed for all time and length scales 19,22 . Most importantly, this approach is difficult to generalize to arbitrary geometries or materials.
In this work, we present a comprehensive set of dynamic EUV scatterometry measurements of non-diffusive heat flow away from 1D-and 2D-confined nanostructures on bulk silicon. We use this data to validate and generalize the Kinetic Collective Model (KCM) 37,38 , which is a mesoscopic model which uses a hydrodynamic-like heat transport equation 16 with ab initio parameters. Contrary to conventional understanding, we show that heat transport away from nanoscale sources on bulk silicon can be predicted by the hydrodynamic equation. This generalizes the hydrodynamic framework to situations where phonon momentum is conserved, which applies not only when normal collisions dominate, but in regions with size comparable to the average resistive phonon mean free path near heat sources and system boundaries 22,38,39 . We also experimentally observe for the first time that closely-spaced 2D-confined (nanodots) on a bulk silicon substrate cool faster than widelyspaced ones, and that this effect is larger in 2D-confined than in 1D-confined (nanoline) sources observed by previous works 8,9 . Moreover, we demonstrate that KCM both fully predicts the heat transport over a wide range of length-scales and time-scales from 1D-and 2D-confined heat sources on a silicon substrate-including the counter-intuitive behavior of the closely-spaced geometry-and captures the full thermomechanical response to the system, which is beyond the capabilities of microscopic models. Our mesoscopic hydrodynamic model also provides insight into the fundamental transport behavior. KCM allows us to identify the time scales over which two different transport mechanisms are dominant: one characteristic time dominated by the thermal boundary resistance and another regime that is dominated by hydrodynamic heat transport.
The latter mechanism is responsible for the slow thermal decay of small heat sources, and consequently its reduction is responsible for the increased dissipation of close-packed nanoheaters. Furthermore, we develop a two-box model, derived from the hydrodynamic equation, which provides a physical interpretation and specific expressions for the two characteristic dissipation mechanisms. We confirm these findings by comparing our models to both past 1D-confined and new 1D-and 2D-confined experimental data. We conclude that KCM-involving only a few parameters-provides a predictive description of the thermal and mechanical response in these complex systems with highly non-diffusive behavior and has specific advantages over the traditional effective Fourier model. This work thus represents a significant advance in both experimental and modeling capabilities opening the door to improved thermal management in iterative nanoscale device design, including possible routes to increase clock rates in nanoelectronics by surpassing what has been called the "thermal wall" 40,41 .

II. EXPERIMENTAL DESIGN AND THEORETICAL MODEL
We measure the heat dissipation away from 1D-confined periodic nanoline heat sources on a silicon substrate using dynamic EUV scatterometry similar to that of Refs. 7-9 , but with significantly improved signal-to-noise ratio-by nearly two orders of magnitude.
These improvements allow us to perform new measurements on 2D-confined periodic nanodot heat sources on a silicon substrate, which were fabricated under identical conditions as the nanoline arrays. Our time-resolved measurements use an ultrafast infrared pump laser pulse to rapidly excite thermal heating and expansion in the metallic structures. The resulting thermal and elastic surface deformation is monitored by measuring the change in diffraction efficiency of an ultrashort EUV probe pulse, as depicted in Figure 1 (see Methods). Using this technique, we observe the heat dissipation from nanodot arrays in general geometries without complex fabrication, for the first time, and of nanoline arrays down to 20nm in size (L) and 80nm in spacing (P).
To interpret the experimental data, we implement a mesoscopic model using KCM and a thermoelastic set of equations (see Methods). For heat transport, we use Fourier's law for the metal sources (which is dominated by electrons), and the Guyer and Krumhansl transport equation 16 for the substrate (silicon), which is the material where non-Fourier behavior is expected: where  is the bulk thermal conductivity of substrate, the relaxation time of flux q, and ℓ the nonlocal length-that can be microscopically interpreted as a weighted average phonon mean free path. All these parameters are intrinsic properties of the substrate.

III. RESULTS
We first study effectively isolated heat sources for both 1D-confined (nanolines) and 2D-confined (nanodots) of different sizes and periodicities. Figure 2 compares the experimental results on nanolines and nanodots with theoretical KCM solutions obtained using COMSOL. We compare both inertial solutions, that include elastic waves generated by the impulsive pump laser excitation, and quasi-static solutions without elastic waves to isolate the effects of the heat flow (see Supplementary Section 1). We use ab initio calculations to compute the intrinsic parameters of bulk silicon at T=300K: = 145 W/mK, = 50ps, and ℓ= 176nm 21,37  For the thermal boundary resistance, which is an intrinsic property that depends only on the materials and the fabrication process, we use R1=2.25 nKm 2 /W for all nanostructure geometries (see Methods). The excellent agreement in Figure   2 between experiment and theory demonstrates an advance in modeling as the nanoline thermal decay has already been shown to be highly non-diffusive 8,9 and KCM-which is based on only a few key parameters-accurately predicts the thermal transport and elastic waves in both nanolines and nanodots without any geometry-dependent fit parameters, which is beyond the current capabilities of microscopic descriptions.
In the close-packed situation, ( − ) < 2ℓ, nonlocal effects are expected to yield interaction between heaters, as phonons from a given source are able to reach neighboring sources before scattering. In this case, one does not expect Eq. (1) to be applicable since higher order derivatives should be included in the transport equation 35 . To keep the model as simple as possible, we propose that the effects of these higher order terms can be absorbed into a geometry-defined value ℓ eff , where Eq. (1) is still sufficient to describe the system. We propose the simplest expression that satisfies limiting cases: ℓ eff = ( − )/2 (< ℓ). For this expression, when the period P tends to the linewidth L, ℓ eff → 0. In this limit, the grating tends to a line of infinite linewidth, thus viscous effects should vanish. In the other limit, if ( − ) → 2ℓ we recover ℓ eff → ℓ as constructed. Using this expression for ℓ eff , we compare KCM predictions with experimental results for close-packed nanoline(nanodot) heaters in Figure 1c(d). The model predicts that closely-spaced heat sources cool faster than widelyspaced ones, as uncovered in previous experiments 6,8,9 . We also experimentally demonstrate, for the first time, that this same counter-intuitive behavior observed in nanoline arrays is universal and manifests in nanodot arrays, since the L = 50nm with P = 200nm nanodot signal is relaxed at 800ps while L = 50nm with P = 400nm is not. The excellent agreement between the KCM prediction and experimental results for the close-packed cases shows that KCM can model this behavior with a simple expression for ℓ eff (without fitting) while the other parameters used are the same used in the isolated cases. In summary, both nanoline and nanodot experiments can be predicted by the KCM using the intrinsic value ℓ=176nm when sources are separated a distance larger than 2ℓ (effectively isolated sources), and a geometrydefined effective value when distances are smaller (close-packed sources). This modification of ℓ for a specific situation allows us to retain both the predictive capability and simplicity of the model.
Using our model, we interpret the behavior of the effectively isolated sources from a hydrodynamic viewpoint and compare it to the close-packed sources. For effectively isolated sources, hydrodynamic effects become relevant when linewidth L is on the same scale as the phonon mean free paths ~ℓ; thus, the non-diffusive terms in Eq. (1) reduce the heat flux, compared to Fourier's law, in agreement with experiments [5][6][7][8][9]12,44 . This phenomenon is analogous to a friction that arises from the large gradients in heat flux that impedes heat flow, referred to as a viscous resistance 38 . In other words, when linewidth L is on the same scale as ~ℓ, there is not enough resistive phonon collisions to scatter the heat outward in all directions as diffusion assumes. Instead, the thermal energy is forced straight downward into the substrate over a distance related to ~ℓ before enough resistive phonon collisions occur to dissipate energy in all directions, shown schematically in Figure 1. These hydrodynamic-like friction effects resulting from a lack of resistive collisions have been described in other formalisms albeit with different interpretations. For example, models using a phonon suppression function predict heat flow that is less efficient than Fourier's law when linewidth L is on the same scale as ~ℓ, similar to our hydrodynamic model; however, this phenomenon is interpreted as a reduced number of carriers due to ballistically traveling phonons 14,34 .
Additionally, models incorporating anisotropic behavior of thermal conductivity are parallel to the downward flux forcing predicted by our hydrodynamic model 10  Fourier's law, hydrodynamic effects might be interpreted either as an increase of the thermal boundary resistance 8 or as a reduction of the thermal conductivity near the heater 5,32 . In the effectively isolated case with L=30nm and P=600nm, this region has a size of order ℓ ≃ 200 , while in the close-packed case (P=120nm), it is much smaller and of order ℓ eff ≃ 50 . Therefore, we hypothesize that the interaction of the nearby heat sources in the closepacked scenario reduces the non-local length, decreasing viscous effects, allowing the system to cool more efficiently than with isolated heaters. The microscopic description of this effect is the subject of future work.
To demonstrate the advantages of our hydrodynamic model over the traditional effective Fourier model with a best-fit boundary resistance, we compare the two theoretical predictions to experimental data for the isolated 250nm linewidth case in However, our EUV probe is only sensitive to surface deformations (see Methods), and thus precisely captures two distinct timescales-a signature of non-Fourier heat transport.
A distinct advantage of KCM is that we can gain deeper insight into the two time scales of thermal relaxation by investigating the role played by hydrodynamics. To do this, we analytically solve the thermal equations in the heater and the substrate for the case < ℓ.
In this range, hydrodynamic effects are dominant: the q term in Eq. (1) can be neglected compared to the Laplacian term and the heat flux obeys the (linear) Navier-Stokes equation.
The system of equations obtained is: where 1 is the heater temperature, 2 the average temperature of the substrate at the interface, and 2 ℓ the average substrate temperature in the outer part of the hydrodynamic region, i.e. at a depth of order ℓ below the heater. 1 = ℎ ℎ denotes the heat capacity of the heater per unit surface, with ch and h the specific heat and height of the heater, respectively.

=
(1 + )/ is a heat capacity per unit surface characterizing the substrate, with cs the substrate specific heat, and B a calculated geometric coefficient that for nanolines is 3.0.
1 is the thermal boundary resistance between the metal and the substrate, and 2 = ℓ 2 is a size-dependent thermal resistance due to viscous effects (details in Supplementary Section 2). At short times, 2 ℓ is close to 2 ∞ as heat has not reached this region and Eq. (2) becomes a linear system with a double-exponential decay: with and the characteristic times and weights, which are determined by 1 , 2 , 1 and 2 . Therefore, KCM provides two characteristic times with specific expressions in terms of the physical properties of the system. Equation (2) can be interpreted intuitively as a two-box model as seen in Figure 5.
One box represents the heater, while the other box is a region of order L in the substrate below the heater (referred to as the dam region). The thermal response of the system begins when the heater is filled with thermal energy from the laser pulse. At short times after the laser pulse, the heater releases the energy into the dam region, which retains the energy and rapidly increases in temperature. The initial rate of this energy transfer is dominated by the intrinsic thermal boundary resistance between the heater and substrate. At larger times, when the dam region has equilibrated with the heater, the dissipation of the thermal energy is dominated by the rate of energy transfer out of the dam region into the rest of the substrate. Therefore, the substrate plays two roles in the thermal response of the system: it acts both as an energy reservoir with heat capacity C2 and as a thermal resistance R2. The rate of energy transfer in these later times is controlled by the viscous resistance, i.e. hydrodynamic effects. The thermal relaxation of the heaters can be described by an equivalent circuit ( Figure 5a) and illustrated by a fluid analog (Figure 5c). The predicted temperature evolution of the system as a function of time and position are shown in Figure 5b and Supplementary Section 2.
For small isolated sources, we find simple expressions for the characteristic times, nm, these expressions yield 1 = 50 and 2 = 1050 , thus 2 is an order of magnitude larger than 1 . In this limit, 1 depends on the thermal boundary resistance, while the viscous time scale 2 does not depend on the thermal boundary resistance, but mainly on the nonlocal length ℓ and geometry: Therefore, for small isolated sources, KCM can provide simple analytical expressions for the two different time scales of the heat transfer, each one associated with a different resistive mechanism. This allows accurate experimental validation of the non-local length value for silicon at room temperature (a sensitivity analysis of various KCM parameters is provided in Supplementary Section 2). Additionally, the two-box model Eqs. (2) can also be applied to close-packed experiments by substituting ℓ by ℓ eff ; however, the simple expression of Eq. (4) cannot be used in this case (see Supplementary Section 2).
Although the two-box model has been derived at small sizes, it also characterizes the non-Fourier behavior for all experimental sizes. To validate the intuition provided by the twobox model, we fit a double-exponential decay (Eq. (3)) to each of our experimental measurements, as shown in Figure 6a. We compare the fits of experiments to fits of numerical KCM simulations and the analytical two-box model in Figure 6b-d. We find that the experimental fit results agree well with both KCM numerical and analytical calculations.
Additionally, we confirm the existence of a short time scale ( 1~1 00ps) which is dominated by the intrinsic thermal boundary resistance in Figure 6b and a longer time scale ( 2~1 ns) which is dominated by the hydrodynamic effects in Figure 6c. Figure 6c also displays the splitting of the decay times between effectively isolated and close-packed experiments, i.e.
the increase in dissipation efficiency for close-packed heat sources. In Figure 6d, we plot the weight of the hydrodynamic dominated decay, 2 in Eq.
(3), which shows a transition from a primarily hydrodynamic decay for small heaters, to a decay ruled by the thermal boundary resistance at large sizes. This is expected as large sizes should converge to the Fourier prediction, which contains a single time scale. Therefore, the size-dependent effective boundary resistance extracted by the effective Fourier model in Refs. 8,9 can be re-interpreted as capturing the weighted average of the time-scales ( 1 , 2 ) generated by a size-independent boundary resistance and size-dependent localized hydrodynamic effects.

IV. CONCLUSIONS
In conclusion, we have shown that by adding a hydrodynamic heat transport term, we can explain the thermal transport behavior of nanoscale metal-semiconductor samples with

Comparison between experiment and model predictions:
Since KCM consists of a linear set of partial differential equations, the surface deformation and the predicted diffraction efficiency linearly depends on the amount of energy deposited in the heater by the laser pulse. In the simulations, a uniform energy density of 1W/m 3 with a duration of <2.5ps is introduced in the heater. To compare the model predictions and the experiments, the diffraction efficiency obtained in KCM inertial simulation is scaled by a factor to match the first experimental peak. This is equivalent to scaling the simulated energy density by this factor and this same scaling is used to normalize the quasi-static simulations. This procedure is also applied to the effective Fourier simulations in order to compare Fourier and experiments in Figure 4. Note that a slight correction factor has been added to scaling of the quasi-static simulation of L = 30nm with P = 120nm in Figure 2d due to a small numerical error in the first few picoseconds of this inertial simulation.
The quasi-static solutions are obtained by removing inertial elastic effects, i.e. dynamic equilibrium is imposed during all the simulation (see Supplementary Section 1). These solutions capture the deformations just due to thermal expansion and hence can be used to track the temperature evolution of the system (see ref. 9 ). Note that the initial peak obtained in the quasi-static simulations is not observed in experiments because the system needs a finite time to expand.
Double Exponential Fitting to Experimental and Numerical Data: For Figure 6, we performed double-exponential fits to both the experimental data and the numerical simulations. The quasi-static numerical solutions can be easily fit to a double-exponential using non-linear least squares; however, due to the noise and inertial elastic effects, a doubleexponential function with four free parameters is too unconstrained to reliably fit to the experimental data. Therefore, we constrain the number of free parameters in a single fit while still independently extracting the four parameters of the double-exponential. We achieve this by fitting the data in several different steps. We determine a cut time, tc, to divide the experimental trace into two parts in time to separately fit the two exponentials. We define the cut time as the time when the ratio between the two exponentials is 1% and compute tc using the fit values from the numerical solution. We fit the experimental data for times > to extract the longer decay time exponential; however, the functional form of the decay for large is not purely single exponential. Because diffusive transport occurs far from the heat source at large , the decay has a power law component super-imposed on the exponential. To mitigate the effects of this power law on the extraction of the decay constant, we fit an effective Fourier model-with only two free parameters of effective thermal boundary resistance ( ) and the overall normalization ( 2 )-for > and truncate the fit at roughly two times the expected decay constant. We can then convert to a decay time, To extract the other exponential, we correctly set the overall normalization by accounting for the inertial elastic effects in the experimental data. To do this, we fit the KCM quasi-static simulation to the experimental data with the acoustics waves subtracted for times      Refs. 8,9 and the current study.

Supplementary Information
Thermal transport from 1D-and 2D-confined nanostructures on silicon probed using coherent extreme UV light: General and predictive model yields new understanding

Thermoelastic Model
The model equations can be solved using Finite Elements methods to determine the evolution of the displacement vector u, the temperature T and the heat flux q in the metalic domains and the substrate, which allow comparision with the experimental measurements. The set of equations is the following:

Linear Elastic Equations
To describe the mechanic evolution of the system we use the elastic equation including inertial effects [15] ρ ∂ 2 u where ρ is the density and σ is the stress tensor of the material.
For the heaters we use the bulk nickel stress tensor (isotropic) with linear thermal expansion: where α Ni , K Ni and µ Ni are the nickel coefficient of thermal expansion, the compressibility modulus and the shear modulus, respectively, and T 0 = 300K is the ambient temperature.
For the silicon stress tensor, we also assume linear thermal expansion and we use an anisotropic form to include the effects coming from the stress generated by the nanolines on the substrate top surface [20]. The characterization of this stress tensor has been reported elsewhere [16] and it is necessary to reproduce the frequency of the mechanic oscillations: where α Si is the silicon coefficient of thermal expansion, I is the identity matrix and D is the anisotropic elasticity matrix.

Heat Transport Equations
Heat transport is modeled including the corrections of the phonon hydrodynamic model in the region where heat is carried mainly by phonons i.e. the silicon substrate.

Nickel Nanostructures (Fourier).
Heat transport in metal lines or dots (height h = 11.5nm and width L) is carried by electrons. As their mean free paths are much shorter than the sizes of the nanostructures, nonlocal effects are not expected. Consequently, Fourier's law is valid on these domains, in which we denote the heat flux and the temperature with subindex 1.
We include the linear thermoelastic coupling in the energy conservation equation to model the transfer between thermal and elastic energy: where c Ni is the specific heat of Nickel and p Ni = α Ni K Ni T 0 . The heating energy density Q is 1W/m 3 for time t < 2.5ps and 0 otherwise. The normalized thermal response of the system does not depend on the used value for Q due to the linearity of the model equations. The temporal duration of the pulse is estimated using a two-temperature model [14,2]. Moreover, we confirmed that the exact length of the heat pulse (0-10ps duration) does not significantly affect the results.
For the transport equation we use the Fourier's law where κ Ni is the bulk thermal conductivity of Nickel.
Silicon substrate (Hydrodynamic Transport). We use the energy conservation equation with linear thermoelastic coupling and without the source term, where c Si is the specific heat of silicon and p Si = α Si K Si T 0 .
For the heat transport we use the GKE including non-local and memory effects, where κ Si is the silicon bulk thermal conductivity, is the non local length, τ is the heat flux relaxation time and α is a dimensionless parameter. Equation (7) can be derived from the BTE assuming an averaged phonon mode relaxation time with α = 1/3 [7] and for the collective regime (normal dominant phonon collisions) with α = 2 [8]. In consistency with [21,18,3,4], here we propose the use of an effective form of equation (7) by considering interpolated parameters κ, , τ between the kinetic (resistive dominant phonon collisions) and the collective situations at the reference temperature T 0 , and α = 1/3. Details of the interpolation and the ab initio calculation of these parameters can be found elsewhere [19,18].

Boundary Conditions
Nanostructure Free surfaces Thermal insulation is imposed by setting to zero the normal heat flux component where n is the boundary normal vector.
Interface On one hand we impose continuity of the normal component of the heat flux going through the interface: where n is the interface normal vector pointing towards the semiconductor.
On the other hand, we impose the temperature jump boundary condition that accounts for the nonequilibrium effects introduced by the interface (assuming diffusive reflections): where β, χ are ab initio calculated characteristic lengths. Derivation of this boundary condition and the explicit microscopic expression for these parameters by imposing energy balance restrictions with the use of the specific non-equilibrium distribution function of each domain can be found elsewhere [4]. The tensor χ is diagonal. A lower bound for the thermal resistance value R 1 is also obtained in [4] by assuming phonon diffusive boundary scattering and perfect interface (no reduction of the contact area between nickel and silicon due to fabrication defects), where γ Si 0 = 1 4 c Si v Si and γ Ni 0 = 1 4 c Ni v Ni (with v the average phonon group velocity for each material respectively). The actual thermal resistance value R 1 is unknown due to the lack of knowledge of the interface defects. Therefore we fit a correction factor R 1 /R min 1 = 3.11 to the data obtained from the largest experimentally available gratings L > 750nm, in which hydrodynamic effects are not relevant (i.e. Fourier and KCM predictions coincide). Hence, this correction does not depend on the model used. The second and third terms in the right-hand side of expression (10) are required to properly describe the non-local effects induced by the interface in the substrate heat flux and temperature profiles. In the present work only cause minor corrections for the smallest linewidths L < 50nm. The same correction factor has been applied to the resistance weighting the non-local term γ Si 0 /γ Si = 3.11.

Slip boundary conditions
For silicon, where a second order derivative of the heat flux is included on the thermal transport equation, a boundary condition for the tangential flux q t is required. We use the slip boundary condition both in the silicon interface and free surfaces [1,3] q t = C ∇ q t · n (12) where p is the fraction of specularly reflected phonons in the boundaries (details can be found in [3]). In consistency with (10) here we consider diffusive reflections p = 0, i.e. C = 1. In the silicon free surface, insulation is also imposed Substrate Only a geometry periodically repeated unit is simulated by imposing periodic boundary conditions. The reference temperature T 0 = 300K and null displacement vector u = 0 is fixed in the substrate base.

Nickel
Silicon

Parameter Values
Regarding the interface boundary condition, the calculated conductance bounds are γ Si 0 = 1068MW/m 2 K and γ Ni 0 = 1960MW/m 2 K, and the thermal boundary resistance used is R 1 = 2.25nK m 2 /W. The tensor χ is diagonal with χ xx = −31nm, χ yy = χ zz = −16nm where x denote the normal direction pointing towards the semiconductor. The length β = −21nm.
Regarding the silicon stress tensor, we use an anisotropic model in Voigt notation {yy, xx, zz, xz, yz, zx} with the following rotated tensor The rest of parameter values used in the model for each material can be found in the Table 1.

Thermal Decay Analysis
In this section, we compare the non-equilibrium evolution of the system according to the Fourier model and to the Kinetic Collective Model (KCM) and we derive the two-box model equations from the KCM analytical solutions in the Stokes regime (L < ).

Fourier Model
First consider diffusive heat transport both in the Si substrate and in the Ni heater along with a Kapitza interface boundary condition with resistance R. We use the bulk thermal properties from Table 1 and we denote the thermal diffusivity of Nickel χ Ni = κ Ni /c Ni = 2.2·10 −5 m 2 s −1 and of silicon χ Si = κ Si /c Si = 9·10 −5 m 2 s −1 . For illustration purposes, we discuss this benchmark model considering the spe-cific case of heater height h = 10nm and width L = 20nm, with R = 1nK m 2 /W. In Figure 1 we show the corresponding temperature evolution of the heater obtained with COMSOL Multiphysics. The time scale of the thermal evolution in the heater is extremely fast h 2 /χ Ni = 4.5ps and hence the temperature in the heater is almost uniform within the time scale of the experiment. In the substrate, at time t, diffusion has penetrated a region of size √ tχ Si . We can quantify an effective thermal resistance due to diffusion r(t) = √ tχ Si /κ Si . At early times R > r(t), so the thermal decay is dominated by the interface. At times larger than R 2 κ Si c Si = 232ps, we have r(t) > R thus the thermal decay is dominated by substrate diffusion.
For t < R 2 κ Si c Si : The heat flux in the interface is | q| = ∆T /R, where ∆T is the temperature difference between the heater and the substrate across the interface. Moreover, the heat flux leaving the heater can be estimated as | q| ∼ −c N i h dT dt , and hence ∆T ∼ −c N i hR dT dt . Therefore, the temperature evolution of the heater is an exponential with characteristic time τ F = c N i hR = 40ps.
For t > R 2 κ Si c Si : Substrate diffusion has no characteristic time scale (infinite substrate) and hence the thermal decay follows a power law with an exponent depending on the space dimensionality. The thermal evolution in the region below the heater is instantaneous L 2 /χ Si = 4.4ps and hence we can consider that the heated region is a point and the temperature evolution is 2D (exponent -1).
In summary, the Fourier model predicts an initial exponential thermal decay with characteristic time τ F = c N i hR = 40ps followed by a power law decay with exponent -1. The transition between both decays is estimated to be at time R 2 κ Si c Si = 232ps. Figure 2: Fourier fits to experimental data restricted to t > 500ps for a heater line of 250nm. The TBR fit (R = 19nK m 2 /W) using the intrinsic value for the substrate conductivity underpredicts the temperature for t < 500ps, whereas the conductivity fit (κ = 0.4κ bulk ) using the intrinsic value for the TBR overpredicts it. In both models, the predicted fraction of energy evacuated from the heater is substantially distorted.
As shown in Figure 4 of the main text, the functional form of the thermal decay according to Fourier's law is not consistent with the experimental decay from EUV scatterometry measurements, where a double exponential decay with two distinct characteristic times is observed (see also Figure 6 of the main text). Care must be taken when comparing our results with those from time-domain thermal reflectance (TDTR) [12], which are two completely different techniques.
One key difference is that visible-based probe experiments often need to limit the analysis to the region t > 500ps due to the challenges separating the contributions of out-of-equilibrium electrons and thermal decay. However, our EUV probe does not suffer from this limitation as we measure the surface deformation through diffraction, which is not altered by the presence of nonequilibrium electrons. If we take EUV scatterometry data and exclude from our analysis the measurements for t < 500ps, one of the decay times (∼100ps) reported in this work cannot be observed. Consequently, a diffusive model could fit the experimental data by excluding the initial timescales, as done in previous works like [12]. In Figure 2, we show two different fits using Fourier with an effective TBR and an effective substrate conductivity, respectively, restricted to t > 500ps. Both approaches can reproduce the tail of the decay but fail to reproduce the initial system response. After the first nanosecond, these models are able to reproduce just the last 20% of the signal amplitude but cannot predict the other 80%, i.e. the largest part of the energy dissipation from heaters which is also the most important for applications.

Kinetic Collective Model
We assume now diffusive heat transport in the heater and hydrodynamic heat transport in the substrate using the Kinetic Collective Model as explained in Supplementary Section 1. In figure 3, we show an example of the temperature evolution of the heater for L = 30nm and P = 400nm obtanied with COMSOL. In contrast to the Fourier model, the KCM predicts a slower thermal decay that can be fitted using a double exponential within the time scale of the experiment. We denote as τ 1 , τ 2 , a 1 , a 2 the characteristic times and weights of the first and the second exponentials, respectively. The geometrical dependencies of these coefficients are represented in Figure 6 of the main text. Figure 6

and sensitivity of the KCM parameters
The information condensed in Figure 6 of the main article (reproduced here as Figure 4) is more easily understandable by identifying all the double exponential decay coefficients for an specific sample and analyzing the sensitivity of the decay to changes in the KCM parameters values. Here we consider two examples that represent two different situations. The selected cases are represented in blue and red in Figure 4. The blue circles indicate the coefficients of the double exponential decay for isolated small heaters with size L = 30nm and periodicity P = 400nm (τ 1 = 68ps, τ 2 = 1470ps, a 2 = 1 − a 1 =0.7). The thermal evolution for this sample is repre-sented with blue lines in the two top plots of Figure 5. The influence of modifying the non-local length or the boundary resistance value R 1 is displayed in the left and right plots, respectively. Note that the smaller time τ 1 is determined by R 1 and the larger one τ 2 by , so each decay is associated to a different mechanism. In this case the weight a 2 > 0.5 and the hydrodynamic time is larger than the TBR dominated time (τ 2 > τ 1 ). Therefore, the heater thermal evolution is more sensitive to the value than to the interfacial resistance value R 1 as shown in Figure 5. Left Plots: Thermal decay using the values given by our KCM-ab initio model (blue line) in comparison with the same system using a non-local length multiplied by two (green) and divided by two (red) times. Right Plots: Analogous comparison is provided but changing the used value of the resistance R 1 . Top: Small heater size (L = 30nm and P = 400nm). Bottom: Large heater size (L = 1µm and P = 4µm). Notice that the boundary resistance controls the initial decay τ 1 and the non-Fourier conduction controls the decay at larger times τ 2 . Moreover, the weight of the second exponential with time τ 2 increases by decreasing the line-width.
The red circles in Figure 4 identify the coefficients for large heaters with size L = 1µm (τ 1 = 139ps, τ 2 = 1840ps, a 2 = 0.09). In this sample the situation is the contrary with respect to the previous case. The hydrodynamic weight is small (a 2 < 0.1) and the timescale τ 1 determined by the TBR dominates. In the bottom plots of Figure 5 the decay for this sample is represented in blue, along with the same sensitivity analysis of the KCM parameters. Notice that the decay is mainly influenced by R 1 while the change of does not have much effect. This is expected since for large heater sizes the decay is well described by Fourier's law along a Kapitza interfacial resistance.
Therefore, there is a transition from the thermal relaxation dominated by hydrodynamic effects observed for small heater sizes, to the evolution dominated by the interfacial resistance for large sizes. For intermediate sizes, the two mechanisms are important and a double exponential decay is evident (see Figure 4 of the main text). The experimental validation of this transition is displayed in Figure 6d of the main text.
The Two-Box Model: Analytical Derivation of double exponential thermal decay.
Here we derive the analytical expression for the parameters of the double exponential thermal decay predicted by KCM in the case L < . We denote x the cross-plane direction towards the substrate and y the in-plane direction. The origin of coordinates is the center of the interface.
In KCM, the heat flux in the substrate is described through the hydrodynamic heat transport equation along with the energy conservation equation We neglect the term τ ∂ q ∂t because it does not play a significant role in the present experimental conditions. We also neglect the thermo-elastic coupling in the energy conservation for simplicity. We consider the Stokes regime L < (heat transport dominated by viscosity); then one can neglect the term q in (15) close to the heater since we expect 2 ∇ 2 q ∼ 2 L 2 q. Is worth to note that the heat flux profile saturates after a fast transient and viscosity remains constant during the rest of the experiment. The heat flux profile saturation time can be estimated as 2 /χ Si = 341ps. To illustrate this, in Figure 6 we show the time evolution of the different terms in equation (15) at x = L and y = 0 with L = 20nm.
Now we perform integration in the region dominated by viscous effects x < 2 : κ Si where T 2 is the substrate temperature at the interface x = 0 and T 2 2 is the substrate temperature at x = 2 . At the time scales considered in the experiment T 2 2 is constant and close to the initial temperature T ∞ 2 = T 0 . Moreover, during the experimental time scale, the heat flux and its derivatives at x = 2 are neglectable in front of the heat flux and its derivatives at the interface x = 0.
whereT 2 ,q x are the average temperature and heat flux in the interface, respectively, and is a geometric parameter related with the average heat flux profile in the hydrodynamic region. This parameter saturates because the heat flux profile reach a stationary situation as can be seen in Figure 3. The saturated value can be estimated from the COMSOL simulations. We obtain a constant value for L < in the 1D geometry: B = 3.
Diffusion in the heater region is extremely fast and hence the temperature of the heater T 1 is uniform within the time scale of the experiment. Therefore, using the averaged form of a Kapitza interface boundary conditionq (22) we obtain the following evolution equation where and is the viscous resistance.
Consider now the energy conservation in the heater c Ni By performing volume integration of this equation with using the Kapitza interface boundary condition and the insulation condition for the other boundaries, we obtain an independent evolution equation Notice that, in the derivation of the evolution equations (24,27), a simplified Kapitza interface condition have been used i.e. the hydrodynamic contributions to the interface boundary condition (10) have been neglected for simplicity. The inclusion of the hydrodynamic contributions to the boundary condition (which cause only small deviations for extremely small L < 50nm) can be found in the ending notes of this section.
The heater temperature T 1 can be obtained from the system of partial differential equations (24,27): In the case of isolated heaters (R 1 < R 2 ), the system (24,27) can be simplified and we obtain being C 1 = c Ni h, C 2 = τ S /R 2 the heat capacities of the heater and the dam region, respectively. With these definitions, the system (24,27) is equivalent to (2) in the main text.
In the case considered here (τ 1 < τ 2 i.e. R 1 < R 2 ), If R 1 ∼ R 2 , then the system (24,27) needs to be solved with no approximations. This is the case of interest for small L and P = 4L (close-packed situation) in which we use a reduced non-local length eff (reduced R 2 ). See the ending notes of this section for details.
For illustration, in Figure 7 we show the comparison between the analytical prediction (28) for the heater temperature evolution T 1 (t) with L = 20nm and P = 800nm compared with the Finite Elements calculation. Note I: Solutions of the system of partial differential equations (24,27) Here we revisit the solutions of the system of partial differential equations (24,27). Exponential solutions exp(tw) satisfy being w the roots of the system characteristic polynomial. There are always two real negative roots w 1 = −1/τ 1 and w 2 = −1/τ 2 . In order to find the expressions (29,30) for the characteristic times we assumed R 2 > R 1 and hence we simplified the previous equation to Regarding the weights of the double exponential, in general we have which in the case τ 1 < τ 2 simplify to equations (31,32).
Notice that equations (33,35,36) are the general system of equations for the parameters τ 1 , τ 2 , a 1 , a 2 . Now consider the two geometrical regimes:  and an effective Fourier model for small heat source size and spacing (close-packed). The experimental change in diffraction efficiency for L = 30nm with P = 120nm is shown in black and error is embodied in grey shading. The KCM prediction is shown in blue including the full inertial calculations (solid) and the quasi-static (dashed). An effective Fourier model with a fitted thermal boundary resistance is shown in green. As predicted by the two-box model, we recover a Fourier-like behavior in this extreme close-packed situation.
(ii) Close-packed situation (P = 4L and small L): In this case eff = (P − L)/2 so that R 1 ∼ R 2 . In this case the use of the approximated equation (34) is not acceptable and hence we use the exact solutions of the general equation (33) to compare with experiments. It is easy to show that by reducing L with P = 4L, a 2 goes to 1 and a 1 goes to 0; τ 2 goes to the Fourier decay time τ F = c Ni hR 1 and τ 1 tends to zero. Therefore, we recover Fourier in this limit and we don't expect a clearly observable double exponential decay. In particular, for L < 50nm, P = 4L we expect a single exponential decay as in the Fourier-based description. This limit is consistent with experimental observations as shown in Figure 8.
Note II: Inclusion of the hydrodynamic contributions to the interface boundary condition (10) in the two box model.
Consider the hydrodynamic interface boundary condition (10) instead of the simplified Kapitza condition. We average it along the interface with expressing ∇ · q and ∇ q as quantities proportional to the interface normal heat fluxq x (this is possible because of the saturation of the heat flux profile close to the interface): where are dimensionless parameters estimated from the corresponding saturated value in the COMSOL simulations.
Now we can input this averaged boundary condition to theT 2 evolution equation (22) and the heater energy conservation equation to obtain the two box model system of equations. This system is exactly the same as the one previously obtained (24,27) with a redefined interface boundary resistance R 1 : Therefore, by including the hydrodynamic contributions to the boundary condition we obtain a larger thermal resistance which become explicitly geometry dependent. This correction slightly increase the τ 1 value for small sizes reported in Figure 6 of the main text with respect to the value obtained using the simplified Kapitza boundary resistance. However, for L > 50nm, we have R 1 ∼ R 1 .

Matrix Pencil Method
We utilize the matrix pencil method (MPM) to filter the elastic wave oscillations out of the experimental change in diffraction efficiency signal. Our MPM algorithm begins by forming the time-lag co-variance matrix from an experimental time signal (i.e. the diffraction efficiency signal as a function of pump-probe time delay). We then perform the singular value decomposition of this matrix and create a scree plot of the singular values. From this plot, we separate the relevant time-correlated patterns in the data from the random noise. Next, we compute the complex exponential (e (α+iβ)t ) that best represents each relevant component. The result is a decomposition of the experimental signal into a few complex exponentials and random noise shown in Figure 9. We remove the oscillatory exponentials (exponentials with complex arguments) since they are not of interest in this study, then sum the remaining exponentials and noise. The final result is an experimental time signal without oscillations from elastic waves also shown in Figure 9. In essence, MPM is similar to a robust least squares fitting of multiple complex exponentials to the experimental data. In the presence of noise, MPM can more precisely and accurately extract the damped oscillations than a Fourier transform. A simple Fourier transform (which assumes stationary oscillations) is a sufficient choice for higher signal-to-noise, long lived periodic signals. A simple example is shown in Figure 10. Mathematical details of MPM along with validation and comparison to a Fourier transform can be found in [13,17,5]. An example application can be found in [10] and the numerical algorithm used in this work can be found in [9].   For any nano-grating geometries with no L or P AFM values in Table 2, we use nominal L and P with the average value for h in the numerical KCM simulations.

More Details on Extraction of Decay Times
In Figure 6 of the article, the agreement between numerical predictions and experimental measurements for τ 1 and τ 2 is quite good. However, there exists disagreement between numeric simulations and experiments for the interface decay (τ 1 ) for small linewidth L and for the hydrodynamic decay (τ 2 ) for large linewidth L. For small L, the amplitude of the interface decay (τ 1 ) is much smaller than the amplitude of the hydrodynamic decay (τ 2 ). Thus, the experimental extraction of τ 1 becomes more challenging. For large L, the interface decay dominates over the hydrodynamic decay making it more difficult to experimentally observe the hydrodynamic decay.
When extracting τ 1 from the experimental measurements for small L, we encounter additional challenges since not only is the amplitude of the τ 1 exponential small but also τ 1 is small. When τ 1 is small, it is on a similar scale as the finite response time of the nanostructures and an elastic wave oscillation. Although the pump laser pulse has an ultrashort duration (≈ 25fs), the deformation of the --11.5±1 ---- Table 2: Table of compiled AFM measurements of nano-grating geometries. L is linewidth, P is period, and h is height. Uncertainty is calculated by standard deviation of measurements and AFM image pixel size. We use WSxM to analyze the raw AFM images [11].
system responds on the time scale of 10ps (we measure surface deformation with the probe beam). These inertial effects can inhibit our ability to extract τ 1 It is challenging to extract τ 1 for small L, close-packed (small period P ) nanostructure arrays. In the extreme case where L = 20nm in the close-packed regime, the predicted τ 1 ≈ 10ps while the experimentally observed time required for the nanostructure to reach maximum thermal expansion is ≈ 15ps. Therefore, we cannot extract τ 1 for very small, close-packed nanostructure arrays. For slightly larger L, such as L = 60nm, we can extract τ 1 ; however, the oscillations from elastic waves, the finite response time of the surface deformation, and the chosen experimental time sampling can artificially increase the extracted value of τ 1 . In the example shown in Figure 12, the τ 1 extracted from experimental data is much different than that predicted by quasi-static numerical simulations. However, we also find deviations in the value extracted from inertial numerical simulations compared to the original value predicted by numerical quasi-static simulation, even though the both simulations are in agreement as shown in the article. Since τ 1 is on the similar time scale as one elastic wave oscillation, we cannot extract τ 1 with high accuracy. In Figure 12, sampling the numerical inertial simulation to match the experimental time resolution and then filtering the elastic wave oscillations does not result in a more accurate value.This example suggests that inertial effects, e.g. elastic waves, can inhibit the ability to observe τ 1 in this regime.
It is also challenging to extract τ 1 for small L, effectively isolated (large period P ) nanostructure arrays. Because the change in diffraction efficiency signal is dominated by nanostructure expansion at small L [6], the signal-to-noise ratio is related to the duty cycle, i.e. the ratio L/P . For small L and large P , the signal-to-noise ratio is low and thus extraction of τ 1 is challenging. As shown Figure 12: Example of inertial effects inhibiting observation of τ 1 where linewidth L = 60nm and period P = 240nm. Extracting τ 1 from numerical quasi-static simulations results in a value of ≈ 80ps as seen in upper left. This τ 1 value is much different than what we experimentally observe, shown in the lower left. However, the values of τ 1 extracted from a numerical inertial simulation are also larger than predicted by quasistatic simulations, shown in the top right. This example suggests that inertial effects, e.g. elastic waves, can inhibit the ability to observe τ 1 in this regime. With data sampled at the experimental resolution, filtering the acoustic oscillations does not result in a more accurate value for τ 1 shown in bottom right.
in Figure 13, the τ 1 extracted from experimental data is almost 4 times higher than that predicted by quasi-static numerical simulations. The extraction of τ 1 from the numerical inertial simulation is about 30% higher than the quasi-static prediction. This discrepancy between the two simulations is most likely due to inertial effects as previously described. However, if we sample the inertial simulation at the experimental time resolution, the extracted τ 1 is now three times higher than the quasi-static prediction. Although our experimental capabilities allow us to capture < 10fs dynamics, observing the entire thermal decay (∼ ns) with fs resolution creates massive experimental data sets. Therefore, we chose the experimental time sampling to be on the order of 10ps, even though it is far from technique limits. If we add Gaussian noise (same signal-to-noise ratio as experiment) to the inertial simulation sampled at experimental resolution, our fit algorithm cannot obtain reasonable values for a fitted τ 1 . These results suggests that inertial effects, sampling, and noise inhibit our ability to accurately Figure 13: Example of inertial effects, sampling, and noise inhibiting observation of τ 1 where linewidth L = 20nm and period P = 400nm. The value of τ 1 extracted from experimental data (bottom left) is almost 4 times higher than that predicted by numerical quasi-static simulations (top left). The extracted values of τ 1 from the numerical inertial simulations is about 30% higher than the quasi-static prediction (top right). However, if we sample the inertial simulation at the experimental time resolution and add Gaussian noise, we cannot extract a reasonable value for τ 1 (bottom right). Experimental data is from [6] extract τ 1 for extremely small L. Note, that the error bars in Figure 6 in the article are computed from the standard deviation of the extracted τ 1 values; therefore, the error bars embody the precision of the experimental data and not necessarily the accuracy of the fit procedure. Note that the quasi-static predictions shown in Figures 12 and 13 use different intrinsic parameters than the results in the article.