In this study, we use transient thermal gratings—a non-contact, laser-based thermal metrology technique with intrinsically high accuracy—to investigate room-temperature phonon-mediated thermal transport in two nanoporous holey silicon membranes with limiting dimensions of 120 nm and 250 nm, respectively. We compare the experimental results with *ab initio* calculations of phonon-mediated thermal transport according to the phonon Boltzmann transport equation (BTE) using two different computational techniques. We find that the calculations conducted within the Casimir framework, i.e., based on the BTE with the bulk phonon dispersion and diffuse scattering from surfaces, are in quantitative agreement with the experimental data and thus conclude that this framework is adequate for describing phonon-mediated thermal transport in silicon nanostructures with feature sizes of the order of 100 nm.

## I. INTRODUCTION

Nanoscale thermal transport has become a topic of much recent interest due to the novel transport phenomena that emerge at the micro- and nanoscale^{1,2} and their relevance to technological fields such as microelectronics and thermoelectrics.^{3,4} In semiconductor systems with feature sizes comparable to the phonon mean free path (MFP), size effects can lead to strong reductions in thermal conductivity—making thermal management in microelectronic devices a significant engineering challenge.^{5} In the field of thermoelectrics, nanostructuring has emerged as a key strategy for enhancing the thermoelectric figure of merit $ZT$ by reducing the thermal conductivity without significantly affecting the electronic properties of the material.^{4,6} Traditionally overlooked for thermoelectric applications, silicon has generated recent interest as a material for thermoelectric devices due to the strongly reduced thermal conductivity achievable through nanostructuring.^{7} Experimental results on silicon nanowires have shown thermal conductivity values two orders of magnitude lower than the bulk value and $ZT$ values approaching unity.^{8–10} Two-dimensional “holey silicon” nanostructures—suspended silicon membranes with a periodic array of nanopores—have exhibited thermal conductivity reductions comparable to nanowires^{11–16} while retaining superior relative mechanical strength. Such nanostructures hold great promise for thermoelectric applications due to the wide variety of well-established and scalable fabrication and manufacturing techniques available for silicon.

Thermal transport at the nanoscale differs significantly from macroscopic, diffusive thermal transport. In structures with feature sizes comparable to the MFP of heat-carrying phonons, thermal transport no longer obeys the heat diffusion equation.^{1} One of the earliest attempts to account for non-Fourier phonon-mediated thermal transport in nanostructures was by Casimir,^{17} whose model featured particle-like phonon transport with diffuse scattering at boundaries. Although Casimir’s original model was concerned with thermal transport in rods, the broader formalism of semiclassical particle-like phonon transport with diffuse boundary scattering is expected to be valid for any nanostructure for which $\lambda th\u226a\u2113$ and $\lambda th/2\pi \u2272R$ , where $R$ is the surface roughness, $\lambda th$ is the representative wavelength of heat-carrying phonons, and $\u2113$ is the limiting dimension of the nanostructure. Heat-carrying phonons at room temperature in silicon have single-digit nanometer wavelengths,^{18} which is of the order of lithographically realistic surface roughnesses.^{19,20} Thus, silicon nanostructures with feature sizes $\u2113>10$ nm should be well described by the Casimir formulation of thermal transport—that is, particle-like phonon transport according to the phonon Boltzmann transport equation (BTE) with diffuse scattering from surfaces. Studies comparing experimental results with *ab initio* theory based on the BTE have shown that the Casimir formulation is indeed valid for nanoscale silicon membranes^{21} and silicon nanobeams.^{22} However, there have been highly conflicting reports regarding the validity of the Casimir formulation for thermal transport in nanoporous holey silicon membranes.^{23} Several studies have reported room-temperature effective thermal conductivities reduced by up to an order of magnitude relative to the Casimir formulation predictions for such structures,^{12,13,24} while others have found good agreement between the Casimir formulation and experiment.^{25–27} In some cases, measurements showing deviations from the Casimir formulation predictions for holey silicon nanostructures have been invoked as evidence of “coherent” thermal transport effects at room temperature.^{11,13,28} This notion, however, has been challenged by recent experimental and theoretical works,^{29–31} in which no effect of nanopore lattice disorder on the room-temperature thermal transport was found. It should be noted that reports of “below Casimir” thermal conductivity rely on the measurements of the absolute values of thermal conductivity, which are challenging even for bulk samples.^{32} If far-reaching conclusions are to be drawn from the absolute value of thermal conductivity, then a technique with high absolute accuracy is desirable.

