Unifying first principle theoretical predictions and experimental measurements of size effects on thermal transport in SiGe alloys

TL;DR

This study combines first principles calculations and experimental measurements to accurately predict size effects on thermal conductivity in SiGe alloys, demonstrating a reliable framework for thermal transport analysis.

Contribution

It introduces a unified approach linking DFT-based phonon properties with experimental TTG measurements to predict size effects in SiGe alloys.

Findings

Theoretical predictions match experimental data within 25% reduction in thermal conductivity.

The variational solution to the phonon BTE accurately models size effects.

Agreement validated against Monte Carlo simulations.

Abstract

In this work, we demonstrate the correspondence between first principle calculations and experimental measurements of size effects on thermal transport in SiGe alloys. Transient thermal grating (TTG) is used to measure the effective thermal conductivity. The virtual crystal approximation under the density functional theory (DFT) framework combined with impurity scattering is used to determine the phonon properties for the exact alloy composition of the measured samples. With these properties, classical size effects are calculated for the experimental geometry of reflection mode TTG using the recently-developed variational solution to the phonon Boltzmann transport equation (BTE), which is verified against established Monte Carlo simulations. We find agreement between theoretical predictions and experimental measurements in the reduction of thermal conductivity (as much as $\sim$ 25\% of the bulk value) across grating periods spanning one order of magnitude. This work provides a framework for the tabletop study of size effects on thermal transport.

Equations

$$A_{VCA}=x{A}^{Si}+(1-x)A^{Ge}$$

$$g=\sum _{i\in Si,Ge}{f}_i(1-\frac{m_i}{m_{vca}}{)}^2$$

$$Q=\delta (t){\mathrm{e}}^{iqx}{U}_0\beta {\mathrm{e}}^{-\beta z}$$

$$\frac{\partial h}{\partial t}=-\alpha q^{2}h+\alpha \frac{{\partial }^{2}h}{\partial z^2}+\frac{\beta U_0}{C{T}_0}{\mathrm{e}}^{-\beta z}\delta (t)$$

$$\begin{array}{rl}h(z,t=0^{-})& =0\\ \frac{\partial h}{\partial z}{|}_{z=0}& =0\\ h(z\to \infty ,t)& =0\end{array}$$

$$\hat{h}(z,s)=\frac{\frac{\beta U_0}{C{T}_0}}{s+\alpha (q^2-{\beta }^2)}\left({\mathrm{e}}^{-\beta z}-\frac{\beta }{\sqrt{q^2+\frac{s}{a}}}{\mathrm{e}}^{-z\sqrt{q^2+\frac{s}{a}}}\right)$$

$$h(z,t)=\frac{\beta U_0}{2C{T}_0}{\mathrm{e}}^{-\alpha t(q^2-{\beta }^2)}\left({\mathrm{e}}^{\beta z}\mathrm{erfc}(\beta \sqrt{\alpha t}+\frac{z}{2\sqrt{\alpha t}})+{\mathrm{e}}^{-\beta z}\mathrm{erfc}(\beta \sqrt{\alpha t}-\frac{z}{2\sqrt{\alpha t}})\right)$$

$$h(z=0,t)=\frac{\beta U_0}{C{T}_0}{\mathrm{e}}^{-\alpha t(q^2-{\beta }^2)}\mathrm{erfc}(\beta \sqrt{\alpha t})$$

$$\frac{\partial g_{\omega }}{\partial t}+{\mathbf{v}}_{\omega }\cdot \nabla g_{\omega }=\frac{g_0-g_{\omega }}{{\tau }_{\omega }}+\frac{Q_{\omega }}{4\pi }$$

$$\frac{\partial {\hat{\tilde{g}}}_{\omega }}{\partial z}+{\hat{\tilde{g}}}_{\omega }\frac{1+s{\tau }_{\omega }+i{\eta }_{\omega }{\mu }_x}{{\Lambda }_{\omega }{\mu }_z}=\frac{{\hat{\tilde{g}}}_0+{\tau }_{\omega }\frac{\widetilde{Q_{\omega }}}{4\pi }}{{\Lambda }_{\omega }{\mu }_z}$$

$${\hat{\tilde{g}}}_{\omega }(z,s,{\mu }_x,{\mu }_z)={\mathrm{e}}^{-Vz}{\hat{\tilde{g}}}_{\omega }(z=0,s,{\mu }_x,{\mu }_z)+{\int }_0^{z}d{z}^{\prime }{\mathrm{e}}^{-V(z-z^{\prime })}\frac{{\hat{\tilde{g}}}_0(z^{\prime },s)+{\tau }_{\omega }\frac{\widetilde{Q_{\omega }}}{4\pi }}{{\Lambda }_{\omega }{\mu }_z}$$

