Modeling quasi-ballistic transient thermal transport with spatially sinusoidal heating: a McKelvey-Shockley flux approach

TL;DR

This paper extends the McKelvey-Shockley flux method to model transient ballistic phonon transport under sinusoidal heating, accurately capturing non-diffusive effects and aligning well with experimental and Boltzmann solutions.

Contribution

The work introduces an extended McKelvey-Shockley flux approach that includes inelastic scattering and supports full phonon dispersion, providing a simple yet accurate model for ballistic to diffusive heat transport.

Findings

The extended model agrees with experimental data.

It accurately reproduces Boltzmann transport solutions.

Inelastic scattering is crucial for the diffusive limit.

Abstract

Ballistic phonon effects, arising on length scales comparable to the mean-free-path, result in non-diffusive heat flow and alter the thermal properties of materials. Simple theoretical models that accurately capture non-diffusive transport physics are valuable for experimental analysis, technology design, and providing physical insight. In this work, we utilize and extend the McKelvey-Shockley (McK-S) flux method, a simple and accurate framework, to investigate ballistic effects in transient phonon transport submitted to a spatially sinusoidal heating profile, simulating a transient thermal grating. We begin by extending a previous McK-S formulation to include inelastic scattering, then obtain an analytical solution in the single phonon energy case (gray approximation), and after show how this approach can readily support a full phonon dispersion and mean-free-path distribution. The results agree with experimental data and compare very well to solutions of the phonon Boltzmann transport equation in the diffusive and weakly quasi-ballistic transport regimes. We discuss the role of ballistic and non-equilibrium physics, and show that inelastic scattering is key to retrieving the heat equation solution in the diffusive limit. Overall the McK-S flux method, which takes the form of a diffusion-like equation, proves to be a simple and accurate framework that is applicable from the ballistic to diffusive transport regime.

Equations

$$\frac{1}{v_x^{+}(ϵ)}\frac{\partial I_Q^{+}(ϵ)}{\partial t}+\frac{\partial I_Q^{+}(ϵ)}{\partial x}$$

$$\frac{1}{v_x^{+}(ϵ)}\frac{\partial I_Q^{-}(ϵ)}{\partial t}-\frac{\partial I_Q^{-}(ϵ)}{\partial x}$$

$$T(x,t)=\frac{{\int }_0^{\infty }T(x,t,ϵ)C_V(ϵ)\mathrm{d}ϵ}{{\int }_0^{\infty }{C}_V(ϵ)\mathrm{d}ϵ},$$

$$I_Q(x,t)={\int }_0^{\infty }{I}_Q(x,t,ϵ)\mathrm{d}ϵ.$$

$$I_Q(ϵ)+{\tau }_Q(ϵ)\frac{\partial I_Q(ϵ)}{\partial t}=-\kappa (ϵ)\frac{\partial T(ϵ)}{\partial x},$$

$${\tau }_Q(ϵ)\frac{{\partial }^{2}T(ϵ)}{\partial t^2}+2\frac{\partial T(ϵ)}{\partial t}+\frac{T(ϵ)}{{\tau }_Q(ϵ)}-\left[\frac{\partial T_0}{\partial t}+\frac{T_0}{{\tau }_Q(ϵ)}\right]$$

$$=D_{\mathrm{ph}}(ϵ)\frac{{\partial }^{2}T(ϵ)}{\partial x^2},$$

$${\int }_0^{\infty }\left[C_V(ϵ)\frac{\partial T(ϵ)}{\partial t}+\frac{\partial I_Q(ϵ)}{\partial x}\right]\mathrm{d}ϵ=0.$$

$$T_0=\frac{{\int }_0^{\infty }T(ϵ)C_V(ϵ)/{\tau }_Q(ϵ)\mathrm{d}ϵ}{{\int }_0^{\infty }{C}_V(ϵ)/{\tau }_Q(ϵ)\mathrm{d}ϵ},$$

$$T(x,t=0)=\Delta T\cos (2\pi x/L),$$

$$T(x,t)=\Delta T\cos (qx)e^{-q^2{D}_{\mathrm{ph}}t},$$

$${\tau }_Q\frac{{\partial }^{2}T}{\partial t^2}+\frac{\partial T}{\partial t}=D_{\mathrm{ph}}\frac{{\partial }^{2}T}{\partial x^2}.$$

