Spectrally Resolved Specular Reflections of Thermal Phonons from Atomically Rough Surfaces

The reflection of waves from rough surfaces is a fundamental process that plays a role in diverse fields such as optics, acoustics, and seismology. While a quantitative understanding of the reflection process has long been established for many types of waves, the precise manner in which thermal phonons of specific wavelengths reflect from atomically rough surfaces remains unclear owing to limited control over terahertzfrequency phonon generation and detection. Knowledge of these processes is critical for many applications, however, and is particularly important for recent attempts to create novel materials by coherently interfering thermal phonons. Here, we report measurements of a key property for these efforts, the phononwavelength-dependent specularity parameter, which describes the probability of specular reflections of thermal phonons at a surface. Our experiments show evidence of specular surface reflections of terahertz thermal phonons in our samples around room temperature and indicate a sensitivity of these reflections to surface imperfections on the scale of just 2–3 atomic planes. Our work demonstrates a general route to probe the microscopic interactions of thermal phonons with surfaces that are typically inaccessible with traditional experiments.


S1. ADDITIONAL TEM IMAGES
In the main manuscript, we showed the representative transmission electron microscope (TEM) images for the membranes M1 and M2. Figure S1 shows the TEM images for M1 and M2 at a few other locations, demonstrating that the surface of M1 is smoother than M2.
FIG. S1. Cross-sectional TEM images of the surface of membranes M1 (A) and M2 (B). The RMS roughness on the crystalline silicon membrane surface is estimated to be 2.5 ± 0.5Å for membrane M1 and 7 ± 0.5Å for membrane M2. Scale bar: 20Å.

S2. AB-INITIO PROCEDURE TO EXTRACT SPECULARITY PARAMETER
In this section, we describe the ab-initio modeling procedure to extract the phonon specularity parameter from grating period-and temperature-dependent thermal conductivity measurements. As described in the main manuscript, the measured grating period-dependent thermal conductivity of the free-standing membranes at a given temperature in the quasiballistic regime is obtained by solving the phonon Boltzmann transport equation as, where q = 2π/δ is the grating wavevector corresponding to a grating period δ, d is the thickness of the membrane, Λ ζ is the mean free path, p λ is the wavelength-dependent specularity parameter, C ζ is the volumetric mode specific heat and v ζ is the group velocity of a phonon mode ζ ∼ (k, j) denoting a particular phonon wavevector k, polarization j and wavelength λ = 2π/ |k|. The semi-analytical form of S (qΛ ζ , Λ ζ /d, p λ ) is described in ref.
[47] of the main manuscript. To interpret our experiments, the solutions of the BTE (equation S1) requires ab-initio phonon properties of silicon projected onto an equivalent isotropic crystal. We use a Gaussian kernel-based regression method to obtain the isotropic phonon properties for silicon from a complete set of phonon properties in reciprocal space as described below (full Brillouin zone phonon properties provided by Dr. Lucas Lindsay).

A. Gaussian Kernel Regression for Isotropic Phonon Properties
In the equivalent isotropic dataset, all phonon properties are represented only as functions of phonon frequency. In particular, for any property I n (k), where k is the wave vector of the reciprocal space and n is the mode number, we are interested in obtaining: where V is the volume of the supercell and w (k) are the weights on the discrete Brillouin zone grid. To evaluate the sum/integral in eq. S2, we use the following approximation for the δ-function: with Ω being a smearing parameter. Note that if Ω frequency grid spacing ∆ω k , then the δ-function is poorly approximated and the isotropic phonon properties will produce very different macroscopic thermal properties upon evaluation. If Ω ∆ω k , then there will be a number of zeros in the evaluation of I n (ω) resulting in a number of unrealistic jumps in the isotropic phonon properties. To overcome these problems, we follow the adaptive broadening scheme for k-space integration introduced in ref. 1 to choose the parameter Ω. In this technique, the parameter Ω is adaptively chosen according to, where a is a constant on the order of 1 and v g (n, k) is the group velocity of the phonon mode n at a wave vector k. Here, Ω is a function of the phonon mode number n and the wave vector k. Intuitively, when the phonon group velocity is small, a large number of points gets clustered in the dispersion. The adaptive broadening scheme balances this effect by including a fewer number of points in the average (equation S2).

B. Smoothness constraint for specularity parameter
To derive a prior distribution function that reflects the smoothness of the specularity parameter, we define smoothness of the specularity profile p λ (where λ = 1, 2, . . . represents the different phonon modes) as, To admit some uncertainty in our prior knowledge, we add perturbations to the definition of smoothness as, In matrix formulation, we get where, From equation S7, the distribution of Lp is the same as the distribution of the random variable W . Since W is normally distributed, the prior distribution of the specularity profile p λ is given by,