$$\begin{array}{rl}{\hat{\tilde{g}}}_{\omega }(z=L,s,{\mu }_x,{\mu }_z<0)& =0\\ {\hat{\tilde{g}}}_{\omega }(z=0,s,{\mu }_x,{\mu }_z>0)& =\sigma \end{array}$$

$$\begin{array}{rl}{\hat{\tilde{g}}}_{\omega }(z,s,{\mu }_x,{\mu }_z)& =-\Theta (-{\mu }_z){\int }_z^{\infty }{\mathrm{e}}^{-V(z-z^{\prime })}\frac{{\hat{\tilde{g}}}_0(z^{\prime },s)+{\tau }_{\omega }\frac{\widetilde{Q_{\omega }}}{4\pi }}{{\Lambda }_{\omega }{\mu }_z}\\ & +\Theta ({\mu }_z)\left({\int }_0^z{\mathrm{e}}^{-V(z-z^{\prime })}\frac{{\hat{\tilde{g}}}_0(z^{\prime },s)+{\tau }_{\omega }\frac{\widetilde{Q_{\omega }}}{4\pi }}{{\Lambda }_{\omega }{\mu }_z}+{\int }_0^{\infty }2{\mathrm{e}}^{-V(z)}{F}_2(z^{\prime })\frac{{\hat{\tilde{g}}}_0(z^{\prime },s)+{\tau }_{\omega }\frac{\widetilde{Q_{\omega }}}{4\pi }}{{\Lambda }_{\omega }}\right)\end{array}$$

$$F_n(z)=\frac{1}{2\pi }\int d\Omega \Theta ({\mu }_z){\mu }_z^{n-2}(e{)}^{-Vz}$$

$$4\pi \int d\omega \frac{1}{{\tau }_{\omega }}{\hat{\tilde{g}}}_0(z,s)=\int d\omega \frac{1}{{\tau }_{\omega }}\int d\Omega {\hat{\tilde{g}}}_{\omega }(z,s,{\mu }_x,{\mu }_z)$$

$$\hat{h}(z,s)\int d\omega \frac{C_{\omega }}{{\tau }_{\omega }}=\int d\omega \frac{C_{\omega }}{2{\Lambda }_{\omega }{\tau }_{\omega }}{\int }_0^{\infty }d{z}^{\prime }\left(\hat{h}(z^{\prime },s)+\frac{{\tau }_{\omega }\widetilde{Q_{\omega }}(z^{\prime })}{C_{\omega }{T}_0}\right)\left(F_1(|z-z^{\prime }|)+2F_2(z)F_2(z^{\prime })\right)$$

$${\tilde{Q}}_{\omega }(z)=\frac{C_{\omega }}{C}{U}_0\beta {\mathrm{e}}^{-\beta z}$$

$$U_0\frac{\lambda }{\pi }=2i{\int }_0^{\infty }dz{\int }_0^{\infty }dt{\tilde{q}}_x(z,t)$$

$${\hat{\tilde{q}}}_x(z,s)=\frac{T_0}{2}\int d\omega \frac{C_{\omega }{v}_{\omega }}{{\Lambda }_{\omega }}{\int }_0^{\infty }d{z}^{\prime }\left(\hat{h}(z^{\prime },s)+\frac{{\tau }_{\omega }\widetilde{Q_{\omega }}(z^{\prime })}{C_{\omega }{T}_0}\right)\left(G_1(|z-z^{\prime }|)+2G_2(z)F_2(z^{\prime })\right)$$

$$G_n(z)=\frac{1}{2\pi }\int d\Omega \Theta ({\mu }_z){\mu }_z^{n-2}{\mu }_x(e{)}^{-Vz}$$

$$k=\frac{\frac{1}{3}\int d\omega C_{\omega }{v}_{\omega }{\Lambda }_{\omega }f({\eta }_{\omega },{\mathrm{Kn}}_{\omega })}{\frac{1}{C}\int d\omega C_{\omega }g({\eta }_{\omega },{\mathrm{Kn}}_{\omega })}$$

$$\begin{array}{rl}f({\eta }_{\omega },{\mathrm{Kn}}_{\omega })& =\frac{3}{{\eta }_{\omega }^2}\left(1-\frac{1}{{\eta }_{\omega }}\arctan ({\eta }_{\omega })+\frac{{\eta }_{\omega }^2\Psi ({\eta }_{\omega },{\mathrm{Kn}}_{\omega })-{\mathrm{Kn}}_{\omega }^2\Psi ({\eta }_{\omega },{\mathrm{Kn}}_{\omega })}{{\eta }_{\omega }^2-{\mathrm{Kn}}_{\omega }^2}\right)\\ g({\eta }_{\omega },{\mathrm{Kn}}_{\omega })& =\frac{1}{{\eta }_{\omega }}\arctan ({\eta }_{\omega })+\Psi ({\eta }_{\omega },{\mathrm{Kn}}_{\omega })\end{array}$$