$$T(x,t)=\Delta T\cos (qx)e^{-t^{\prime }}\left[\cosh (\beta t^{\prime })+\frac{\sinh (\beta t^{\prime })}{\beta }\right],$$

$$T(x,t)=\frac{\Delta T\cos (qx)}{2}\times$$

$$\left[\left(1+\frac{1}{\beta }\right)e^{-(1-\beta )t^{\prime }}+\left(1-\frac{1}{\beta }\right)e^{-(1+\beta )t^{\prime }}\right],$$

$$T(x,t)=\Delta T\cos (qx)e^{-q^2{D}_{\mathrm{ph}}t},$$

$$T(x,t)=\Delta T\cos (qx)e^{-\frac{v_x^{+}}{\lambda }t}\cos (q{v}_x^{+}t).$$

$$T(x,t)=\Delta T\cos (qx)\cos (q{v}_x^{+}t).$$

$${\tau }_Q(ϵ)\frac{{\partial }^{2}T(t,ϵ)}{\partial t^2}$$

$$=\left[\frac{\partial T_0(t)}{\partial t}+\frac{T_0(t)}{{\tau }_Q(ϵ)}\right],$$

$$T_{i+1}(ϵ)=$$

$$+$$

$$v_x^{+}(ϵ)$$

$$\lambda (ϵ)$$

Full text

I Introduction

Understanding the detailed thermal transport properties of materials is important for many applications, such as thermoelectric energy conversion and heating in nanoelectronic devices. Phonons, the dominant heat carrier in semiconductors and insulators, can have mean-free-path (MFP) values ranging over orders of magnitude, from nanometers to micrometers Henry2008 ; Esfarjani2011 ; Li2012a ; Li2012b ; Jain2015 . When the characteristic length scale of the system is comparable to the phonon MFP ballistic effects appear. This can modify the thermal transport properties, such as a reduction in thermal conductivity. Recent experiments have exploited this phenomenon to probe deviations from expected classical results and gain insight into the fundamental phonon properties of materials, such as the MFP distribution Koh2007 ; Siemens2010 ; Minnich2011 ; Regner2013 ; Johnson2013 ; Wilson2014 ; Hu2015 ; Hoogeboom-Pot2015 ; Zeng2015 ; Johnson2015 ; Cuffe2015 .

One such experiment, known as transient thermal grating (TTG), utilizes short, interfering laser pulses to create a 1D spatially-periodic, sinusoidal heating pattern with period $L$ (see Fig. 1) Johnson2013 ; Johnson2015 ; Cuffe2015 . By monitoring and analyzing the temperature decay versus time, using a diffracted probe laser, the thermal diffusivity (or thermal conductivity) can be extracted for different values of $L$. TTG has the advantages of being contactless, meaning there is no thermal transport across heterojunctions, and amenable for theoretical treatment given the relatively simple heating pattern. When $L$ is much longer than the phonon MFP transport is diffusive and the traditional heat equation is valid. When $L$ is similar to the phonon MFP transport is quasi-ballistic and a more rigorous approach that captures ballistic and non-equilibrium effects is needed to model and analyze experiments.

Several theoretical models have been developed to understand and analyze TTG experiments Maznev2011 ; Collins2013 ; Hua2014a ; Ramu2014 ; Minnich2015 ; Vermeersch2016b ; Chiloyan2016 ; Zeng2016 , and to relate raw experimental data to physically meaningful properties. Most models are based on the phonon Boltzmann transport equation (BTE), which is a rigorous semi-classical transport framework. Solving the BTE can, in general, be computationally challenging, however more compact solutions can be obtained for specific problems. Simple models are important for routine experimental analysis and technology design, as well as providing physical insight. The McKelvey-Shockley (McK-S) flux method Mckelvey1961 ; Shockley1962 was previously demonstrated to efficiently and accurately treat steady-state and transient phonon transport in all transport regimes from ballistic to diffusive Maassen2015a ; Maassen2015b ; Maassen2016 ; Kaiser2017 . The McK-S framework can be viewed as a simple version of the rigorous BTE, which has the benefits of physical transparency and computational efficiency, and could prove useful for treating complicated problems. Interestingly, the McK-S equations were shown to be mathematically equivalent to well-known diffusion equations Maassen2015a ; Maassen2015b ; Maassen2016 ; Kaiser2017 , which when solved with the appropriate physical boundary conditions provide surprisingly good agreement with more rigorous approaches. It is important to explore the accuracy of McK-S under different conditions, by benchmarking this approach against the BTE, to understand its limits and provide improvements.