Transient thermal gratings (TTGs) is a non-contact optical technique that measures the time evolution of an impulsively generated sinusoidal temperature profile.^{33,34} The experimental observable is the amplitude of this sinusoidal temperature profile, which decays as heat spreads from the peaks to the nulls of the grating. For a one-dimensional TTG, the amplitude of the thermal profile and, therefore, the intensity of the heterodyned TTG signal is given by

where $\tau \u22611/\alpha q2$, $\alpha $ is the thermal diffusivity, $q\u22612\pi /L$ is the transient grating wavevector, and $L$ is the transient grating period. The only parameter other than $\alpha $ that affects the decay rate is $L$, which can be measured with high accuracy. Thus, the thermal diffusivity can be determined to high accuracy from the decay rate of the TTG signal. Furthermore, TTG’s non-contact nature reduces additional sources of error due to the absence of any interfaces with metrological structures.

In this paper, two 250 nm-thick holey silicon membrane nanostructures are investigated with the TTG technique. The experimental results from TTG measurements are compared to the results of two *ab initio* numerical Boltzmann transport techniques: the OpenBTE computational framework developed by Romano and Grossman^{35} and the energy-based deviational Monte Carlo BTE simulation technique developed by Péraud and Hadjiconstantinou.^{36,37} Quantitative agreement between numerical calculations and experiment is found for both the unpatterned silicon membrane and holey silicon structures, confirming the validity of the Casimir formulation for room temperature heat transport in silicon nanostructures with feature sizes on the order of 100 nm.

## II. EXPERIMENTAL

### A. Sample fabrication

The holey silicon structures were fabricated using electron beam lithography (EBL) and reactive ion etching (RIE) of a 250 nm-thick freestanding silicon membrane $3.2\xd73.2$ $\mu $m window area, obtained from Norcada Inc.^{38} Each of the two structures was a 100 $\mu $m-diameter region of the freestanding membrane patterned with a square lattice of nanopores. SEM micrographs of the regions are shown in Fig. 1. “Region A” had a pitch size (nanopore periodicity) of 400 nm and a nanopore diameter of 280 nm, and “region B” had a pitch size of 500 nm and a nanopore diameter of 250 nm.

### B. Transient thermal grating (TTG) measurements

As shown in Fig. 2(a), two “pump” laser pulses are crossed at the sample, where optical interference and subsequent absorption lead to the establishment of a transient sinusoidal temperature profile with spatial period $L=\lambda /2sin\u2061(\theta /2)$, where $\lambda $ is the pump wavelength and $\theta $ is the crossing angle for the two pump beams. Through the temperature dependence of the material’s complex index of refraction $n~\u2261n+ik$—where $n$ is the real index of refraction and $k$ is the absorption coefficient—this sinusoidal temperature profile is accompanied by a spatially sinusoidal modulation in $n~$ as well. A quasi-continuous “probe” beam then impinges on the sample, diffracting from this transient optical grating. As the amplitude of the temperature grating diminishes due to heat transport from the peaks to the troughs, the amplitude of the grating in $n~$—and thus the amplitude of the diffracted signal—diminishes accordingly. In this way, the time dependence of the diffracted signal can be directly related to the thermal diffusivity according to Eq. (1). TTG measures the thermal transport dynamics over a length scale set by the period of the transient grating, which can be tuned by changing the crossing angle of the pump beams. Further details regarding this technique can be found in Ref. 34.

The pump beams were derived from a 515 nm source with a 60 ps pulse duration and 1 kHz repetition rate, and the probe beam was derived from a continuous-wave 532 nm source. A “reference” beam was derived from the same source as the probe beam, and the relative phase between the two was controlled by tilting a highly parallel optical flat through which the probe beam passes to achieve heterodyne detection.^{39} At the sample, the probe beam diffracts from the transient grating and becomes superposed with the transmitted reference beam, and the combined heterodyned signal is collected by a fast photodiode detector and recorded on an oscilloscope. The $1/e2$-intensity radius of the pump and probe beams were 100 $\mu $m and 40 $\mu $m, respectively. While the pump spot size is commensurate with the patterned regions, the probe spot size is much smaller. Thus, although our pump may be exciting a grating pattern that extends somewhat outside of the patterned region, our experiment is only sensitive to the transport dynamics within the region bounded by the much smaller probe spot. The pump pulse energies ranged from 170 to 340 nJ, and the instantaneous power of the probe beam at the sample ranged from 0.8 to 1.6 mW. The probe beam was shuttered by an electro-optic modulator with a duty cycle of 5% to prevent steady-state heating of the sample.