$$\begin{array}{rl}\Psi (x,z)& =\frac{1}{2}{\psi }_1(x,z)-\frac{1}{1+\sqrt{1+x^2}}{\psi }_0(x,z)\\ {\psi }_n(x,z)& =\frac{1}{2\pi }\int d\Omega \Theta ({\mu }_z)\frac{z{\mu }_z}{(1+ix{\mu }_x{)}^n(1+z{\mu }_z+ix{\mu }_x)}\end{array}$$

$$\begin{array}{rl}f({\eta }_{\omega },{\mathrm{Kn}}_{\omega }\to \infty )& =\frac{3}{2{\eta }_{\omega }^2}\left(1-\frac{1}{{\eta }_{\omega }}\arctan ({\eta }_{\omega })\right)-\frac{1}{{\eta }^3(1+\sqrt{1+{\eta }^2})}\left((1+{\eta }_{\omega }{)}^{\frac{3}{2}}-\frac{3}{2}{\eta }_{\omega }^2-{\eta }_{\omega }^3-1\right)\\ g({\eta }_{\omega },{\mathrm{Kn}}_{\omega }\to \infty )& =\frac{1}{2{\eta }_{\omega }}\arctan ({\eta }_{\omega })+\frac{1}{1+\sqrt{1+{\eta }^2}}\end{array}$$

Full text

I Introduction

Silicon-Germanium (SiGe) alloys are the canonical example for the study of thermal transport in a mass-disordered, yet crystalline system as evidenced by the plethora of work, dating back to the original work by Stohr Stohr and Klemm (1939) and Abeles Abeles et al. (1962); Abeles (1963), where it was noted that the the mass-disorder scatters short-wavelength phonons consequently shifting the dominant contribution to thermal conductivity to long wavelength phonons.

On the theoretical side, earlier works relied on empirical models. Skye and Schelling used molecular dynamics to study the relative contribution between mass and bond disorder, finding larger resistivity than experiment Skye and Schelling (2008). Bera et al. estimated the mean free paths (MFP) of phonons in SiGe using a model based on the expected scalings of the phonon lifetimes Bera et al. (2010). Significant progress took place when Garg et al. demonstrated the viability of first principle approaches to estimate the bulk thermal conductivity of SiGe Garg et al. (2011). Recently, Iskandar modeled thin SiGe films to include the effect of boundary modes Iskandar et al. (2015).

On the experimental side, Koh et al. reported a modulation frequency dependent estimate of thermal conductivity under the Fourier model of the experimental geometry of time domain thermoreflectance (TDTR) Koh and Cahill (2007). The authors proceeded to argue that the frequency dependence corresponds to suppression of phonons with MFP greater than the thermal penetration depth. This result led to a series of theoretical explanations. da Cruz et al. presented a framework to divide thermal transport into harmonic and anharmonic channels da Cruz et al. (2012). Vermeersch et al. explained the divergence from bulk using truncated Levy walks Vermeersch et al. (2015a, b). Recently, Hua et al. have found that a microscopic model of the interface between the aluminum transducer and the substrate can be used to explain away this frequency dependence and recover a bulk value of SiGe Hua et al. (2016). Wilson et al. reported decreasing thermal conductivity and increasing interface conductance with increasing modulation frequency and argue that the interplay between the contributions of long wavevector phonons and interfacial scattering is responsible for this observation Wilson and Cahill (2014).

The objective of this work is to provide a framework upon which size effects on thermal transport can be probed and understood with tabletop experiments. This is accomplished using a bottom-up theoretical approach combined with a simple experimental geometry that is not obfuscated by the interface present in TDTR; the transient thermal grating (TTG) technique. In doing so, we are able to unify the pictures obtained from the macroscopic observables of experiment to the microscopic properties from theory.

The structure of the paper is as follows. In Section II we present the phonon properties obtained using density functional theory. In Section II.2, the variational solution to the phonon BTE for the TTG experimental geometry is developed. In Section  III, results obtained from TTG are presented and compared with our BTE-based predictions. Finally, we close with a look towards future work in Section IV.

II Theory

II.1 First Principle Calculations

We follow the general procedure established by Broido Broido et al. (2005, 2007) and Esfarjani Esfarjani et al. (2011), to obtain the phonon properties for SiGe. While the details can be found in these works, an outline of the procedure is included for the sake of completeness.

For a non-alloy system, the harmonic phonon properties (dispersion, heat capacity, group velocity) are obtained using density-functional perturbation theory (DFPT). The underlying premise is to treat the mechanical displacement corresponding the wavevector of a phonon as a linear perturbation to the electronic Hamiltonian, from which atomic forces can be calculated under the self-consistent criteria of DFT. These forces are then converted into harmonic force constants and used to construct the dynamical matrix for the perturbing wavevector, which can then be diagionalized to obtain the corresponding frequencies.