C. Metropolis Hastings Markov Chain Monte Carlo algorithm
The next step is to sample different specularity profiles from the prior probability distribution (equation S9) and compute the posterior probability estimate from Bayes theorem, as described in the main manuscript. In this work, we use a Metropolis Hastings Markov Chain Monte Carlo (MH-MCMC) algorithm to sample trial specularity parameters p tr λ from the prior probability distribution (equation S9). The MH-MCMC algorithm algorithm proceeds as follows: Draw a sample m from the proposal density q (p k , m)
In this method, the parameter γ represents an artificial temperature which controls the level of perturbation on a profile p k of the k th step to obtain a sample m at the (k + 1) th step. Since we only deal with ratios of π (p)'s in this technique, there is no need to know the proportionality constant in equation S9. Moreover, since the samples m are generated from a proposal distribution q(m, p) for an independent and identically distributed normal random variable, standard sampling techniques developed for a multivariate normal distribution can be employed. The information about the coupling between different components of p λ are included in the expression for the acceptance probability α (p k , m).
Once we have these random trial samples (p tr λ ) drawn according the prior distribution function (equation S9), we can compute the residual function N q,T k expt − k BTE (p tr λ ) , σ 2 I by solving the BTE and finally invoke the Bayes theorem to obtain the posterior probability density (π posterior (p tr λ )) for every sample trial specularity profile p tr λ . While plotting the specularity profiles p tr λ as in fig. 4(E) in the main manuscript, we use this posterior probability density (π posterior (p tr λ )) as the intensity of each of the sampled trial specularity profiles p tr λ .

D. Convergence of Bayesian inference
To start the Bayesian inference scheme, we choose an initial guess p 0 that best represents the experimental data points by trial and error. This process reduces the search space for possible specularity profiles that fit the measurements, thereby reducing computational load. To confirm uniqueness of the final solution, we perform convergence studies by varying the artificial temperature parameter γ and the number of samples, N. Figure S2 shows the converged specularity profiles for three different values of γ. When γ is too small, as in fig. S2 (A), the perturbations in the Bayesian sampling are not sufficient to sample the entire specularity parameter solution space. As γ is increased, the sampling procedure spans the entire solution space, resulting in the posterior probability density smoothly dropping down to 0, as shown in fig. S2 (C) by the low intensity specularity profiles towards the edges of the gray region. The specularity profiles outside this converged region of solution space cannot explain the experimental measurements adequately. This convergence behavior demonstrates the uniqueness of the solution obtained using the Bayesian inference technique. We have also tested the convergence of the specularity parameter with respect to the number of samples in the Bayesian inference scheme as shown in fig. S3, and we require about 10 5 accepted specularity samples to span the entire solution space of the specularity profiles, as shown in fig. S3 (C).

S3. RESULTS FOR MEMBRANE M3 (590 NM THICK)
In this section, we summarize the results for the third membrane M3 (590 nm). Figure S4 (A) and (B) show the measured thermal conductivity on sample M3 using the TG technique at 125 K and 400 K. Similar to the membrane M1 (1145 nm) in the main manuscript, the measured thermal conductivity lies far away from the diffuse limit from 100 K up to 400 K for membrane M3. We applied the Bayesian inference technique to thermal conductivity measurements on M3, finding that the specularity profile as a function of phonon wavelength is slightly shifted to the right compared to membrane M1 ( fig.4(E) in the main manuscript), but the wavelengths that are reflected specularly are essentially the same. It is to be noted that the membrane M3 is similar in thickness as the membrane M2 (515 nm) in the main manuscript but with a higher thermal conductivity, indicating the presence of additional specular reflections compared to those in M2.  fig. 3 in the main manuscript, the gray region represents the thermal conductivity bounds enclosing the experimental measurements obtained using Bayesian inference. (C) Specularity parameter versus phonon wavelength obtained from Bayesian inference. The gray region is an intensity plot of the posterior probability distribution for the specularity parameter with the dashed lines indicating a 95% credible interval. For membrane M3, phonons with wavelength less than 35Å are nearly entirely diffusely reflected, while phonons with wavelengths longer than 80Å are reflected nearly completely specularly. Figure S5 (A) explains the challenges of extracting the specularity parameter from thickness-dependent data available in the literature [2][3][4][5][6][7][8][9] . We emphasize that this plot in itself is already an advance compared to the literature because the curves are computed from ab-initio phonon properties and rigorous solutions of the Boltzmann transport equation that were not available for most of these articles. Figure S5 (A) shows that for films with thickness larger than 100 nm, the reported measurements could be explained well by fully diffuse scattering and by using Ziman's model with RMS roughness down to 1 Angstrom. TEM images were mostly not reported and so attributing this variation to sample differences is problematic. As shown in fig. S5 (B), these extremes in the available data leave ambiguity regarding the value of the specularity parameter for different phonon modes.