The thickness of the membrane was smaller than the optical penetration depth of silicon for the wavelengths involved in the measurements, which permitted measurements in the transmission geometry as shown in Fig. 2(a). The raw TTG data obtained from the two holey regions and the unpatterned silicon membrane at a grating period of 4.25 $\mu $m are shown in Fig. 2(b). Measurements were performed under medium vacuum at a pressure of 1 mbar. The maximum amplitude of the temperature grating was determined to have an upper bound of 35 K. Upper bounds on the average heating of the sample due to the pump and probe beams were determined to be 20 K each.

The TTG signal for a one-dimensional thermal grating exhibiting diffusive thermal transport is given by Eq. (1). $\tau $—or equivalently $\alpha $—is the only free parameter required to model the normalized TTG signal. In addition, the heterodyne detection scheme further yields a gain in signal-to-noise by a factor of the reference field amplitude, which can be increased arbitrarily up to the saturation threshold of the photodetector. The low-dimensionality of the dynamical parameter space and the fact that neither precise knowledge of the magnitude of the temperature variation nor of the magnitude of the heat flux is required in the analysis of the data allow for the determination of the thermal diffusivity with high absolute accuracy. Further discussion regarding the accuracy of transmission-geometry TTG experiments on nanomembranes can be found in Ref. 34. The traces were truncated such that fitting began 5 ns after pump incidence to ensure that the fitted region corresponds only to thermal transport signal without any potential contribution from the fast electronic response shown in the inset of Fig. 2(b). The acquired fits are plotted alongside the raw TTG data in Fig. 2(b). Figure 2(c) shows the measured thermal diffusivity values obtained according to Eq. (1) as a function of TTG period for each of the three regions measured. Each raw TTG trace consisted of 50 000 individual measurements. The statistical error of the measurement was determined by partitioning the data into subsets of 10 000 measurements, fitting each subset to Eq. (1) and taking the standard error of the mean of the resulting distribution of $\tau $ values. In addition to the statistical error of the measurement, the systematic error due to laser heating effects was also considered. The effects of laser heating were determined by performing each measurement three times—once at a baseline set of pump and probe powers, and two additional times at which the pump and probe powers, respectively, were doubled. Linearly extrapolating the measured values of $\tau $ to zero pump and probe laser power allows us to determine the systematic error due to laser heating, which was then added to the appropriate side of the errorbars for each point to account for this systematic heating effect. We note that the upper bounds on laser heating reported above are non-negligible relative to room temperature. However, since the effect of laser heating is experimentally quantified in our error analysis, we can still compare our experimental results with calculations that use room-temperature material properties. Despite the somewhat high upper bounds on laser heating, we nevertheless note that the effect of laser heating on the experimentally determined values of $\alpha $ was generally found to be less than 10%.

For grating periods from 4.25 to 7.5 $\mu $m, we find that the experimental values of thermal diffusivity are independent of $L$ for both the unpatterned membrane and the holey membranes, consistent with preliminary TTG results on holey silicon structures.^{34} The exponential form of the TTG data and the invariance of thermal diffusivity as a function of grating period indicates that the transport kinetics are “effectively diffusive” over the TTG experimental length scales, albeit with “effective” thermal diffusivity values $\alpha eff$ reduced relative to the bulk because of the non-Fourier size effect due to nanostructuring.

It should be noted that occasionally an additional transient with a characteristic timescale much longer than the acquisition timescale (i.e., approximately a constant offset from the pre-pump baseline) was observed in some of the obtained TTG traces. However, we determined that the presence of this contribution to the signal (which is roughly on the timescale that would correspond to thermal diffusion out of the pump spot) was not associated with any change in the $\alpha eff$ value that was calculated from the time constant of the exponentially decaying contribution to the signal observed on the 10s–100s of ns timescale (which we took to be the true TTG signal) that remained after subtracting out this approximately constant offset. This issue is more thoroughly addressed in the supplementary material.

Experimental values of the effective thermal conductivity $\kappa eff$ were calculated from the data in Fig. 2(c) according to

where $\varphi $ is the void fraction of the holey silicon membrane and $cSi$ is the bulk volumetric specific heat of silicon. The resulting experimental values of $\kappa eff$ are shown in Fig. 3, where the effective thermal conductivity values are plotted against the neck width $\u2113n$ (i.e., the difference between the pitch size and the nanopore diameter).