The anharmonic properties, namely lifetimes (but also frequency shifts) can be obtained by extending the perturbation to higher orders. An alternative approach is to construct a symmetry-reduced set (based on the space group of the lattice) of atomic displacements in a supercell, where each member of the set undergoes a standard DFT self-consistent calculation, each yielding the force field for the configuration. With this set of force fields, the third order force constants are extracted (as solutions to a set of linear equations). Phonon lifetimes are related to the third order force constants through the application of Fermi’s golden rule.

Integrating the modal thermal conductivity over the Brouillin zone, under the relaxation time approximation to the phonon Boltzmann transport equation (BTE), yields the lattice thermal conductivity. This full procedure is made concrete with the ShengBTE package Li et al. (2014).

To extend the above procedure to a crystalline alloy, approximations are necessary. Following Garg et al. Garg et al. (2011), we use the virtual crystal approximation (VCA). Within this approximation, two paths can be taken. One can compositionally average the pseudopotentials for the constituent atoms, and then proceed with the usual procedure. Alternatively, one can calculate the harmonic and third order force constants for the unalloyed crystalline versions of the constituent atoms, take the mass normalized compositional average and then proceed to calculate the phonon properties:

[TABLE]

where $A$ is a placeholder for the harmonic force constants, the third order force constants, the atomic masses and the lattice constants Antonius and Louie (2016). We have followed both VCA procedures, and find negligible difference in the phonon properties (see supplementary material).

The penultimate step in the alloy calculation is to include the effect of mass disorder. Again, following Garg’s work, the phonon lifetimes are modified under Matthiessen’s rule using the theory established by Tamura Tamura (1983) to treat isotope scattering as an elastic perturbation through the coupling parameter $g$ defined as:

[TABLE]

where $f_{i}$ is the concentration and the scattering time scales as $\tau_{disorder}~{}\sim\frac{1}{g}$. Garg et al. went a step further to estimate the anharmonic shifts do due disorder through supercell calculations. Feng et al. used molecular dynamics to show that the application of Matthiessen’s rule leads to an overestimation of thermal conductivity in SiGe due to neglecting four and five-phonon processes Feng et al. (2015). Our experimental results will show that the harmonic mass disorder approximation under Matthiessen’s rule is a sufficient approximation. We note that this procedure will not capture the frequency shifts that can be observed in the SiGe Raman spectra Feldman et al. (1966); Sui et al. (1993) (see the supplementary material). It is expected that these modes do not significantly contribute to thermal conductivity, as their group velocities are small and their lifetimes have been reduced by mass disorder scattering. The virtual crystal approximation is expected to break down when the mass disorder takes on a correlation length on the order of the probing length scale Mendoza et al. (2015).

The DFT calculation parameters used in this work are the following: for the DFPT portion, a 16 $\times$ 16 $\times$ 16 Monkhorst-Pack $k$ mesh with a kinetic energy cutoff of 50 Ry and a convergence criteria of 1E-12 Ry is used. For the supercell calculations, a 4 $\times$ 4 $\times$ 4 supercell was used such that third order force constants up to the fifth nearest neighbor could be obtained and only wavefunctions at the gamma point were calculated. Both (Si,Ge).pz-bhs.UPF and (Si,Ge).pz-n-nc.UPF pseudopotentials were tested yielding a negligible difference between thermal conductivity estimates (see supplementary material). The DFPT calculations were done with a 6 $\times$ 6 $\times$ 6 $q$ mesh. Interpolation was done on a 48 $\times$ 48 $\times$ 48 $q$ mesh with a Gaussian smearing parameter of 0.1 for the Kronecker delta approximation to yield convergence of the thermal conductivity. All calculations were done with the quantum-ESPRESSO package Giannozzi et al. (2009). The input files and the properties are available in the supplementary material.

II.2 Solving the Boltzmann Transport Equation

Given the bulk phonon properties of Si93.4Ge6.6, we now turn to the study of classical size effects on thermal transport in the reflection mode TTG geometry. The diffusive temperature profile has been solved previously in order to analyze the temperature signal using TTG for opaque materials Johnson et al. (2012). For the experimental conditions of a spatially periodic heat source defined by wavevector $q=\frac{2\pi}{\lambda}$, the temperature is given by $T(x,z,t)=T_{0}+T_{0}e^{iqx}h(z,t)$ in complex form, and this serves as a definition of the non-dimensional temperature $h$. The temperature $T_{0}$ is the background equilibrium temperature of the system, for example the room temperature. The heating by the laser is incorporated with a volumetric heat generation term, given by the functional form:

[TABLE]

where $U_{0}$ represents the energy per unit area deposited into the substrate by the pulse, and $\beta$ is the inverse penetration depth of the heating profile. The derivation found in Johnson et al. (2012) takes into consideration different in-plane and cross-plane thermal conductivities, however the experimental signal is sensitive to the in-plane thermal conductivity. For simplicity, we show the derivation for an isotropic system, where the Fourier heat conduction equation simplifies to:

[TABLE]

with the initial and boundary conditions given by:

[TABLE]

which assumes an adiabatic surface at $z=0$, and that the system starts at equilibrium prior to the energy deposited by the laser. We present the solution in the Laplace transformed domain for convenience:

[TABLE]

We intend to utilize this Fourier heat conduction temperature profile in our variational solution of the BTE. Taking the inverse Laplace transform of this yields the temperature as a function of the depth into the substrate and time:

[TABLE]

where the surface heating profile is:

[TABLE]

II.2.1 Temperature integral equation

We begin with the spectral Boltzmann transport equation under the relaxation time approximation (RTA):

[TABLE]

where $g_{\omega}$ is the phonon energy density per unit frequency interval per unit solid angle above the reference background energy, related to the distribution function as $g_{\omega}=\frac{\hbar\omega D(\omega)}{4\pi}(f_{\omega}-f_{0}(T_{0}))$. $\mathbf{v_{\omega}}$ is the group velocity, $\tau_{\omega}$ is the relaxation time, and $g_{0}$ is the equilibrium energy density, given by $g_{0}\approx\frac{1}{4\pi}C_{\omega}(T-T_{0})$ in the linear response regime. The sinusoidal heating profile in the $x$-direction (in-plane), given by the pulse form $Q_{\omega}(x,z,t)=\delta(t)\mathrm{e}^{iqx}\tilde{Q_{\omega}}(z)$, means we can expect that the spectral and equilibrium energy densities to also obey a sinusoidal profile $g_{\omega}=\mathrm{e}^{iqx}\tilde{g_{\omega}}$ and the equilibrium distribution will simplify accordingly to $\tilde{g_{0}}=\frac{C_{\omega}T_{0}}{4\pi}h(z,t)$. By inputting this in-plane sinusoidal profile and utilizing the Laplace transform (denoted by the $\hat{}$ symbol) in the time domain, the BTE simplifies to:

[TABLE]

where we have defined $\eta_{\omega}=q\Lambda_{\omega}$. For convenience, we define the parameter $V=\frac{1+s\tau_{\omega}+i\eta_{\omega}\mu_{x}}{\Lambda_{\omega}\mu_{z}}$ to group the variables in a compact form for the following solution of the BTE:

[TABLE]

The boundary conditions are taken to be:

[TABLE]

The first boundary condition takes an imaginary blackbody wall at length $L$ into the substrate at the background temperature to account for the semi-infinite substrate, where this length will limit to infinity. The second boundary condition provides the adiabatic boundary condition with diffuse scattering, where $\sigma=\frac{1}{\pi}\int d\Omega\Theta(\mu_{z})\mu_{z}\hat{\tilde{g}}_{\omega}(z=0,s,\mu_{x},-\mu_{z})$, which is proportional to the specular heat flux approaching the surface. We have utilized the Heaviside step function to reduce the integration over the solid angle only to consider phonons approaching the surface. Solving the boundary conditions, and taking the artificial length $L$ to infinity yields the formal solution to the BTE for the spectral energy density in terms of the equilibrium energy density:

[TABLE]

where we have defined the following solid angle integral function:

[TABLE]

The first term represents phonons moving towards the surface of heating at $z=0$, whereas the second term represents phonons moving away from the surface.

The temperature equation can be derived by utilizing the equilibrium condition obtained by integrating Eq. 13 with respect to frequency and the solid angle Majumdar (1993). The equilibrium condition in this case can be expressed as:

[TABLE]

Performing the solid angle integral, and inputting the expression for the non-dimensional temperature expression $\hat{\tilde{g_{0}}}=\frac{C_{\omega}T_{0}}{4\pi}\hat{h}(z,s)$, we obtain the integral equation for the temperature distribution:

[TABLE]

This is an integral equation in the spatial variable $z$ for the non-dimensional temperature in the Laplace domain, which after solving, requires an inverse Laplace transform in order to obtain the full temperature solution in the time domain. For the thermal distribution, the spectral heat generation takes the form:

[TABLE]

Note that $\frac{C_{\omega}}{C}$ is a weighting of the contribution of a given mode to heat generation under the assumption of thermalized distribution and is different than the form found in Hua and Minnich (2014). While other distributions can be taken, we utilize this form in order to compare to the Fourier heat conduction solution.

II.2.2 Variational solution