In this paper, we apply the McK-S flux method to investigate the transient temperature decay when subjected to an initial spatially sinusoidal heating profile (i.e. TTG setup). We extend the McK-S equations to include inelastic scattering, which is later shown to be important for transient, spectral (non-gray or energy-resolved) problems. An analytical temperature solution is obtained in the single phonon energy (gray) approximation. We then calculate temperature using a detailed material model, including a full phonon dispersion and MFP distribution, which is compared to experiment and the BTE. The McK-S flux method is found to agree well with measurements and rigorous calculations, particularly when transport is diffusive or weakly quasi-ballistic.

The paper is structured as follows. In Sections II and III we present the TTG model structure and the McK-S flux method. Section IV presents the temperature solutions along with comparisons to experiment and the phonon BTE. In Section V we discuss our results and compare theoretical approaches. Lastly, Section VI summarizes our findings.

II Model Structure

We model the TTG temperature profile by considering a semiconductor/insulator film initially submitted to a spatially sinusoidal temperature pattern along the $x$ direction, with heating period $L$ (see Fig. 1). The film has a thickness $l$ along the $z$ direction, and is infinite in the $x$-$y$ plane. We assume that the initial heating and temperature is uniform along the $y$ and $z$ directions. The thickness of the film, $l$, enters the calculations through surface scattering, arising when phonons impinge on the top/bottom film boundaries leading to enhanced thermal resistance (discussed later). After the initial heating profile is imposed on the sample, we calculate the transient thermal response using the McK-S flux method described in the next section.

III Theoretical Approach

We employ the McKelvey-Shockley flux method Mckelvey1961 ; Shockley1962 to treat phonon transport, which was previously shown to apply from the ballistic to the diffusive transport regimes Maassen2015a ; Maassen2015b ; Maassen2016 ; Kaiser2017 . McK-S discretizes the phonon distribution into forward and reverse moving streams (relative to the transport direction, $x$). The governing equations are:

[TABLE]

where $I_{Q}^{\pm}(x,t,\epsilon)$ are the directed heat fluxes of the forward ($+$) and reverse ($-$) moving phonons, evaluated at a specific phonon energy $\epsilon$. $I_{Q,0}^{+}(x,t,\epsilon)$=$I_{Q,0}^{-}(x,t,\epsilon)$ are the fluxes of a (fictitious) local equilibrium distribution characterized by the temperature $T_{0}(x,t)$, which comes from a relaxation time approximation treatment of scattering. $v_{x}^{+}(\epsilon)$ is the angle-averaged $x$-projected velocity of phonons at energy $\epsilon$. $\lambda(\epsilon)$ is the mean-free-path for backscattering, defined as the average distance traveling along $x$ before backscattering (scattering from a forward to a reverse flux or vice versa).

Heat current and heat density are obtained from $I_{Q}(x,t,\epsilon)$=I_{Q}^{+}(x,t,\epsilon)$$-$$I_{Q}^{-}(x,t,\epsilon) and $Q(x,t,\epsilon)$=Q^{+}(x,t,\epsilon)$$+$$Q^{-}(x,t,\epsilon), where $Q^{\pm}(x,t,\epsilon)$=$I_{Q}^{\pm}(x,t,\epsilon)/v_{x}^{+}(\epsilon)$. Eqns. (1)-(2) are different from previous formulations of the McK-S approach Maassen2015a ; Maassen2015b ; Maassen2016 ; Kaiser2017 , since we have not yet imposed conservation of energy (will come later). Note that the explicit $x$ and $t$ dependence of the directed fluxes, heat density and heat current is not always shown for clarity, but is implied.

III.1 Energy-integrated temperature and heat current

The total, or energy-integrated, heat density ($Q$) and heat current ($I_{Q}$) are obtained by simply integrating $Q(\epsilon)$ and $I_{Q}(\epsilon)$ over energy. For this study we prefer to deal with temperature instead of heat density, which is possible by assuming all temperature variations are small perturbations relative to a constant background temperature, $T_{\rm ref}$. In this case, temperature is related to heat density via heat capacity ($C_{V}$): $Q=C_{V}T$, where $C_{V}=\int_{0}^{\infty}C_{V}(\epsilon)\,{\rm d}\epsilon$, $C_{V}(\epsilon)=\epsilon\,D(\epsilon)\,\partial n_{\rm BE}(\epsilon,T)/\partial T|_{T_{\rm ref}}$ is the energy-resolved heat capacity, $D(\epsilon)$ is the phonon density of states, and $n_{\rm BE}$ is the Bose-Einstein distribution (note that $T$ is the temperature relative to $T_{\rm ref}$). We can also define an energy-resolved temperature using $Q(\epsilon)=C_{V}(\epsilon)T(\epsilon)$. Using this relation, combined with $Q=\int_{0}^{\infty}Q(\epsilon)\,{\rm d}\epsilon$, we can express the energy-integrated temperature and heat current as:

[TABLE]

Eqns. (3)-(4) represent quantities of interest, which can be compared to measurements.

III.2 Governing equations for heat current and temperature

By subtracting Eqns. (1)-(2), we obtain a constitutive equation for heat current:

[TABLE]

where $\kappa(\epsilon)=C_{V}(\epsilon)D_{\rm ph}(\epsilon)$, $D_{\rm ph}(\epsilon)=\lambda(\epsilon)v_{x}^{+}(\epsilon)/2$ and $\tau_{Q}(\epsilon)=\lambda(\epsilon)/(2v_{x}^{+}(\epsilon))$ are the energy-dependent bulk thermal conductivity, phonon thermal diffusivity and thermal relaxation time, respectively.

By adding Eqns. (1)-(2) and using Eq. (5), we arrive at the following equation for temperature:

[TABLE]

where $T_{0}(x,t)$ is the local equilibrium temperature of the equilibrium flux $I_{Q,0}^{\pm}(x,t,\epsilon)$. Note that, unlike $T$, $T_{0}$ is not energy dependent since it describes an equilibrium distribution. Eqns. (5)-(6), for heat current and temperature, are equivalent to solving Eqns. (1)-(2). Next, we show how to determine $T_{0}$ by imposing conservation of energy.

III.3 Conservation of energy

In the absence of a heat sink or source, conservation of energy can be written as a continuity equation:

[TABLE]

By adding Eqns. (1)-(2) and imposing Eq. (7), we find:

[TABLE]

which is used to determine $T_{0}(x,t)$. Thus, one must solve Eqns. (6)-(8) simultaneously to obtain $T(x,t,\epsilon)$, then heat current is evaluated using Eq. (5), and finally the energy-integrated quantities, $T(x,t)$ and $I_{Q}(x,t)$, are computed with Eqns. (3)-(4). While in general Eqns. (6)-(8) must be solved self-consistently, later we show it is not necessary for this problem.

The McK-S equation for temperature, Eq. (6), captures ballistic and non-equilibrium effects, and can treat transport continuously from the ballistic to the diffusive limit, while taking the form of an efficient diffusion-like equation. A different McK-S formulation for transient heat transport was previously presented Maassen2015b ; Maassen2016 ; Section V.1 discusses their differences and demonstrates why the current approach is advantageous for the treatment of transient, spectral problems.

IV Results

We begin with the classical heat equation, which will serve as a comparison to the McK-S flux method. Then we solve the McK-S equation for temperature considering i) a single phonon energy (gray approximation), and ii) a more realistic detailed phonon model used for comparison with the phonon BTE and experimental data.

IV.1 Heat equation solution

Solving the heat equation, $\partial T/\partial t\,$=$\,D_{\rm ph}\,\partial^{2}T/\partial x^{2}$, with the initial temperature profile,

[TABLE]

we obtain Johnson2013

[TABLE]

where $q=2\pi/L$. Classically, temperature decays exponentially over a characteristic time that scales with $q^{-2}\propto L^{2}$ (experimentally controlled) and $D_{\rm ph}^{-1}$ (a material parameter).

IV.2 McK-S solution: single phonon energy (gray approximation)

Next we solve Eq. (6) for a single phonon energy (gray approximation). Imposing conservation of energy within the gray approximation is equivalent to setting the integrand of Eq. (7) to zero, which gives $T_{0}(x,t)=T(x,t,\epsilon)$. Using this result with Eq. (6), we obtain the hyperbolic heat equation (HHE):

[TABLE]

It was previously shown how the HHE could be derived from the McK-S equations Maassen2015b ; Maassen2016 . Solving the HHE with the initial condition given by Eq. (9) and $\partial T/\partial t|_{t=0}=0$ gives

[TABLE]