## III. COMPARISON TO FIRST-PRINCIPLES NUMERICAL CALCULATIONS

Numerical calculations of the thermal transport through the membranes were performed according to the linearized isotropic phonon Boltzmann transport equation (BTE) under the single-mode relaxation time approximation (RTA), which is given by

where $fkp(r,t)$ is the occupation function for a mode traveling with wavevector $k$ and polarization $p$, $vkp$ is the (isotropic) group velocity, $f0(\u210f\omega ,TL(r,t))$ is the Bose–Einstein distribution, $TL(r,t)$ is the local temperature field defined such that energy is locally conserved, $\u210f\omega $ is the phonon energy, and $\tau kp$ is the (isotropic) single-mode relaxation time (where $k\u2261|k|$).

The simulation domain is one pore-centered unit cell of the nanopore lattice with the cylindrical axis of the pore chosen to be oriented along $z^$. Periodic boundary conditions are applied along both the $x$- and $y$-axes. The phonon group velocities and relaxation times were determined, respectively, from the harmonic and anharmonic force constants, which were obtained from density functional theory (DFT) calculations using the temperature dependent effective potential (TDEP) method.^{40} Naturally occurring isotope disorder was taken into account. Details on the DFT calculations can be found in the supplementary material.

The OpenBTE computational technique of Romano and Grossman^{35} and the energy-based deviational Monte Carlo BTE (MC-BTE) technique of Péraud and Hadjiconstantinou^{36,37} were both used for *ab initio* calculations of $\kappa eff$ for both the holey and unpatterned membranes.

where $s^(\Omega )$ is the unit vector for the propagation direction $\Omega $, $T(r,\Omega ,\Lambda )$ is the “effective temperature” of phonons with MFP $\Lambda $ traveling in direction $\Omega $ (i.e., the sum of their energy densities divided by $cSi$), $K(\Lambda )$ is the bulk MFP distribution (i.e., the derivative of the thermal conductivity accumulation function with respect to $\Lambda $), and $\u27e8x(\Omega )\u27e9\u2261(1/4\pi )\u222b4\pi x\Omega d\Omega $ is the angular average over all propagation directions. Equation (4) is derived by imposing steady-state conditions on Eq. (3), as well as assuming that $\delta T(r)\u2261TL(r)\u2212T0$ (where $T0$ is the reference temperature, which in this study was 300 K) is small such that $f0[\u210f\omega ,TL(r)]$ in Eq. (3) can be expanded to first order in $\delta T(r)$ and any temperature dependence of material properties can be neglected. The advantage of this approach is that the only input required to solve Eq. (4) is the MFP distribution $K(\Lambda )$.

In OpenBTE, a difference of temperature $\Delta T$ is applied at the two opposing faces of the unit cell along the $x$ axis, and the first guess for $TL$ was given by the standard diffusive equation. Diffuse-scattering boundary conditions were imposed on all surfaces of the computational domain. Diffuse scattering at boundaries was modeled in such a way that phonons of a given value of $\Lambda $ were diffusely emitted (i.e., such that the distribution of outgoing phonons next to the surface is isotropic) from the surface with a total energy equal to the total energy of all incident particles with the same value of $\Lambda $. To overcome numerical instability due to small-MFP phonons, OpenBTE switches to a modified Fourier’s law to compute the diffusive component to heat transport^{41} for such modes. Equation (4) was solved by the finite-volume method while a Delaunay mesh was generated for space discretization.^{42,43}

The deviational energy-based MC-BTE technique^{36,37} was also used to calculate $\kappa eff$ for the nanostructures investigated. This technique achieves low statistical variance compared to other Monte Carlo techniques by only simulating the trajectories of “deviational” particles that describe the excess/deficit thermal energy in a given mode relative to equilibrium, and achieves high computational efficiency by performing the calculation in an energy-based BTE formulation that lends itself naturally to energy conservation. The diffuse boundary scattering condition was modeled in the same fashion as in the OpenBTE method described above—namely, deviational particles with a given MFP were diffusely emitted from the surface with a total energy equal to the total energy of all incident particles with the same MFP. Unlike the OpenBTE technique as described in Ref. 35, the MC-BTE solver applies a constant temperature gradient throughout the simulation domain. To assess the impact of this discrepancy in the applied perturbation, we compute $\kappa eff$ with OpenBTE using both approaches on test aligned structures, finding negligible differences in effective thermal conductivity values.