Given the mathematical challenges in finding a closed solution to Eq. 16, we opt for a simpler path. The insight is to take the known Fourier heat conduction solution (Eqs. 6, 7) as a starting point for the variational trial function. The simplest trial function is to take the diffusive temperature profile and allow just the thermal diffusivity to be a variational parameter. In general, the size effects exhibited by the BTE will affect both the temporal as well as the spatial distributions of the temperature. However, the simple variational solution that varies only one parameter, the thermal diffusivity, performs admirably by approximately solving for the thermal decay from the BTE over a broad range of grating period length scales. We proceed by taking the Fourier heat conduction solution of Eq. 6 as a trial function and use the thermal diffusivity as the variational parameter.

To solve for the variational parameter, we can utilize mathematical optimization methods such as least squares on the error residual of the temperature equation Chiloyan et al. (2016a), or impose a physical condition that we wish the trial function to satisfy. Here, we impose that the trial function must satisfy energy conservation taken over the control volume of the semi-infinite substrate over all time, analogous to the condition utilized for the thin film TTG geometry Chiloyan et al. (2016b). This mathematical condition can be obtained by integrating the BTE of Eq. 15 over the solid angle and frequency, and then also over the depth variable $z$ as well as over all time to yield:

[TABLE]

This statement says that the total energy per unit area perpendicular to the $z$-axis deposited in the semi-infinite substrate initially (left hand side of Eq. 17) must be equal to the total energy that moves away in the in-plane direction. The in-plane heat flux is obtained by utilizing the spectral energy density of Eq. 13, and integrating over the frequency and solid angle $\hat{\tilde{q}}_{x}(z,s)=\int d\omega\int d\Omega\Theta v_{\omega}\mu_{x}\hat{\tilde{g}}_{\omega}(z,s,\mu_{x},\mu_{z})$ to obtain the in-plane heat flux:

[TABLE]

where we have defined the solid angle integral function:

[TABLE]

Inserting the heat flux expression of Eq. 19 into the energy conservation statement of Eq. 18, and inputting the variational trial function of the Fourier heat conduction solution of Eq. 6 as well as the thermal distribution for the heat generation rate, we can solve for the effective thermal conductivity after cleaning up some of the solid angle integrals. We obtain a form similar in structure to the results from the thin film TTG Chiloyan et al. (2016b) and the one-dimensional limit of the TTG Chiloyan et al. (2016a):

[TABLE]

where $\mathrm{Kn}_{\omega}=\Lambda_{\omega}\beta$ and $f$ and $g$ are the kernels that weigh a given mode’s contribution to effective thermal conductivity under the imposed size effects, explicitly given as

[TABLE]

We have defined the following solid angle integral functions:

[TABLE]

If we take the limit of $\mathrm{Kn}_{\omega}\to 0$, i.e. the case of very long penetration depth, the solid angle integrals vanish as $\psi_{n}(\eta_{\omega},\mathrm{Kn}_{\omega}\to 0)\propto\mathrm{Kn}_{\omega}$, and we recover the one-dimensional TTG limit as in this case the substrate essentially starts off at a uniform temperature, and we recover the previously derived effective thermal conductivity Chiloyan et al. (2016a). Note that information concerning the spectral contribution to heat capacity in needed, unlike prior work Yang and Dames (2013), in the equation for effective thermal conductivity. The more interesting case for this problem is the reduction to surface heating, i.e. $\mathrm{Kn}_{\omega}\to\infty$. In this case, the kernel functions simplify to:

[TABLE]

For the general case of arbitrary penetration depth, the solid angle integral functions can be calculated analytically, which allows for a fully analytical effective thermal conductivity for any penetration depth into the substrate.

II.2.3 Comparison between the Variational Solution and Monte Carlo Simulations

To study the effect of the optical penetration depth in the case of a diffuse surface boundary condition, we first plot the kernels $f$ and $g$ as a function of $\eta$ for the extremal limits of $\mathrm{Kn}_{\omega}$. The one dimensional limit of $\mathrm{Kn}_{\omega}\to 0$ and the surface heating limit of $\mathrm{Kn}_{\omega}\to\infty$ define the envelope of curves for which the kernels for arbitrary values of the penetration depth must lie between. As the Knudsen number increases, the size effect due to the optical penetration depth increases, which physically results in a decrease of the effective thermal conductivity. This occurs due to the decrease in the numerator kernel $f$, and the increase of the denominator kernel $g$. Figure 2 shows that one dimensional limit and the surface heating limit are practically indistinguishable, indicating that the effective thermal conductivity due to a diffuse boundary experiences weak effects from the optical penetration depth.

Utilizing the derived kernels to calculate the effective thermal conductivity for Si93.4Ge6.6, we show in Figure 3 the effective thermal conductivity in the various limits. Note that the effective thermal conductivity is quite similar in the one dimensional limit and in the surface heating limit. As expected, when the thermal grating period is much smaller than the optical penetration depth, the effective thermal conductivity takes on values of the one dimensional limit, as the transport is mostly in-plane. In the opposite case, when the grating period is much larger than the optical penetration depth, the effective thermal conductivity limits to the surface heating limit.