where $t^{\prime}=t/(2\tau_{Q})$ and $\beta=\sqrt{1-q^{2}\lambda^{2}}$. $q\lambda$, or equivalently $\lambda/L$, controls the phonon transport regime. When $q\lambda\rightarrow 0$ ($L\gg\lambda$) transport is diffusive, when $q\lambda\sim 1$ ($L\sim\lambda$) transport is quasi-ballistic, and when $q\lambda\rightarrow\infty$ transport is ballistic.

Fig. 2a shows the TTG temperature decay curves calculated with McK-S (Eq. (12)) and the heat equation (Eq. (10)) for different $q\lambda$. For $q\lambda<1$ the solutions are monotonically decaying, and there is close agreement between both approaches. In the diffusive limit, as $q\lambda\rightarrow 0$, both methods are identical (see Appendix A). For $q\lambda>1$ notable differences appear, in particular the McK-S solution shows oscillatory behavior. These oscillations can be understood by considering that at $t$=[math] a sinusoidal phonon distribution is created, with half the phonons moving to the left (right) with velocity $-v_{x}^{+}$ ($+v_{x}^{+}$). When $L$ is comparable to $\lambda$ scattering does not fully randomize the phonon distribution, which retains some of its traveling wave / oscillatory behavior.

In Fig. 2b, we compare the temperatures obtained with McK-S and the more rigorous phonon BTE Collins2013 . Both approaches agree for small $q\lambda$, when approaching the diffusive regime. For larger $q\lambda$ the BTE also presents temperature oscillations, however much less pronounced than with McK-S. This difference stems from an assumption in the McK-S flux method that all phonons (at a given energy $\epsilon$) travel at the same angle-averaged $x$-projected velocity ($v_{x}^{+}$), which leads to stronger oscillations. With the BTE phonons can travel at any angle giving a range of $x$-projected velocities, which effectively washes out the oscillations. This is clear in the ballistic limit: the McK-S solution is $T\propto\cos(qv_{x}^{+}t)$ (see Appendix A), while the BTE solution is $T\propto\sin(2qv_{x}^{+}t)/(2qv_{x}^{+}t)$ Collins2013 . Thus, we find McK-S performs well in the diffusive and weakly quasi-ballistic regime, and deviates from the BTE in the strongly quasi-ballistic and ballistic regimes.

Interestingly, at early times both McK-S and BTE agree and show a slower temperature decay compared to the heat equation, consistent with experimental observations of a decreased diffusivity / thermal conductivity. The gray approximation is good for illustrating the transport physics, however a quantitative comparison to experiment requires a more realistic phonon model.

IV.3 McK-S solution: detailed phonon model (Si film)

Next we solve the McK-S expression for temperature (Eqns. (6) and (8)) using a realistic phonon model to investigate the TTG temperature decay in a 400 nm-thick Si film.

Our Si model uses a full phonon dispersion calculated from density functional theory and a mean-free-path distribution obtained by fitting phenomenological scattering models for boundary, impurity and Umklapp phonon-phonon scattering to experimental data. The details of this model, which can reproduce the measured thermal conductivity of bulk Si within 15% from 5 K to 300 K, are provided in Ref. Maassen2015a . From this detailed phonon information we extract $v_{x}^{+}(\epsilon)$ and $\lambda(\epsilon)$ (see Def_parameters for definitions of these quantities), which are then converted to $\tau_{Q}(\epsilon)=\lambda(\epsilon)/(2v_{x}^{+}(\epsilon))$ and $D_{\rm ph}(\epsilon)=\lambda(\epsilon)v_{x}^{+}(\epsilon)/2$ used in Eqns. (6) and (8). To treat a Si film, we include surface scattering using a simple expression which reduces the bulk mean-free-path for backscattering: $1/\lambda_{\rm film}(\epsilon)$=$1/\lambda(\epsilon)$+$1/(\beta l)$, where $l=400$ nm is the thickness of the film and $\beta=2.21$ is a parameter fitted to give the measured film thermal conductivity, which is $\approx\,$62.5% the bulk thermal conductivity Johnson2013 .

We note that a rigorous treatment of TTG in a film would require solving a 2D heat transport problem that considers temperature variations along both $x$ and $z$ (cross-plane direction), and explicitly includes the film thickness. Here we consider 1D transport along the grating direction, $x$, which is a common approximation Collins2013 ; Hua2014a , and capture the effect of film thickness via surface scattering. The extension of the McK-S method to higher dimensions will be the focus of future work.