For both computational techniques, the conductance of one unit cell was calculated by dividing the total heat flux through one end of the simulation domain by $\Delta T$. The effective thermal conductivity $\kappa eff$ was then obtained by dividing this conductance value by the rectangular cross-sectional area of the unit cell normal to $x^$ and by multiplying by the unit cell length along $x^$.^{44} Our results for both computational techniques are shown in Fig. 3 for comparison to the experimental TTG results for both holey silicon structures and the unpatterned membrane. A discussion of the uncertainties associated with the computational methods is provided in the supplementary material.

In a previous paper,^{34} preliminary TTG results investigating thermal transport in a similar holey silicon membrane were compared to the values of $\kappa eff$ obtained from *ab initio* MC-BTE simulations. Agreement between experiment and theory was found to be within $\u223c20%$. However, the sample in that study was patterned over the entirety of the suspended membrane, and as such comparison with an unpatterned region to ensure the intrinsic quality of the silicon membrane was not possible. There, thus, remained an ambiguity in the previous study as to whether the discrepancy between theory and experiment was due to deviations from the Casimir formulation or due simply to material quality effects. Our computational results for $\kappa eff$ are shown alongside our experimental results in Fig. 3 as well as the values of $\kappa eff$ for the structures calculated using the Fourier law with the bulk silicon thermal conductivity value of 143 W/m K. We see that the size effect associated with the thickness of the unpatterned membrane alone reduces $\kappa eff$ by nearly a factor of two relative to the value obtained from the Fourier law (which is simply the value for bulk silicon in the case of the unpatterned membrane), in good agreement with previous measurements on nanoscale silicon membranes.^{21} A further reduction of $\kappa eff$ is observed due to the nanopore superlattice patterning, resulting in a reduction of $\kappa eff$ by a factor of 3 relative to the Fourier law prediction for region A and a near order of magnitude reduction in $\kappa eff$ for region A relative to bulk silicon. We observe very good agreement between the experimental and calculated results (with the two computational techniques being in excellent agreement with each other), which indicates that the Casimir formulation is valid for nanostructures of this kind.

## IV. CONCLUSIONS

We have used the non-contact optical TTG method to investigate thermal transport in holey regions and an unpatterned region of the same silicon membrane. We observe effectively diffusive transport at grating periods larger than 4 $\mu $m and a reduction in effective thermal conductivity by nearly an order of magnitude relative to the bulk value. Two *ab initio* numerical techniques simulating transport according to the semiclassical phonon Boltzmann transport equation yielded excellent agreement with the measurements. Our results indicate that the Casimir framework of semiclassical particle-like phonon-mediated thermal transport with diffuse boundary scattering is adequate for describing thermal transport in holey silicon structures with limiting dimensions of $\u223c100$ nm.

## SUPPLEMENTARY MATERIAL

See the supplementary material for discussion regarding the long-time contribution to the obtained TTG signals that appears in some of the measurements, details regarding the density functional theory calculations, and discussion of computational uncertainties.

## ACKNOWLEDGMENTS

We would like to thank Charles Shi and Jonas Rajagopal for assisting with the TTG measurements. R.A.D. and K.A.N. acknowledge support from the NSF EFRI 2-DARE (Grant No. EFMA-1542864). A.A.M. and K.A.N. acknowledge support by the Solid State Solar-Thermal Energy Conversion Center (S3TEC), an Energy Frontier Research Center funded by the U.S. Department of Energy (DOE), Office of Science, Office of Basic Science, under Award No. DE-SC0001299. The ICN2 is funded by the CERCA program/Generalitat de Catalunya and supported by the Severo Ochoa Centres of Excellence program, funded by the Spanish Research Agency (AEI, Grant No. SEV-2017-0706). M.S. and C.M.S.T. acknowledge support from the Spanish National Project SIP (No. PGC2018-101743-B-100) and from AGAUR (Grant No. 2017SGR806). O.H. gratefully acknowledges financial support from the VINN Excellence Center for Functional Nanoscale Materials (FunMat-2) (Grant No. 2016-05156) and the Knut and Alice Wallenberg Foundation through Wallenberg Scholar (Grant No. 2018.0194). This work has been partly supported by the La Caixa Foundation MISTI Global Seed Fund program (No. LCF/PR/MIT18/11830008).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.