Figure 3 demonstrates that the variational technique predicts that transport has a weak dependence on the optical penetration depth, a consequence of the kernels’ weak dependence on optical penetration depth. As such, the one-dimensional limit of the TTG Chiloyan et al. (2016a) geometry is sufficient to characterize the dependence of effective thermal conductivity on grating period.

In the limit of large grating periods, the thin film TTG limits to the Fuchs-Sondheimer Sondheimer (1952) problem of in-plane transport, and there is still a reduction of the effective thermal conductivity due to the finite size of the membrane. In contrast, for the reflection TTG, the limit of large grating period yields the bulk thermal conductivity, regardless of the optical penetration depth. Thus a modified Fourier approach will fail to capture the details of a thermal decay due to a localized heat source (i.e. delta function in space and time in a semi-infinite geometry). In this case, the transport at short times (on the order of the dominant relaxation times) is initially ballistic and given sufficient scattering, the transport becomes diffusive. The variational method, using the Fourier temperature profile as input, reveals that the thermal conductivity that best recovers this behavior is bulk. This can be understood as a consequence of the constraint imposed by the equilibrium condition of Eq. 18, which dictates the behavior of the variational temperature profile in the large time limit where transport is diffusive. An example of this limitation is presented in the supplementary material. To ensure that this limitation is not present in the current experimental study, we compare against established Monte Carlo simulations of the RTA-BTE Péraud and Hadjiconstantinou (2011, 2012). As is seen in Figure 4, agreement at a grating period of 100 nm and an optical penetration depth of 10 nm and for a grating period of 10 um and an optical penetration depth of 1 um is observed.

As our experimental conditions operate at penetration depths on the order of 1 um for Si93.4Ge6.6 Braunstein et al. (1958), with grating periods of between 1 and 13.5 um, by the pigeonhole principle, we can move forward with our variational solutions.

III Experiment

III.1 Sample specifications

The deposition of the SiGe sample was done by metal-organic chemical vapor deposition (MOCVD). Briefly, SiH4 and GeH4 enter the reactor, which break up into Si, Ge, and H2 from exposure to high temperatures (750-800C). The composition is controlled by tunning the flow rates of SiH4 and GeH4. A single crystal sample consisting of 93.4% Si, 6.6% Ge with a thickness of 6 um on a Si wafer was used for this work.

III.2 Transient Thermal Grating

Transient grating spectroscopy is a variant on four-wave-mixing spectroscopic techniques that can measure thermal transport dynamics over a well-defined in-plane length scale. In this technique, two pump laser pulses (515 nm, 60 ps FWHM) are crossed at the surface of the sample, where they interfere to yield a sinusoidal intensity pattern. Absorption by the sample creates a matching temperature profile, which evolves as a function of time through in-plane and cross-plane transport. The time dynamics of this “transient grating” are measured by the diffraction of a quasi-continuous probe beam (532 nm), and phase-specific information is extracted through heterodyned detection of the TTG signal by superposition of the diffracted signal with a reference beam (local oscillator) derived from the probe beam source. The signal is detected using a fast photodiode (Hamamatsu C5658, 1 GHz bandwidth) and recorded on an oscilloscope (Tektronix TDS 7404, 4GHz bandwidth). Specific details of the optical setup can be found elsewhere Maznev et al. (1998); Johnson et al. (2012); Vega-Flick et al. (2015).

The TTG signal will in principle have both real and imaginary field contributions due to “amplitude-grating” and “phase-grating” responses, respectively. The phase grating contributions contains decay components that correspond to thermal expansion and the imaginary part of the thermoreflectance and acoustic oscillations corresponding to the impulsive stimulation of surface acoustic waves (SAWs), whereas the amplitude-grating response only contains one term corresponding to the real part of the thermoreflectance Maznev et al. (1998). Analysis of the amplitude-grating contribution is simpler due to the single contribution, and so this term was isolated during the measurements by optimizing the heterodyne phase to minimize the SAW signal which only appears in the phase-grating response.

III.3 Results

All measurements of the Si93.4Ge6.6 sample were conducted at room temperature. Figure 5 depicts the TTG measurements alongside the prediction from the variational solution using properties obtained from first principle calculations following Section II. We have used an optical penetration depth of 1500 nm, according to Braunstein et al. (1958). The effect of uncertainty in the penetration depth is discussed in the supplementary material. The agreement is remarkable, considering we are simply fitting the TTG measurements to the Fourier-based temperature profiles (Eq. 8) to extract effective thermal conductivity. This agreement persists for a range of grating periods, from $\sim$ 13.5 to 1 um. Example fits of the TTG data with comparisons to the variational predictions are found in the supplementary material.

IV Discussion and Outlook