Fig. 3a presents the TTG temperature profiles of a 400 nm-thick Si film, calculated with McK-S (solid lines) and the heat equation (dashed lines), with varying heating period $L$. Details of the numerical solution of the McK-S temperature equation are provided in Appendix B. $T(t)$ decays more rapidly with decreasing $L$, as expected classically. In fact, the McK-S temperature solution, $T_{\rm McK-S}$, appears to decay only slightly slower than the heat equation solution, $T_{\rm HE}$. To better observe the differences between both methods, Fig. 3b shows the relative difference in temperature. As $L$ decreases and ballistic effects appear, $T_{\rm McK-S}$ deviates further from $T_{\rm HE}$ ($T_{\rm McK-S}\rightarrow T_{\rm HE}$ when $L\rightarrow\infty$). No temperature oscillations are observed with the shortest heating period $L=600$ nm, because the effective contribution of phonons at different energies leads to a monotonically decaying curve (note that $\lambda_{\rm film}$ varies from $\approx 1$ nm to $\approx 1$ $\mu$m).

By fitting the heat equation solution (Eq. (10)) to our McK-S temperature profiles, we extract an effective diffusivity or thermal conductivity, $\kappa_{\rm eff}$, as is done experimentally Johnson2013 . Fig. 4 presents the calculated (blue line) and measured (markers) $\kappa_{\rm eff}/\kappa_{\rm bulk}$ versus heating period $L$ (experimental data taken from Ref. Johnson2013 ). For large $L$, when transport is diffusive, $\kappa_{\rm eff}$ approaches $\kappa_{\rm film}$ (which is less than $\kappa_{\rm bulk}$ because of surface scattering). For $L<5$ $\mu$m, $\kappa_{\rm eff}$ rolls off as ballistic effects become prominent. Quasi-ballistic phonons decay more slowly than diffusive phonons due to fewer scattering events, which results in a reduced effective thermal conductivity. Note that the heat equation predicts a constant $\kappa_{\rm eff}=\kappa_{\rm film}$ for all $L$.

Overall McK-S generally agrees with the measured data, but shows a $\kappa_{\rm eff}$ roll-off at smaller $L$ compared to experiment. This can be the result of our adopted Si film model or a limitation of the McK-S flux method. To test the latter we utilize a rigorous BTE solution to the TTG problem developed by Hua and Minnich Hua2014a , implemented with the same Si film model. Comparing the McK-S (solid line) and BTE (dashed line) solutions in Fig. 4 we observe only small differences. This indicates that our Si film model, calibrated to bulk Si and using a simple expression for surface scattering, may not accurately represent the samples in Ref. Johnson2013 and is the reason for the discrepancy with the measured data. A comparison of McK-S and BTE temperature profiles is presented in Fig. 5. We find that the simple McK-S flux method compares remarkably well to the more rigorous BTE.

In this study, we have calculated the TTG temperature decay and effective thermal conductivity, given a particular phonon dispersion and MFP. We note, however, that it is possible to solve the inverse problem to reconstruct the MFP distribution given the measured $\kappa_{\rm eff}$ Minnich2012 ; Minnich2015 ; Cuffe2015 . This commonly involves calculating the so-called suppression function from the model, which dictates how the phonons with different MFP contribute to $\kappa_{\rm eff}$, and depends on the specific experimental setup and the characteristic length scale (i.e. the heating length, $L$, in the case of TTG). The McK-S approach could also be used to reconstruct the MFP distribution from measured $\kappa_{\rm eff}$, and may be particularly useful for complicated setups that cannot be easily solved with the BTE.

V Discussion

V.1 McK-S flux method: inelastic versus elastic scattering

Previous studies showed how the McK-S equations could be equivalently expressed as the hyperbolic heat equation (HHE), which when solved with the correct physical boundary conditions captures ballistic and non-equilibrium effect Maassen2015b ; Maassen2016 . To obtain the HHE from Eqns. (1)-(2), conservation of energy is imposed at every $\epsilon$, meaning that the integrand of Eq. (7) is set to zero. Consequently all energies are decoupled and phonon scattering is elastic, as collisions attempt to restore equilibrium independently at each $\epsilon$.

