Experimental demonstration of a generalized Fourier's Law for non-diffusive thermal transport

TL;DR

This paper derives and experimentally validates a generalized Fourier's law that accurately describes non-diffusive thermal transport from ballistic to diffusive regimes, highlighting the importance of nonlocal effects and inhomogeneous source terms.

Contribution

The paper introduces a generalized Fourier's law derived from the Peierls-Boltzmann equation, accounting for nonlocal effects and inhomogeneous sources in thermal transport.

Findings

The generalized Fourier's law matches TDTR measurements in the non-diffusive regime.

Neglecting the generalized law leads to inaccurate thermal property measurements.

The law applies across ballistic to diffusive regimes for various geometries.

Abstract

Phonon heat conduction over length scales comparable to their mean free paths is a topic of considerable interest for basic science and thermal management technologies. Although the failure of Fourier's law beyond the diffusive regime is well understood, debate exists over the proper physical description of thermal transport in the ballistic to diffusive crossover. Here, we derive a generalized Fourier's law that links the heat flux and temperature fields, valid from ballistic to diffusive regimes and for general geometries, using the Peierls-Boltzmann transport equation within the relaxation time approximation. This generalized Fourier's law predicts that thermal conductivity not only becomes nonlocal at length scales smaller than phonon mean free paths, but also requires the inclusion of an inhomogeneous nonlocal source term that has been previously neglected. We demonstrate the validity of this generalized Fourier's law through direct comparison with time-domain thermoreflectance (TDTR) measurements in the nondiffusive regime without adjustable parameters. Furthermore, we show that interpreting experimental data without this generalized Fourier's law leads to inaccurate measurement of thermal transport properties.

Equations

$$\frac{\partial g_{\mu }(\mathbf{x},t)}{\partial t}+{\mathbf{v}}_{\mu }\cdot \nabla g_{\mu }(\mathbf{x},t)=-\frac{g_{\mu }(\mathbf{x},t)-g_0(T(\mathbf{x},t))}{{\tau }_{\mu }}+{\dot{Q}}_{\mu }(\mathbf{x},t),$$

$$\frac{\partial E(\mathbf{x},t)}{\partial t}+\nabla \cdot \mathbf{q}(\mathbf{x},t)=\dot{Q}(\mathbf{x},t),$$

$${\Lambda }_{\mu x}\frac{\partial {\tilde{g}}_{\mu }}{\partial x}+{\Lambda }_{\mu y}\frac{\partial {\tilde{g}}_{\mu }}{\partial y}+{\Lambda }_{\mu z}\frac{\partial {\tilde{g}}_{\mu }}{\partial z}+(1+i\eta {\tau }_{\mu }){\tilde{g}}_{\mu }=C_{\mu }\Delta \tilde{T}+{\tilde{Q}}_{\mu }{\tau }_{\mu },$$

$$\xi$$

$$\rho$$

$$\zeta$$

$${\Lambda }_{\mu \xi }\frac{\partial {\tilde{g}}_{\mu }}{\partial \xi }+{\alpha }_{\mu }{\tilde{g}}_{\mu }=C_{\mu }\Delta \tilde{T}+{\tilde{Q}}_{\mu }{\tau }_{\mu },$$

$${\tilde{g}}_{\mu }(\xi ,\rho ,\zeta ,\eta )$$

$${\tilde{g}}_{\mu }(\xi ,\rho ,\zeta ,\eta )=P_{\mu }{e}^{-{\alpha }_{\mu }\frac{\xi }{{\Lambda }_{\mu \xi }}}+{\int }_{\Gamma }\frac{C_{\mu }\Delta \tilde{T}+{\tilde{Q}}_{\mu }{\tau }_{\mu }}{|{\Lambda }_{\mu \xi }|}{e}^{-{\alpha }_{\mu }\left|\frac{{\xi }^{\prime }-\xi }{{\Lambda }_{\mu \xi }}\right|}d{\xi }^{\prime },$$