Before delving into the nuances of the work, a quick reminder of what we have done. We calculated the first principles phonon properties to match the exact composition of the sample studied experimentally. We then used these properties and the variational solution to the BTE to predict (without any fitting required) the recorded observable of TTG experiments, the temperature decay. In doing so, we uncover excellent agreement between the effective thermal conductivities of theory and experiment.

One of the first explanations of size effects in SiGe grew out of the observation of frequency dependence in TDTR measurements Koh and Cahill (2007). This explanation relied on the application of thermal penetration depth, $L_{tpd}\sim\sqrt{\frac{\alpha_{bulk}}{\omega}}$, as a heuristic approximation to estimate the magnitude of the deviation from a bulk thermal conductivity. For Si93.4Ge6.6, $\alpha_{bulk}=1.2358\text{E-5}~{}\text{ m}^{2}\text{/s}$, with 10 MHz, yields a $L_{tpd}\sim 1$ um. Under this approximation, we can take the MFP thermal conductivity accumulation function at 1 um, yielding $0.4k_{bulk}=7.3$ W/mK. From our TTG results, using $L_{tpd}~{}\sim\lambda$, we find $0.25k_{bulk}=4.5$ W/mK. By this same argument, frequency dependence should also be observed in silicon with $\alpha_{bulk}=8.8\text{E-5}~{}\text{ m}^{2}\text{/s}$, which at 10 MHz, yields a $L_{tpd}\sim 3$ um and $\sim 0.7k_{bulk}=98$ W/mK from the MFP accumulation function Esfarjani et al. (2011), but $k_{exp}>120$ W/mK for the same frequency range is often reported Wilson and Cahill (2014, 2015). The reason for this discrepancy has not been satisfactorily resolved Ding et al. (2014); Vermeersch and Shakouri (2016). For example, the results of Hua et al. Hua et al. (2016) and Wilson et al. Wilson and Cahill (2014) suggest that the reported thermal conductivity obtained from a TDTR measurement is dependent upon the interface conductance, indicating that this thermal conductivity can no longer be interpreted as an intrinsic property of the material. Meanwhile, the penetration depth argument has been used to interpret frequency dependence in BB-FDTR measurements Regner et al. (2013), suggesting that this tool could be used for phonon MFP spectroscopy. The next natural step in the interpretations of deviations from bulk required theory to move beyond the Heaviside cutoff of the thermal penetration depth and obtain a gray suppression function from solving the gray BTE Minnich (2012); Regner et al. (2014); Zeng et al. (2014). This function is then used as a kernel in the effective thermal conductivity integral. This picture has also turned out to be an oversimplification Collins et al. (2013); Chiloyan et al. (2016a). Here we show that a fully spectral solution to the BTE is required to characterize the effective conductivity. This progression from penetration depth to gray suppression to fully spectral interpretations is shown in Figure 6.

In contrast to the interpretation of thermal penetration depth of TDTR, the length scales in TTG do not depend on the intrinsic value of a material’s transport coefficient, and are therefore physically well-defined independent variables. Although the information concerning the optical penetration depth is required, this is well within current characterization technology Fuyuki et al. (2005). Given that the variational solutions to the 1D and surface heating TTG geometries predict approximately the same effective thermal conductivity dependence on grating period, we have theoretical bounds on the observed experimental decay curves from which the the transport regime can be determined (i.e. purely 1D, finite penetration depth, or surface heating). In doing so, we have presented a theoretical framework that is falsifiable, given that experimental deviations from theory can be understood as departures from the approximations used in this work: the VCA, the RTA-BTE and the specific trial solution for the temperature profile used in the variational method. These approximations can be lifted, as will be shown in future works. With the methodology presented here, the TTG can be used to study in-plane transport in opaque thin films that require a supporting substrate.

As TDTR measurements are sensitive to the cross-plane transport, the TTG provides a complementary tool for measuring in-plane transport. The variational method can be extended to more complicated geometries, such as layered systems with interfaces, ideally suited for providing insight into the interpretations of TDTR and TTG measurements. Such an extension would provide a path towards unifiying the interpretations of the measurements from TDTR and TTG.

V Conclusion

Our TTG experimental results augmented with DFT-based modeling and the variational BTE solution indicate that this experimental geometry is capable of meeting the predictive criteria necessary for studying size effects on thermal transport in complex materials, such as the SiGe alloy studied here. Interesting questions can now be asked, such as in what systems or at what length scales can we expect to find a breakdown of the VCA. This geometry will likely prove useful in the study of systems where the relaxation time approximation fails. The TTG experiment provides a path towards tabletop studies of the microscopic properties of thermal transport.

VI Acknowledgments

We acknowledge Jiawei Zhou and Bai Song for helpful discussions. This work was done as part of 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, Basic Energy Sciences (BES), under Award DE-SC0001299 / DE-FG02-09ER46577.

VII Appendix A

VIII Appendix B

IX Appendix C

X Appendix E

XI Appendix F