In the case of a transient, spectral treatment, for example using a full phonon dispersion and mean-free-path distribution, previous studies Maassen2015b ; Maassen2016 calculated temperature from Eq. (3) with $T(x,t,\epsilon)$ coming from the HHE. It was pointed out that the effective (i.e. energy integrated) heat transport properties calculated from gray approximation solutions (such as the HHE) do not agree with the heat equation in the diffusive transport limit Vermeersch2016a . We note that this problem only occurs with spectral calculations; the individual gray approximation solutions, evaluated at a specific energy, obey the heat equation in the diffusive case (see Appendix A). This issue arises because phonons at each $\epsilon$ decay independently, according to their respective time constant. To address this point, in this work we impose conservation of energy via Eq. (7), which requires solving Eq. (6) instead of the HHE for $T(x,t,\epsilon)$. This effectively couples phonons of all energies, as scattering attempts to bring the whole phonon distribution to equilibrium at temperature $T_{0}$. Since the continuity equation is not imposed for every $\epsilon$, energy is not conserved for each individual phonon channel. As a result, the approach presented in this paper (Eqns. (6) and (7)) includes inelastic scattering.

Fig. 5a shows the TTG temperature decay for $L=20$ $\mu$m (diffusive transport case) calculated with the approach described in this work (“McK-S”, solid line), the elastic scattering version of McK-S (“McK-S (elastic)”, red dashed line), and the heat equation (“Heat eq”, green dashed line). McK-S is found to agree with the heat equation, unlike the elastic scattering version of McK-S. In the diffusive limit the HHE gives $T\propto\exp(-q^{2}D_{\rm ph}(\epsilon)\,t)$, which once integrated over energy gives a temperature profile that is a sum of different exponentials. We determine that when performing transient, spectral calculations inelastic scattering, which couples all phonons, is needed to correctly retrieve the heat equation in the diffusive limit.

V.2 Local temperature $T(x,t)$ versus local equilibrium temperature $T_{0}(x,t)$

When solving the McK-S temperature equation both $T(t)$ and $T_{0}(t)$ are obtained. In the diffusive transport limit $T$ and $T_{0}$ are almost identical (see Fig. 5a), since the phonon distribution is near equilibrium. As $L$ decreases $T$ and $T_{0}$ differ (see Fig. 5b), since quasi-ballistic transport results in a non-equilibrium phonon population. While $T_{0}$ has on occasion been interpreted as the local temperature of the system, we view $T_{0}$ as a mathematical device related to the relaxation time approximation treatment of scattering and conservation of energy.

We consider that the definition of $T(x,t)$, given by Eq. (3), is physically meaningful and valid when the temperature variations are small relative to the background temperature. This definition obeys linear response theory, as one would expect physically. Note that $T(x,t,\epsilon)$, which appears in Eq. (3), only has an energy dependence when the phonon distribution deviates from equilibrium, leading to $T\neq T_{0}$.

VI Summary

The McKelvey-Shockley (McK-S) flux method was employed to investigate ballistic effects in transient phonon transport under a spatially sinusoidal heating condition, similar to the transient thermal grating (TTG) experiment. We developed a more general formulation of the McK-S equations that captures inelastic scattering, extending previous studies Maassen2015b ; Maassen2016 .

In the case of a single phonon energy (gray approximation) the McK-S approach is equivalent to solving the hyperbolic heat equation, for which an analytical temperature solution was derived. When transport is diffusive ($L\gg\lambda$) the McK-S/HHE solution agrees with the heat equation, and deviations appear when transport is quasi-ballistic ($L\sim\lambda$). Interestingly, for $L<2\pi\lambda$ the temperature oscillates versus time, a feature shared with the more rigorous phonon BTE. This behavior is more pronounced with McK-S, which originates from the assumption that all phonons (at a given energy) travel at the same angle-averaged $x$-projected velocity.

A full phonon dispersion and mean-free-path distribution was used to calculate the TTG temperature decay, and extract the effective thermal conductivity, of a 400 nm-thick Si film. In the diffusive limit, the temperature profile agrees with the heat equation. As $L$ is decreased, a slower decay compared to the heat equation is observed, resulting in a reduced effective thermal conductivity. Comparison with experiment shows good agreement, however the measured onset of ballistic behavior appears at slightly larger $L$, which is attributed to the Si film model since the BTE compares very well to McK-S.

We discussed differences between the previous and current formulations of the McK-S flux method, and demonstrated how the latter is better for treating transient, spectral calculations. We argued the physical basis of our definition of temperature, and highlight the difference with the local equilibrium temperature related to our use of the relaxation time approximation.