$$P_{\mu}=\left\{\begin{array}[]{rl}P^{+}_{\mu}e^{\alpha_{\mu}\frac{L_{1}}{\Lambda_{\mu\xi}}}&\mbox{ if $v_{\mu\xi}>0$}\\ P^{-}_{\mu}e^{\alpha_{\mu}\frac{L_{2}}{\Lambda_{\mu\xi}}}&\mbox{ if $v_{\mu\xi}<0$}\end{array}\right.$$

$$\Gamma\in\left\{\begin{array}[]{rl}[L_{1},\xi)&\mbox{ if $v_{\mu\xi}>0$}\\ \left(\xi,L_{2}\right]&\mbox{ if $v_{\mu\xi}<0$}\end{array}\right..$$

$$\sum _{\mu }{v}_{\mu \xi }\frac{\partial {\tilde{g}}_{\mu }}{\partial \xi }+i\eta \sum _{\mu }{\tilde{g}}_{\mu }=\sum _{\mu }{\tilde{Q}}_{\mu },$$

$${\tilde{q}}_{\mu \xi }=P_{\mu }{v}_{\mu \xi }{e}^{-{\alpha }_{\mu }\frac{\xi }{{\Lambda }_{\mu \xi }}}+{\int }_{\Gamma }{\tilde{Q}}_{\mu }({\xi }^{\prime },\rho ,\zeta ,\eta )e^{-{\alpha }_{\mu }\left|\frac{\xi -{\xi }^{\prime }}{{\Lambda }_{\mu \xi }}\right|}d{\xi }^{\prime }+{\int }_{\Gamma }\frac{C_{\mu }{v}_{\mu \xi }}{|{\Lambda }_{\mu \xi }|}\Delta \tilde{T}({\xi }^{\prime },\rho ,\zeta ,\eta )e^{-{\alpha }_{\mu }\left|\frac{\xi -{\xi }^{\prime }}{{\Lambda }_{\mu \xi }}\right|}d{\xi }^{\prime }.$$

$${\tilde{q}}_{\mu \xi }=-{\int }_{\Gamma }{\kappa }_{\mu \xi }(\xi -{\xi }^{\prime })\frac{\partial \tilde{T}}{\partial {\xi }^{\prime }}d{\xi }^{\prime }+B_{\mu }(\xi ,\rho ,\zeta ,\eta ),$$

$$B_{\mu }(\xi ,\rho ,\zeta ,\eta )$$

$${\kappa }_{\mu \xi }(\xi )=C_{\mu }{v}_{\mu \xi }{\Lambda }_{\mu \xi }\frac{e^{-{\alpha }_{\mu }\left|\frac{\xi }{{\Lambda }_{\mu \xi }}\right|}}{{\alpha }_{\mu }|{\Lambda }_{\mu \xi }|}.$$

$${\tilde{q}}_{\mu \xi }=-{\int }_{\Gamma }{C}_{\mu }{v}_{\mu \xi }{\Lambda }_{\mu \xi }\frac{\partial \tilde{T}}{\partial {\xi }^{\prime }}\delta (\xi -{\xi }^{\prime })d{\xi }^{\prime }=-{\kappa }_{\mu \xi }\frac{\partial \tilde{T}}{\partial \xi },$$

$$\underset{\Xi \to 0}{\lim }{\alpha }_{\mu }\approx 1,$$

$$\underset{\Xi ,\text{ Kn}\to 0}{\lim }{e}^{-{\alpha }_{\mu }\frac{\xi }{|{\Lambda }_{\mu \xi }|}}\approx 0,$$

$$\underset{\Xi ,\text{ Kn}\to 0}{\lim }\frac{e^{-{\alpha }_{\mu }\left|\frac{\xi -{\xi }^{\prime }}{{\Lambda }_{\mu \xi }}\right|}}{|{\Lambda }_{\mu x}|}\approx \delta (\xi -{\xi }^{\prime }).$$

$$-\sum _{\mu }{\kappa }_{\mu \xi }\frac{{\partial }^2\tilde{T}}{\partial {\xi }^2}+i\eta \sum _{\mu }{C}_{\mu }\Delta \tilde{T}=\sum _{\mu }{\tilde{Q}}_{\mu }.$$

$${\kappa }_x\frac{{\partial }^2\tilde{T}}{\partial x^2}+{\kappa }_y\frac{{\partial }^2\tilde{T}}{\partial y^2}+{\kappa }_z\frac{{\partial }^2\tilde{T}}{\partial z^2}+i\eta \Delta \tilde{T}\sum _{\mu }{C}_{\mu }=\sum _{\mu }{\tilde{Q}}_{\mu },$$

$${\tilde{q}}_x(x,\eta )=i\beta \tilde{T}(\eta )e^{i\beta x}\sum _{\mu x>0}\frac{{\kappa }_{\mu x}}{{\alpha }_{\mu }^2+{\Lambda }_{\mu x}^2{\beta }^2}+\sum _{\mu x>0}\frac{Q_{\mu }}{\delta }\frac{e^{i\beta x}{\alpha }_{\mu }{\Lambda }_{\mu x}}{{\alpha }_{\mu }^2+{\Lambda }_{\mu x}^2{\beta }^2},$$

$${\tilde{q}}_x(x,\eta )=i\beta \tilde{T}(\eta )e^{i\beta x}\sum _{{\mu }_x>0}\frac{{\kappa }_{\mu x}}{1+{\Lambda }_{\mu x}^2{\beta }^2}+\sum _{{\mu }_x>0}\frac{Q_{\mu }}{\delta }\frac{e^{i\beta x}{\Lambda }_{\mu x}}{1+{\Lambda }_{\mu x}^2{\beta }^2}.$$

$${\kappa }_x=\sum _{{\mu }_x>0}\frac{{\kappa }_{\mu x}}{1+{\Lambda }_{\mu x}^2{\beta }^2},$$

$${\tilde{q}}_x(x,f_y,f_z,\eta )=-{\int }_{L_1}^{L_2}{\kappa }_x(x-x^{\prime },f_y,f_z,\eta )\frac{\partial T}{\partial x^{\prime }}d{x}^{\prime }+\sum _{\mu }{\tilde{B}}_{\mu }(x,f_y,f_z,\eta ),$$

$${\kappa }_x(x,f_y,f_z,\eta )=\sum _{{\mu }_x>0,{\mu }_y,{\mu }_z}{\kappa }_{\mu x}\frac{e^{-\frac{1+i{\Xi }_{\mu }+i{f}_y{\Lambda }_{\mu y}+i{f}_z{\Lambda }_{\mu z}}{|{\Lambda }_{\mu x}|}x}}{(1+i{\Xi }_{\mu }+i{f}_y{\Lambda }_{\mu y}+i{f}_z{\Lambda }_{\mu z})|{\Lambda }_{\mu x}|}.$$

$$\frac{1}{G}=\frac{4}{{\sum }_{\mu }{C}_{\mu }^{Si}{v}_{\mu }^{Si}{T}_{\text{Si}\to \text{Al}}}-\frac{2}{{\sum }_{\mu }{C}_{\mu }^{Al}{v}_{\mu }^{Al}}-\frac{2}{{\sum }_{\mu }{C}_{\mu }^{Si}{v}_{\mu }^{Si}}$$

$${\tilde{q}}_{\mu x}=-{\int }_{\Gamma }{C}_{\mu }{v}_{\mu x}{\Lambda }_{\mu x}\frac{e^{-{\alpha }_{\mu }\left|\frac{x-x^{\prime }}{{\Lambda }_{\mu x}}\right|}}{{\alpha }_{\mu }|{\Lambda }_{\mu x}|}\frac{\partial \tilde{T}}{\partial x^{\prime }}d{x}^{\prime }+{\int }_{\Gamma }{\tilde{Q}}_{\mu }(x^{\prime })e^{-{\alpha }_{\mu }\left|\frac{x-x^{\prime }}{{\Lambda }_{\mu x}}\right|}d{x}^{\prime },$$

$$\Gamma\in\left\{\begin{array}[]{rl}[-\infty,\xi)&\mbox{ if $v_{\mu\xi}>0$}\\ \left(\xi,\infty\right]&\mbox{ if $v_{\mu\xi}<0$}\end{array}\right..$$

$${\tilde{q}}_x=-i\beta \Delta \tilde{T}(\eta )\sum _{{\mu }_x>0}{\int }_{-\infty }^{\infty }{C}_{\mu }{v}_{\mu x}{\Lambda }_{\mu x}\frac{e^{i\beta x^{\prime }}{e}^{-{\alpha }_{\mu }\left|\frac{x-x^{\prime }}{{\Lambda }_{\mu x}}\right|}}{{\alpha }_{\mu }|{\Lambda }_{\mu x}|}d{x}^{\prime }+\sum _{{\mu }_x>0}{\int }_{-\infty }^{\infty }\frac{Q_{\mu }}{\delta }{e}^{i\beta x^{\prime }}{e}^{-{\alpha }_{\mu }\left|\frac{x-x^{\prime }}{{\Lambda }_{\mu x}}\right|}d{x}^{\prime }.$$

Full text

I Introduction

Fourier’s law fails when a temperature gradient exists over a length scale comparable to or smaller than the mean free paths (MFPs) of heat carriers. In this regime, the heat flux and temperature fields may differ from the predictions of heat diffusion theory based on Fourier’s law. These discrepancies have been observed at a localized hotspot created by a doped resistor thermometer in a suspended silicon membranesverdrup_measurement_2001 and more recently in optical pump-probe experiments including soft x-ray diffraction from nanoline arrays,highland_ballistic-phonon_2007 ; siemens_quasi-ballistic_2010 transient grating,johnson_direct_2013 and thermoreflectance methods.cahill_nanoscale_2014 ; koh_frequency_2007 ; minnich_thermal_2011 ; regner_broadband_2012 ; english_mean_2013 ; hu_spectral_2015 ; Hua_2017 ; ravichandran_spectrally_2018 In particular, due to the absence of scattering the transport properties become nonlocal, in contrast to Fourier’s law in which the heat flux at a certain location is determined by the temperature gradient only at that location.

Lattice thermal transport in crystals is generally described by the Peierls-Boltzmann equation (PBE), first derived by Peierls,Peierls1929 from which the thermal conductivity is given in terms of the microscopic properties of phonons.Peierls1929 ; Peierls1955 However, solving the PBE for a general space-time dependent problem remains a challenging task due to the high dimensionality of the integro-differential equation. Thus, most prior works have determined solutions of the PBE in certain limiting cases.

Guyer and Krumhanslguyer_solution_1966 first performed a linear response analysis of the PBE, deriving a space-time-dependent thermal conductivity by assuming the Normal scattering rates were much larger than Umklapp scattering rates, and they applied their solution to develop a phenomenological coupling between phonons and elastic dilatational fields caused by lattice anharmonicity. Hardy and coworkers reported a rigorous quantum-mechanical formulation of the theory of lattice thermal conductivity using a perturbation method that included both anharmonic forces and lattice imperfections.hardy_energy-flux_1963 ; hardy_perturbation_1965 ; hardy_lowestorder_1965 This quantum treatment of lattice dynamics was then verified by both theoretically and experimentally demonstrating the presence of Poiseuille flow and the second sound in a phonon gas at low temperatures when Umklapp processes may be neglected.sussmann_thermal_1963 ; guyer_thermal_1966 ; hardy_phonon_1970 ; jackson_thermal_1971 ; beck_phonon_1974 The variational principle was also used to solve the PBE with Umklapp scattering incorporated.hamilton_variational_1969 ; srivastava_derivation_1976 Levinson developed a nonlocal diffusion theory of thermal conductivity from a solution of the PBE with three-phonon scattering in the low frequency limit.levinson_nonlocal_1980

Advances in computing power have enabled the numerical solution of the PBE with inputs from density functional theory, fully ab initio. For instance, bulk lattice thermal conductivities are now routinely computed from first principles using an iterative solution of the PBEward_ab_2009 ; broido_lattice_2005 ; li_shengbte:_2014 ; carrete_almabte_2017 ; omini_iterative_1995 or from variational approaches.fugallo_ab_2013 Chaputchaput_direct_2013 presented a direct solution to the transient linearized PBE with an imposed constant temperature gradient. Cepellotti and Marzaricepellotti_thermal_2016 introduced the concept of a ”relaxon”, an eigenstate of the symmetrized scattering operator of the PBE first used by Guyer et. al.guyer_solution_1966 and Hardyhardy_phonon_1970 . They applied this treatment to solve steady-state problems in two-dimensional systems with a constant temperature gradient imposed in one directioncepellotti_boltzmann_2017 .

Solving the PBE with the full collision operator, even in its linearized form, is difficult for complicated geometries. Therefore, various theoretical frameworks based on a simplified PBE have been developed to describe nonlocal thermal transport for general problems. Non-diffusive responses observed in experimentsjohnson_direct_2013 ; cahill_nanoscale_2014 ; koh_frequency_2007 ; regner_broadband_2012 ; Wilson2013 ; Dames2015 have been explained using the phonon-hydrodynamic heat equationtorres_emergence_2018 , a truncated Levy formalismVermeersch2015b , a two-channel model in which low and high frequency phonons are described by the PBE and heat equationMaznev2011 , and a Mckelvey-Shockley flux methodMaassen2015 . Methods based on solving the PBE under the relaxation time approximation (RTA), where each phonon mode relaxes towards thermal equilibrium at a characteristic relaxation rate, have been developed to investigate nonlocal transport in an infinite domainMahan_1988 ; hua_analytical_2014 ; allen_temperature_2018 ; Collins2013APL , a finite one-dimensional slabHua_semi_analytical_2015 ; Koh_2014 , and experimental configurations such as transient gratingRamu2014 ; Collins2013APL and thermoreflectance experimentsRegner2014 ; zeng_disparate_2014 ; Vermeersch_2015a . An efficient Monte Carlo scheme has been used to solve the PBE under the RTA in complicated geometries involving multiple boundariesPeraud:2011 ; Peraud:2012 ; hua_importance_2014 .

In the diffusion regime, Fourier’s law is a relation between heat flux and temperature fields, independent of the specific problem. In the nondiffusive regime, obtaining such a relation is more complicated because the transport is inherently nonlocal. The works described above solve the PBE for problems with specific boundary conditions or inputs. Thus, despite these efforts, a description generalizing Fourier’s law for arbitrary geometries and linking the heat flux and temperature fields in all transport regimes is not available.

Here, we derive a generalized Fourier’s law to describe non-diffusive thermal transport for general geometries using the linearized PBE within the RTA. The generalized Fourier’s law requires the inclusion of an inhomogeneous nonlocal term arising from the source or the boundary conditions of the particular problem. By including the inhomogeneous contribution to the heat flux, the space- and time-dependent thermal conductivity is independent of the specific geometry or inputs. This generalized Fourier’s law is validated by favorable comparisons with a series of TDTR measurements in the non-diffusive regime. We also show that neglecting the inhomogeneous contribution to the heat flux leads to inaccurate measurement of thermal transport properties in the non-diffusive regime. Our work provides a unified description of heat transport for a wide range of problems from ballistic to diffusive regimes.

II Theory

II.1 Governing Equation

We begin by briefly reviewing the derivation of the transport solution to the PBE. The mode-dependent PBE under the relaxation time approximation for transport is given by

[TABLE]

where $g_{\mu}(\mathbf{x},t)=\hbar\omega_{\mu}(f_{\mu}(\mathbf{x},t)-f_{0}(T_{0}))$ is the deviational energy distribution function at position $\mathbf{x}$ and time $t$ for phonon states $\mu$ ($\mu\equiv(\mathbf{q},s)$, where $\mathbf{q}$ is the wavevector and $s$ is the phonon branch index). $f_{0}$ is the equilibrium Bose-Einstein distribution, and $g_{0}(T(\mathbf{x},t))=\hbar\omega_{\mu}(f_{0}(T(\mathbf{x},t))-f_{0}(T_{0}))\approx C_{\mu}\Delta T(\mathbf{x},t)$, where $T(\mathbf{x},t)$ is the local temperature, $T_{0}$ is the global equilibrium temperature, $\Delta T(\mathbf{x},t)=T(\mathbf{x},t)-T_{0}$ is the local temperature deviation from the equilibrium value, and $C_{\mu}=\hbar\omega_{\mu}\frac{\partial f_{0}}{\partial T}|_{T_{0}}$ is the mode-dependent specific heat. Here, we assume that $\Delta T(\mathbf{x},t)$ is small such that $g_{0}(T(\mathbf{x},t))$ is approximated to be the first term of its Taylor expansion around $T_{0}$. Finally, $\dot{Q}_{\mu}(\mathbf{x},t)$ is the heat input rate per mode, $\mathbf{v}_{\mu}=(v_{\mu x},v_{\mu y},v_{\mu z})$ is the phonon group velocity vector, and $\tau_{\mu}$ is the phonon relaxation time.

To close the problem, energy conservation is used to relate $g_{\mu}(\mathbf{x},t)$ to $\Delta T(\mathbf{x},t)$ as

[TABLE]

where $E(\mathbf{x},t)=V^{-1}\sum_{\mu}g_{\mu}(\mathbf{x},t)$ is the total volumetric energy, $\mathbf{q}(\mathbf{x},t)=V^{-1}\sum_{\mu}g_{\mu}(\mathbf{x},t)\mathbf{v_{\mu}}$ is the directional heat flux, and $\dot{Q}(\mathbf{x},t)=V^{-1}\sum_{\mu}\dot{Q}_{\mu}(\mathbf{x},t)$ is the volumetric mode-specific heat input rate. Here, the sum over $\mu$ denotes a sum over all phonon modes in the Brillouin zone, and $V$ is the volume of the crystal. The solution of Eq. (1) yields a distribution function, $g_{\mu}(\mathbf{x},t)$, from which temperature and heat flux fields can be obtained using Eq. (2). Like the classical diffusion case, the exact expression of the temperature field varies from problem to problem. However, in a diffusion problem, the constitutive law that links the temperature and heat flux fields is governed by one expression, Fourier’s law. Here, we seek to identify a similar relation that directly links temperature and heat flux fields for non-diffusive transport, regardless of the specific problem.

To obtain this relation, we begin by rearranging Eq. (1) and performing a Fourier transform in time $t$ on Eq. (1), which gives

[TABLE]

where $\eta$ is the Fourier temporal frequency, and $\Lambda_{\mu x}$, $\Lambda_{\mu y}$ and $\Lambda_{\mu z}$ are the directional mean free paths along $x$, $y$, and $z$ directions, respectively. Equation (3) can be solved by defining a new set of independent variables $\xi$, $\rho$, and $\zeta$ such that

[TABLE]

where $\Lambda_{\mu}=\sqrt{\Lambda^{2}_{\mu x}+\Lambda^{2}_{\mu y}+\Lambda^{2}_{\mu z}}$. The Jacobian of this transformation is $\Lambda^{2}_{\mu x}/\Lambda^{2}_{\mu}$, a nonzero value. After changing the coordinates from $(x,y,z)$ to the new coordinate system $(\xi,\rho,\zeta)$, $(v_{\mu\xi}=v_{\mu x},0,0)$ is the set of elements for the velocity vector $\mathbf{v}_{\mu}$ in the new coordinates, and Eq. (3) becomes a first order partial differential equation with only one partial derivative

[TABLE]

where $\alpha_{\mu}=1+i\eta\tau_{\mu}$. Assuming that $\xi\in[L_{1},L_{2}]$, Eq. (5) has the following solution:

[TABLE]

$P^{+}_{\mu}$ and $P^{-}_{\mu}$ are functions of $\rho$, $\zeta$, $\eta$ and are determined by the boundary conditions at $\xi=L_{1}$ and $\xi=L_{2}$, respectively. Using the symmetry of $v_{\mu\xi}$ about the center of the Brillouin zone, i.e., $v_{\mu\xi}=-v_{-\mu\xi}$, Eqs. (6a) & (6b) can be combined into the following form:

[TABLE]

where

[TABLE]

and

[TABLE]

The energy conservation equation becomes

[TABLE]

where $v_{\mu\xi}\tilde{g}_{\mu}$ gives the mode-specific heat flux along the $\xi$ direction expressed as:

[TABLE]

Applying integration by parts to the third term in Eq. (11), we can write the heat flux per mode as:

[TABLE]

where

[TABLE]

is solely determined by the boundary condition and the volumetric heat input rate. $\kappa_{\mu\xi}(\xi)$ is the modal thermal conductivity along the $\xi$ direction given by

[TABLE]

Equation (12) is the primary result of this work. This equation links temperature gradient to the mode-specific heat flux for a general geometry. Since this constitutive equation of heat conduction is valid from ballistic to diffusive regimes, we denote it as a generalized Fourier’s law. It describes that for a specific phonon mode $\mu$, heat only flows in the $\xi$ direction in the new coordinate system ($\xi$, $\rho$, $\zeta$) since the velocities in $\rho$ and $\zeta$ directions are zero. To obtain the total heat flux in the original coordinate system, e.g. $q_{x}$, $q_{y}$, and $q_{z}$ in Cartesian coordinates, all the functions involved in Eq. (12) must be mapped from the coordinate system ($\xi$, $\rho$, $\zeta$) to ($x$, $y$, $z$). Analytical mappings exist for a few special cases that we will discuss shortly.

There are two parts in Eq. (12). The first part represents a convolution between the temperature gradient along the $\xi$ direction and a space-and time-dependent thermal conductivity, $\kappa_{\mu\xi}(\xi)$. As reported previously, this convolution indicates the nonlocality of the thermal conductivity.Mahan_1988 ; Koh_2014 ; allen_temperature_2018 However, a second term exists that is determined by the inhomogeneous term originating from the boundary conditions and source terms. Similar to the first term, the contribution from the heat input to the heat flux, given by $\int_{\Gamma}\tilde{Q}_{\mu}(\xi^{\prime},\rho,\zeta,\eta)e^{-\alpha_{\mu}\left|\frac{\xi-\xi^{\prime}}{\Lambda_{\mu\xi}}\right|}d\xi$, is nonlocal, meaning the contribution at a given point is determined by convolving the heat source function with an exponential decay function with a decay length of $\Lambda_{\mu x}$.

While the nonlocality of thermal conductivity was identified by earlier works on phonon transportguyer_solution_1966 ; guyer_thermal_1966 ; levinson_nonlocal_1980 ; koh_frequency_2007 ; Vermeersch2015b ; Mahan_1988 ; Allen_2018 ; allen_temperature_2018 , the contribution from the inhomogeneous term has been neglected. Recently, Allen and Perebeinosallen_temperature_2018 considered the effects of external heating and derived a thermal susceptibility based on the PBE that links external heat input to temperature response and a thermal conductivity that links temperature response to heat flux. However, their derived thermal susceptibilities and thermal conductivities are subject to the specific choice of the external heat input. In this work, we demonstrate that there exists a general relation between heat flux and temperature distribution without specifying the geometry of the problem. The space- and time-dependent thermal conductivity in the first term of Eq. (12) is independent of boundary conditions and heat input. The dependence of heat flux on the specific problem is accounted for by the inhomogeneous term.

II.2 Diffusive limit

Here we examine some specific limits of Eq. (12). First, in the diffusive regime, the spatial and temporal dependence of thermal conductivity disappears and asymptotically approaches a constant. To demonstrate this limit, we first identify two key nondimensional parameters in Eq. (12): Knudsen number, Knμ $\equiv\Lambda_{\mu\xi}\xi^{-1}$, which compares phonon MFP with a characteristic length, in this case $\xi^{-1}$, and a transient number, $\Xi_{\mu}\equiv\eta\tau_{\mu}$, which compares the phonon relaxation times with a characteristic time, in this case $\eta^{-1}$. In the diffusive limit, both $\Xi$ and Kn are much less than unity. Then, Eq. (7) is simplified to $C_{\mu}\Delta\tilde{T}$, and Eq. (12) becomes

[TABLE]

since in this limit we can perform the following simplifications

[TABLE]

[TABLE]

[TABLE]

The equation of energy conservation becomes

[TABLE]

Since $\frac{\partial}{\partial\xi}=\frac{\partial}{\partial x}+\frac{\partial}{\partial y}\frac{\Lambda_{\mu y}}{\Lambda_{\mu x}}+\frac{\partial}{\partial z}\frac{\Lambda_{\mu z}}{\Lambda_{\mu x}}$, Eq. (19) can be mapped back to Cartesian coordinates as

[TABLE]

where $\kappa_{i}=\sum_{\mu}C_{\mu}v_{\mu i}\Lambda_{\mu i}$ is the thermal conductivity along axis $i=x,y$, or $z$. Here, we assume that the off-diagonal elements of the thermal conductivity tensor are zero, i.e., $\kappa_{ij}=\sum_{\mu}C_{\mu}v_{\mu i}\Lambda_{\mu j}=0$ when $i\neq j$. Equation  (20) is the classical heat diffusion equation.

II.3 Generalized Fourier’s law in a transient grating experiment

We now check another special case of Eq. (12) by applying it to the geometry of a one-dimensional transient grating experiment.johnson_direct_2013 ; Hua_transport_2014 Since it is a 1D problem, $\xi$ in Eq. (12) is equivalent to $x$. In this experiment, the heat input has a spatial profile of $e^{i\beta x}$ in an infinite domain, where $\beta\equiv 2\pi/L$ and $L$ is the grating period. The boundary term vanishes, i.e., $\xi\in(-\infty,\infty)$, and both the distribution function and temperature field exhibit the same spatial dependence. Then, the total heat flux is expressed as

[TABLE]

where the total volumetric energy deposited on a sample is given by $\sum_{\mu}Q_{\mu}$, and the duration of the energy deposition is $\delta$. A derivation of Eq. (21) is given in Appendix A.

The time scale of a typical TG experiment is on the order of a few hundred nanoseconds while relaxation times of phonons are typically less than a nanosecond for many semiconductors at room temperature. Therefore, we assume that $\Xi\ll 1$, and Eq. (21) is simplified to

[TABLE]

which is consistent with what has been derived in our earlier work.Hua_transport_2014 ; Maznev2011 The first part of Eq. (22) represents the conventional understanding of nonlocal thermal transport, a Fourier type relation with a reduced thermal conductivity given by

[TABLE]

while the second part of the equation represents the contribution from the heat source to the total heat flux, which increases as the Knudsen number $\Lambda_{\mu x}/L$ increases. In a TG experiment, the presence of a single spatial frequency simplifies the convolutions in Eq. (12) into products, and the only time dependence comes from the temperature. Therefore, the decay rate of the measured transient temperature profile is directly proportional to the reduced thermal conductivity. In general, the spatial dependency of the temperature field is less complicated in a TG experiment than in other experiments, making the separation of the intrinsic thermal conductivity contribution from the inhomogeneous contribution easier.

II.4 Generalized Fourier’s law with infinite transverse geometries

The third special case considered here is when the $y$ and $z$ directions extend to infinity. The analytical mapping of Eq. (12) to Cartesian coordinates can be completed via Fourier transform in $y$ and $z$. After Fourier transform, Eq. (12) becomes

[TABLE]

where thermal conductivity $\kappa_{x}$ is given by

[TABLE]

$f_{y}$ and $f_{z}$ are the Fourier variables in the $y$ and $z$ directions, correspondingly, and $\tilde{B}_{\mu}(x,f_{y},f_{z},\eta)$ is the Fourier transform of $B_{\mu}(x,\rho,\zeta,\eta)$ with respect to $f_{y}$ and $f_{z}$. The exact expression of $\tilde{B}_{\mu}$ and a derivation of Eq. (24) are given in Appendix B. In this case, both temperature field and the inhomogeneous term have spatial and temporal dependence. Their dependence on the boundary conditions and heat source should be accounted for when extracting the intrinsic thermal conductivity from the observables such as total heat flux or an average temperature.

III Results

We now experimentally validate the generalized Fourier’s law by comparing the predicted and measured surface thermal responses to incident heat fluxes in TDTR experiments. We consider a sample consisting of an aluminum film on a silicon substrate. Phonon dispersions for Al and Si and lifetimes for Si were obtained from first-principles using density functional theory (DFT).lindsay_abinitio_semiconductors_2013 We assumed a constant MFP for all modes in Al; the value $\Lambda_{Al}=$ 60 nm is chosen to yield a lattice thermal conductivity $\kappa\approx$ 123 Wm*-1K-1* so that no size effects occur in the thin film. The justification of such an approach can be found in Ref. Hua_2017 .

In a TDTR experiment, the in-plane directions are regarded as infinite, thus Fourier transforms in the $y$ and $z$ directions are justified. The cross-plane direction in the substrate layer is semi-infinite. Therefore, the cross-plane heat flux in the substrate is described by Eq. (24) with $x\in[0,\infty)$. $P^{+}_{\mu}$ in Eq. (8) is determined by the interface conditions.Hua_2017

In the diffusion regime, energy conservation, Fourier’s law and the boundary conditions can fully describe a transport problem. In the nondiffusive regime, the replacement of Fourier’s law is the generalized Fourier’s law described in this work. If validated, this methodology allows the full prediction of the surface response for the wide variety of parameters employed in TDTR experiments, e.g., heating geometry, modulation frequency, and temperature.

To validate this methodology, we compare the calculated TDTR responses using the generalized Fourier’s law with pump-size-dependent TDTR measurements on the same Al/Si sample measured in Ref. Hua_2017 , where the $1/e^{2}$ diameter $D$ of a Gaussian pump beam was varied between 5 to 60 $\mu$m at different temperatures. As the pump size decreases and becomes comparable to the thermal penetration depth along the cross-plane direction ($x$-axis of the schematic in Fig. 1(a)), in-plane thermal transport is no longer negligible and requires a three-dimensional transport description. $P^{+}_{\mu}(\xi_{y},\xi_{z},\eta)$ is determined from the interface condition, and for a given $\xi_{y}$, $\xi_{z}$, and $\eta$ it is determined by the spectral transmissivity of phonons as given in Ref. Hua_2017 . In the same workHua_2017 , we have already used the PBE within the RTA to model the one-dimensional ($\xi_{y}=\xi_{z}=0$) thermal transport in a TDTR experiment with a uniform film heating and developed a method to extract the spectral transmissivity profile of phonons from measurements. We also provided evidence that elastic transmission of phonons across an interface was the dominant energy transmission mechanism for materials with similar phonon frequencies. Therefore, the measured transmissivity profile at room temperature given in Ref. Hua_2017 should be able to fully describe the interface conditions at other temperatures, and there are no adjustable parameters in the present simulations.

We compared the measured signals directly to predictions from the nonlocal transport governed by the generalized Fourier’s law and the strictly diffusive transport governed by Fourier’s law. To ensure a consistent comparison between the constitutive relations, the thermal conductivity of silicon is obtained using the same DFT calculations, and the interface conductance $G$ is given byzeng_nonequilibrium_2001

[TABLE]

where $T_{\text{Si}\rightarrow\text{Al}}$ is the phonon transmissivity from Si to Al. This expression was first derived by Chen and Zeng, which considers the non-equilibrium nature of phonon transport at the interface within the phonon transmissivity.zeng_nonequilibrium_2001

Figures 1 (a) & (b) show the total signal versus delay time with a pump size of 15 $\mu m$ at room temperature for the experiments and predictions from the generalized Fourier’s law and original Fourier’s law. As in Ref. Hua_2017 , the intensity of the shaded regions correspond to the likelihood of the measured transmissivity profile plotted in the inset of Fig. 1(b). A higher likelihood of a transmissivity profile is indicated by a higher intensity of the shaded area, and thus the PBE simulation using a transmissivity profile with higher likelihood better fits the experimentally measured TDTR signals. Excellent agreement between predictions from the generalized Fourier’s law and experiments are observed. On the other hand, Fourier’s law fails to accurately account for the experimental data, overestimating the phase and underestimating the amplitude after 2 ns in delay time. Note that the different transmissivity profiles in the inset of Fig. 1(b) give a value of interface conductance $G=223\pm 10$ Wm*-2K-1* using Eq. (26), and this deviation in $G$ leads to uncertainties in the TDTR signals that are within the linewidth of the plotted curves.

In Figs. 2(a) & (b), comparisons of phases at different pump sizes between the generalized Fourier’s law, original Fourier’s law and experimental data are given at 300 K. In Figs. 2(c) & (d), we compare the measured phase versus modulation frequency at a fixed pump size to predictions from the generalized Fourier’s law and original Fourier’s law at 100, 150, 200, 250 K. The data are given for two different delay times, 1.5 ns and 6 ns. The figure shows that predictions from the original Fourier’s law do not reproduce the experimental results. The deviation of Fourier’s law predictions from the experimental results becomes larger when the temperature decreases, modulation frequency increases, or pump size decreases, all indicating that non-diffusive effects increase with these changes. On the other hand, predictions from the generalized Fourier’s law agree well with the experimental measurements for the various temperatures, modulation frequencies ,and pump sizes, indicating its validity to describe the nonlocal thermal transport in different regimes.

IV Discussion

All the above comparisons between the simulations and experiments with different heating geometries and at different temperatures provide evidence that the generalized Fourier’s law is an appropriate replacement of Fourier’s law in the nondiffusive regime. We now use this formalism to examine TDTR measurements on boron arsenide (BAs).

Boron arsenide has recently attracted attention because of its ultra-high thermal conductivity determined from measurements based on the TDTR technique and reported by several groups.Lieaat8982 ; Kangeaat5522 ; Tianeaat7932 Moreover, pump-size-dependent measurements have also been conducted in an attempt to access information of the phonon MFPs in BAs.Kangeaat5522 The thermal conductivity measurements are based on interpreting the raw TDTR data as fit to a diffusion model based on Fourier’s law with thermal conductivity of BAs and interface conductance between Al and BAs as two fitting parameters.

However, due to the presence of phonons with long MFPs compared to the TDTR thermal length scale, Fourier’s law is no longer valid at the length scales probed by TDTR. As predicted by DFT-based PBE calculationslindsay_BAs_DFT_2013 ; feng_four-phonon_2017 , more than 70% of phonons in single crystalline BAs have mean free paths longer than 1 $\mu m$. Therefore, properly interpreting the data requires using the generalized Fourier’s law.

Equation (12) needs close examination to understand the microscopic information contained in the surface temperature responses measured in experiments. In a two-layer structure like the one used in TDTR, the second term in Eq. (12) does not vanish and has a non-local nature through the source term. This nonlocality implies that even though only the transient temperature at the metal surface is observed, the measurement contains information from the entire domain. We have demonstrated in Ref. Hua_2018 that the spectral distribution of the source term alters the surface temperature response. In other words, even though the first term in Eq. (12) remains the same, the observable at the surface is altered by the inhomogeneous source term originating from the interface.

To illustrate this point, we choose two transmissivity profiles T${}_{\text{BAs}\rightarrow\text{Al}}$ as shown in Figs. 3 (a) & (b). The profiles share a similar dependence on phonon frequency but with a different magnitude. The nominal interface conductance is calculated to be $115$ and $253$ MWm*-2K-1*, respectively, using Eq. (26). Along with the ab initio properties of BAs at room temperature, we calculate the TDTR responses at the Al surface with different pump sizes at three modulation frequencies. The calculated TDTR responses are then fit to the typical diffusion model based on Fourier’s law to extract the effective thermal conductivity, as was performed in Refs Lieaat8982 ; Kangeaat5522 ; Tianeaat7932 .

The results are shown in Fig. 3(c). The key observation from Fig. 3(c) is that the effects of modulation frequency and pump size on the effective thermal conductivity depend on the transmissivity profiles. A decrease in the effective thermal conductivity is observed in both profiles as the pump size decreases or the modulation frequency increases. However, the magnitude of the reduction and the absolute value compared to the bulk value depend on the transmissivity. While the effective thermal conductivity seems to be approaching the bulk value at a large pump size and low modulation frequency using profile 1, the effective thermal conductivity using profile 2 exceeds the bulk value under the same conditions. The reduction of the effective thermal conductivity using profile 1 as pump size decreases is less than 5% at a given modulation frequency, while the reduction using profile 2 can be as much as 40%.

Our calculations therefore indicate that in the nondiffusive regime, simply interpreting measurements from a method such as TDTR using Fourier’s law with a modified thermal conductivity may yield incorrect measurements. Not only is Fourier’s law unable to describe the nonlocal nature of thermal conductivity, but it also does not include the effects of inhomogeneous terms. Therefore, when interpreting TDTR measurements of high thermal conductivity materials, Fourier’s law is not the appropriate constitutive relation. In contrast, we have provided experimental evidence that the generalized Fourier’s law, Eq. (12), accurately describes thermal transport in the non-diffusive regime.

V Conclusions

In summary, we derived a generalized Fourier’s law using the Peierls-Boltzmann equation under the relaxation time approximation. The new constitutive relation consists of two parts, a convolution between the temperature gradient and a space- and time-dependent thermal conductivity, and an inhomogeneous term determined from boundary conditions and heat sources. By comparing predictions from this new constitution law to a series of time-domain thermorflectance measurements in the nondiffusive regime, we provide experimental evidence that the generalized Fourier’s law more accurately describes thermal transport in a range of transport regimes. We also show that interpreting nonlocal thermal transport using Fourier’s law can lead to erroneous interpretation of measured observables. To correctly extract microscopic phonon information from the observation of nonlocal thermal transport, it is necessary to separate the inhomogeneous contribution from the nonlocal thermal conductivity based on the generalized Fourier’s law developed here.

VI Acknowledgements

C. H. and L. L. acknowledge support from the Laboratory Directed Research and Development Program of Oak Ridge National Laboratory, managed by UT-Battelle, LLC, for the U.S. Department of Energy. A. J. M. acknowledges support from the National Science Foundation under Grant No. CBET CAREER 1254213. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231.

Appendix A Derivation of Eq. (21)

In a one-dimensional (1D) problem, Eq. (12) becomes

[TABLE]

where

[TABLE]

In a 1D transient grating experiment, both temperature profile and mode-specific heat input rate have a spatial dependence of $e^{i\beta x}$, i.e., $\tilde{T}(x,\eta)=e^{i\beta x}\tilde{T}(\eta)$ and $\tilde{Q}_{\mu}(x)=Q_{\mu}\delta^{-1}e^{i\beta x}$. Summing Eq. (27) over all the phonon modes and using the symmetry of $v_{\mu x}$ about the center of the Brillouin zone, i.e., $v_{\mu x}=-v_{-\mu x}$, the total heat flux is expressed as

[TABLE]

Now define $y=x^{\prime}-x$. Then the above equation becomes:

[TABLE]

Appendix B Derivation of Eq. (24)

When the $y$ and $z$ directions can be regarded as infinite, the analytical mapping to the Cartesian coordinates can be completed via Fourier transform in $y$ and $z$. To show it , we first define $g(x,y,z)=f(\xi,\rho,\zeta)$ with the coordinate transform given by Eq. (4). $G$ and $F$ are the functions after Fourier transform in $y$ and $z$. Using the affine theorem of two-dimensional Fourier transform, we obtain

[TABLE]

where $f_{y}$ and $f_{z}$ are the Fourier variables in the $y$ and $z$ directions, respectively.

Applying Eq. (30) to Eq. (11) gives

[TABLE]

where $\gamma_{\mu}=\eta\tau_{\mu}+f_{y}\Lambda_{\mu y}+f_{z}\Lambda_{z}$. Note that

[TABLE]

Applying integration by parts to the third term in Eq. (31) and summing up all the phonon modes gives

[TABLE]

where thermal conductivity $\kappa_{x}$ is given by

[TABLE]

and

[TABLE]