We find the McK-S flux method is a simple, efficient and physically transparent approach that provides similar accuracy to the BTE at a lesser computational demand. Extensions of this method to enable the treatment of a broader class of problems, such as 2D/3D transport, will be the focus of future work.

Acknowledgements.

This work was partially supported by DARPA MATRIX (Award No. HR0011-15-2-0037) and NSERC (Discovery Grant RGPIN-2016-04881).

Appendix A Limiting cases of the HHE solution

Here we verify the diffusive and ballistic limits of the HHE solution (Eq. (12)).

Diffusive transport limit. We begin by rewritting Eq. (12) in terms of exponentials, which gives

[TABLE]

where $t^{\prime}=t/(2\tau_{Q})$, $\beta=\sqrt{1-q^{2}\lambda^{2}}$, and $q=2\pi/L$. In the diffusive transport limit, we have $\lambda/L\rightarrow 0$, $q\lambda\rightarrow 0$, and $\beta\approx(1-q^{2}\lambda^{2}/2)$. Using this in Eq. (13), we find

[TABLE]

where we used $\tau_{Q}=\lambda/(2v_{x}^{+})$ and $D_{\rm ph}=\lambda v_{x}^{+}/2$. Thus, in the diffusive limit, we retrieve the heat equation solution (Eq. (10)).

Ballistic transport limit. In the ballistic limit $\lambda/L\rightarrow\infty$, $q\lambda\rightarrow\infty$, and $\beta\rightarrow iq\lambda$. Using this in Eq. (12), along with the identities $\cosh(ix)=\cos(x)$ and $\sinh(ix)=\sin(x)$, we obtain

[TABLE]

Interestingly, taking the limit $L\rightarrow 0$ gives Eq. (15), which is different from assuming $\lambda\rightarrow\infty$, which gives the ballistic solution:

[TABLE]

In the ballistic limit the temperature is oscillatory. This occurs since both forward/reverse-moving phonons are initially excited, with the profile $\cos(qx)$, and propagate at a constant velocity without decaying which creates a standing wave.

Appendix B Numerical solution of the McK-S temperature equation

Since the initial applied heating profile is spatially sinusoidal, given by Eq. (9), the temperature solution will also be spatially sinusoidal: $T(x,t)=T(t)\,\cos(qx)$. Inserting this into Eq. (6) and dividing by $\cos(qx)$, we obtain:

[TABLE]

where $T_{0}(t)$ is related to $T(t,\epsilon)$ through Eq. (8). To numerically solve Eq. (17) we discretize time, $t_{i}$, and use the finite difference formulas ${\rm d}T/{\rm d}t|_{t_{i}}\approx(T_{i}-T_{i-1})/\Delta t$ and ${\rm d}^{2}T/{\rm d}t^{2}|_{t_{i}}\approx(T_{i+1}-2T_{i}+T_{i-1})/\Delta t^{2}$, where $T_{i}\equiv T(t_{i})$ and $\Delta t=t_{i}-t_{i-1}$ is the uniform grid spacing. Finite difference discretization of Eq. (17) gives:

[TABLE]

To calculate the temperature at the next time step, $T(t_{i+1})$, we need the temperature at the previous two time steps ($T(t_{i})$ and $T(t_{i-1})$) and, importantly, also the local equilibrium temperature ($T_{0}(t_{i})$ and $T_{0}(t_{i-1})$).

The initial conditions are $T(t=0)=\Delta T$, $\partial T/\partial t|_{t=0}=0$, $T_{0}(t=0)=\Delta T$, and $\partial T_{0}/\partial t|_{t=0}=0$, which once discretized become: $T(t_{1})=\Delta T$, $T(t_{0})=\Delta T$, $T_{0}(t_{1})=\Delta T$, and $T_{0}(t_{0})=\Delta T$, where $t_{1}$ corresponds to zero time and the first recorded time step. To solve the McK-S equation for temperature, the procedure is as follows: 1) With the initial conditions calculate $T(t_{2},\epsilon)$ using Eq. (18). 2) Use Eq. (8) to obtain $T_{0}(t_{2})$ by performing a numerical integration over $\epsilon$. 3) Calculate the energy-integrated temperature $T(t_{2})$ using Eq. (3). 4) Advance to the next time step ($t_{3},\,t_{4},\,\ldots$) and repeat this procedure using information from the previous two time steps.