Relativistic stable processes in quasi-ballistic heat conduction in thin film semiconductors

TL;DR

This paper introduces a relativistic alpha stable process model to explain quasi-ballistic heat conduction in thin semiconductor films, connecting experimental ultrafast laser heating results with advanced stochastic process theory.

Contribution

It applies relativistic alpha stable processes and Feynman-Kac formalism to model heat conduction, bridging short-scale Levy processes and large-scale Brownian motion.

Findings

Model fits experimental ultrafast laser heating data

Derives sharp bounds for transition kernels using Feynman-Kac formalism

Establishes a theoretical connection between Levy processes and Brownian motion

Abstract

In this article, we show how relativistic alpha stable processes can be used to explain quasi-ballistic heat conduction in semiconductors. This is a method that can fit experimental results of ultrafast laser heating in alloys. It also provides a connection to a rich literature on Feynman-Kac formalism and random processes that transition from a stable L\'evy process on short time and length scales to the Brownian motion at larger scales. This transition was captured by a heuristic truncated L\'evy distribution in earlier papers. The rigorous Feynman-Kac approach is used to derive sharp bounds for the transition kernel. Future directions are briefly discussed.

Equations

$$T(t,x)=\mathbf{E}\left[T_0(x+X_t)\right],$$

$${\varphi }_M(\xi )\equiv \mathbb{E}\left[\exp \left(ı\xi \cdot (X_t^M-X_0^M)\right)\right]=\exp \left(-t\left({\left(|\xi {|}^2+M^{2/\alpha }\right)}^{\alpha /2}-M\right)\right).$$

$${\partial }_{t}T(t,x)={\mathcal{L}}^{M}T(t,x),\text{ with }T(0,x)=T_0(x).$$

$$\mathbb{P}\left(X_{s+t}^M\in A|X_s=x\right)={\int }_A{p}_t^M(x,y)dy,$$

$$p_t^M(x,y)={\mathcal{F}}^{-1}{\varphi }_M(x-y)=\frac{1}{(2\pi {)}^d}{\int }_{{\mathbb{R}}^d}{e}^{ı(x-y)\cdot \xi }{e}^{-t\{(M^{2/\alpha }+|\xi {|}^2{)}^{\alpha /2}-M\}}d\xi .$$

$$p_t^M(x,y)\le c_1(M^{d/\alpha -d/2}{t}^{-d/2}+t^{-d/\alpha }).$$

$$z^{\beta }-\beta z\le 1-\beta .$$

$$\frac{M}{(|\xi {|}^2+M^{2/\alpha }{)}^{\alpha /2}}-\frac{\alpha M^{2/\alpha }}{2(|\xi {|}^2+M^{2/\alpha })}\le 1-\frac{\alpha }{2},$$

$${\left(|\xi {|}^2+M^{2/\alpha }\right)}^{\alpha /2}-M\ge \frac{\alpha |\xi {|}^2}{2(|\xi {|}^2+M^{2/\alpha }{)}^{1-\alpha /2}}.$$

$$A_1=\{\xi :|\xi {|}^2\le M^{2/\alpha }\},\text{and}{A}_2=\{\xi :|\xi {|}^2>M^{2/\alpha }\}.$$

$$(|\xi {|}^2+M^{2/\alpha }{)}^{\alpha /2}-M\ge \{\begin{array}{ll}\frac{\alpha }{2(2M^{2/\alpha }{)}^{1-\alpha /2}}|\xi {|}^2,& \text{when }\xi \in A_1\\ \frac{\alpha }{2^{2-\alpha /2}}|\xi {|}^{\alpha },& \text{when }\xi \in A_2.\end{array}$$

$$p_s^M(x,y)=\frac{1}{(2\pi {)}^d}{\int }_{{\mathbb{R}}^d}{e}^{ı(x-y)\cdot \xi }{e}^{-t\{(M^{2/\alpha }+|\xi {|}^2{)}^{\alpha /2}-M\}}d\xi .$$

$$\begin{gathered} p_t^M(x,y)\le \frac{1}{(2\pi {)}^d}{\int }_{A_1}{e}^{-\frac{\alpha t}{2(2M^{2/\alpha }{)}^{1-\alpha /2}}|\xi {|}^2}d\xi \\ +\frac{1}{(2\pi {)}^d}{\int }_{A_2}{e}^{-\frac{\alpha t}{2^{2-\alpha /2}}|\xi {|}^{\alpha }}d\xi \le I_t^1+I_t^2, \end{gathered}$$

$$I_t^1$$

$$I_t^2$$

$$I_t^1=\frac{c_2}{t^{d/2}},\text{and}{I}_t^2=\frac{c_3}{t^{d/\alpha }}.$$

$$C_1^{-1}\left(t^{-d/\alpha }\wedge \frac{t{e}^{-M^{1/\alpha }|x-y|}}{|x-y{|}^{(d+\alpha +1)/2}}\right)\le p_t^M(x,y)\le C_1\left(t^{-d/\alpha }\wedge \frac{t{e}^{-M^{1/\alpha }|x-y|}}{|x-y{|}^{(d+\alpha +1)/2}}\right).$$

$$C_2^{-1}{t}^{-d/2}\le p_t^M(x,y)\le C_2{t}^{-d/2}.$$

$$P(|\xi |,t)=\exp \left(-t\psi (|\xi |)\right),$$

$$\psi (|\xi |)=D_{\alpha }\left[{\left(|\xi {|}^2+M^{2/\alpha }\right)}^{\alpha /2}-M\right].$$

$$x_0=|{\xi }_0{|}^{-1}=M^{-1/\alpha }.$$

$$\psi (|\xi |)$$

$$=D_{\alpha }|{\xi }_0{|}^{\alpha }\left[{\left({\tilde{\xi }}^2+1\right)}^{\alpha /2}-1\right],$$

$$\text{alpha-stable regime }\tilde{\xi }\gg 1$$

$$\text{Brownian regime }\tilde{\xi }\ll 1$$

$$\frac{\alpha D_{\alpha }}{2|{\xi }_0{|}^{2-\alpha }}=D⟹D_{\alpha }=\frac{2D}{\alpha x_0^{2-\alpha }}.$$

$$\psi (|\xi |)=\frac{2D}{\alpha x_0^2}\left[{\left(1+x_0^2|\xi {|}^2\right)}^{\alpha /2}-1\right].$$

Full text

S. Tindel is supported by the NSF grant DMS-1613163

1. Introduction

Standard heat flow, at a macroscopic level, is modeled by the random erratic movements of Brownian motions starting at the source of heat. However, this diffusive nature of heat flow predicted by Brownian motion is not seen in certain materials (semiconductors, dielectric solids) over short length and time scales [18]. Experimental data portraying the non-diffusive behavior of heat flow has been observed for transient thermal grating (TTG) [12, 13], time domain thermoreflectance (TDTR) [26] and others [9, 22], by altering the physical size of the heat source. The thermal transport in such materials is more akin to a superdiffusive heat flow, and necessitates the need for processes beyond Brownian motion to capture this heavy tail phenomenon. Recent works [7, 8, 10, 11, 14, 17, 19] try to explain the physics behind the quasiballistic heat dynamics. But these methods, driven mostly by the Boltzmann transport equation, are infeasible for processing experimental data. Some more recent studies [15, 16] try to explain the non-diffusive heat flow through hyperbolic diffusion equations, however, closer investigation shows that these methods fail to capture the inherent onset of nondiffusive dynamics at short length scales in periodic heating regimes.

The attempts mentioned above fail to provide a stochastic process that would explain the heat dynamics under short length-time regimes. The most natural stochastic process to explain a superdiffusive behavior is an alpha-stable Lévy process [1]. Alpha-stable Lévy processes differ from Brownian motion in that its movements are governed by stable distributions as compared to Gaussian distributions for the latter. In this context, some of us have tried to explain the heat flow dynamics through a "truncated Lévy distribution" approach [24, 25], where it has been possible to extract the value of the Lévy superdiffusive coefficient $\alpha$ that regulates the alloy’s quasiballistic heat dynamics.

The current contribution can be seen as a further step in this direction. Specifically, let $T(t,x)$ designate the temperature of a semiconductor or dielectric solid in the experimental settings alluded to above, with initial condition $T_{0}(x)$. Then we shall describe $T$ through the following Feynman-Kac formula (see Section 2.1 for more details about Feynman-Kac representations):

[TABLE]

where $X$ is a well-defined Lévy process that captures the observed quasiballistic heat dynamics, in addition to being a good candidate for explaining the usual diffusive nature under non-special large length-time regimes. We shall see that such a process $X$ can be chosen as a so-called relativistic stable process (see [20], and [3] for properties related to the relativistic Schrödinger operator). It possesses the remarkable property of behaving like an alpha-stable process under short length-time scales while being closer to Brownian motion otherwise. This is reflected in the estimates of the transition density, provided below in Section 2.2. Summarizing, our result lays the mathematical foundations of heat flow modeling on short time scales by means of stochastic processes. In addition, in spite of the fact that our computations are mostly one-dimensional, the model we propose allows natural generalizations to multidimensional and multilayer settings.

2. Relativistic stable process: a primer

In this section we give a short introduction on relativistic stable processes. We first recall the definition of this family of processes. Then we will give some kernel bounds indicating how relativistic processes transition, as $t$ increases, from an $\alpha$-stable behavior to a Brownian type behavior (this property being crucial to model quasiballistic heat dynamics in semiconductors).

2.1. Characteristics and Feynman-Kac formula

For any $M>0$ and $\alpha\in(0,2]$, a relativistic $\alpha$-stable process $X^{M}$ on $\mathbb{R}^{d}$ with a positive weight $M$ is a Lévy process whose characteristic function is given, for any $t\geq 0$ and $\xi\in\mathbb{R}^{d}$, by:

[TABLE]

When $M=0$, $X^{M}$ is simply a rotationally symmetric $\alpha$-stable process in $\mathbb{R}^{d}$. When $\alpha=2$, regardless of the value of $M$, the process is a Brownian motion. The infinitesimal generator of $X^{M}$ is $\mathcal{L}^{M}=M-(-\Delta+M^{2/\alpha})^{\alpha/2}$. When $\alpha=1$, this reduces to the free relativistic Hamiltonian $M-\sqrt{-\Delta+M^{2}}$, which explains the name of the process. An explicit expression for the Lévy measure of $X^{M}$ can be found in [5, 6, 4]. We omit this formula for sake of conciseness, since it will not be used in the remainder of the paper.

Lévy processes like $X^{M}$ are classically used in order to represent solutions of deterministic PDEs. In our case, consider the following equation governing the temperature $T$ in our material:

[TABLE]

Then it is a well known fact (see [1]) that the solution $T$ to (3) can be represented by the Feynman-Kac formula (1), where the process $X^{M}$ is our relativistic $\alpha$-stable process. The Feynman-Kac representation is crucial in order to get equation (13) below.

2.2. Transition kernel estimates

In this subsection, we identify the behavior of a relativistic stable process with a stable process on short time scales and a Brownian motion on larger time scales. As mentioned above, this will be achieved by observing the patterns exhibited by the transition kernel of $X_{t}^{m}$. Some results will be stated without formal proof, and interested readers are referred to [5, 6, 4] for more details.

Since $X^{M}$ is a Levy process, it is also a Markov process. As such it admits a transition kernel $p_{t}^{M}$, defined by:

[TABLE]

for all $x\in\mathbb{R}^{d}$ and $A\subset\mathbb{R}^{d}$. Notice that $p^{M}$ is related to the function $\phi_{M}$ (see definition (2)) as follows:

[TABLE]

We start by a simple bound on $p^{M}_{t}$, exhibiting the stable behavior for small times and the Brownian behavior for large times. Observe that this bound does not depend on the space variables $x,y$. We include its proof, which is based on elementary considerations involving Fourier transforms, for sake of completeness.

Theorem 2.1**.**

Consider a relativistic $\alpha$-stable process $X^{M}$, and let $p^{M}$ be its transition kernel. Then there exists $c_{1}=c_{1}(\alpha)>0$ such that for all $t>0$ and all $x,y\in\mathbb{R}^{d}$:

[TABLE]

Proof.

The strategy of our proof is based on the fact that the characteristic function $\phi_{M}$ defined by (2) behaves like a Gaussian characteristic function for low frequencies, and like an $\alpha$-stable characteristic function for high frequencies. We shall quantify this statement below.

Step 1: Elementary inequalities: Let $\beta=\frac{\alpha}{2}$. The following inequality, valid for for $0\leq z\leq 1$, is readily checked:

[TABLE]

Substituting $z=\frac{M^{2/\alpha}}{|\xi|^{2}+M^{2/\alpha}}$ in (6), we thus have:

[TABLE]

which yields:

[TABLE]

Relation (7) prompts us to split the frequency domain in two sets:

[TABLE]

Accordingly, we get the following lower bounds:

[TABLE]

This relation summarizes the separation between an $\alpha$-stable and a Gaussian regime alluded to above.

Step 2: Consequence for the transition kernel. Recall relation (4) for $p^{M}$, that is:

[TABLE]

In the integral above, we simply bound $|e^{\imath(x-y)\cdot\xi}|$ by 1 and split the integration domain $\mathbb{R}^{d}$ into $A_{1}\cup A_{2}$. Taking our relation (7) into account, this yields:

[TABLE]

where

[TABLE]

It is now easily checked that

[TABLE]

Plugging this information into (10), our claim (5) follows. ∎

The upper bound (5) already captures a lot of the information we need on relativistic stable processes. Invoking sophisticated arguments based on stopping times and Dirichlet forms, one can get upper and lower bounds on the transition kernel $p^{M}$ involving some exponential decay in the space variables $x,y$. We summarize those refinements in the following theorem.

Theorem 2.2**.**

Let $p^{M}$ be the transition kernel defined by (4). Then the following estimates hold true.

(i) Small time estimates.* Let $T>0$ be a fixed time horizon. Then there exists $C_{1}>0$ such that for all $t\in(0,T]$ and $x,y\in\mathbb{R}^{d}$,*

[TABLE]

(ii) Large time estimates.* There exists $C_{2}\geq 1$ such that for every $t\in[1,\infty)$ and $x,y\in\mathbb{R}^{d}$,*

[TABLE]

3. Application of relativistic stable processes to thermal modelling

In this section we show how to apply the mathematical formalism of Section 2 to our concrete physical setting. More specifically, in Section 3.1 we shall introduce length scales in our Lévy exponent (2). Then Section 3.2 is devoted to a description of our experimental setting, and also relates our measurements to the Fourier exponents we have put forward.

3.1. Formulation in terms of material thermal properties

The aforementioned evolution of relativistic processes, from alpha-stable behaviour at short length and time scales to regular Brownian motion at longer scales, renders them suitable to describing quasiballistic thermal transport in semiconductor alloys. Let us consider such a material, having nominal thermal diffusivity $D=\kappa/C_{v}$ with $\kappa$ being the thermal conductivity and $C_{v}$ the volumetric heat capacity (in ${\rm Jm}^{-3}{\rm K}^{-1}$). The physical quantity we have access to is a slight variation of the function $T(t,x)$ defined by (1). Specifically the single pulse response for the d-dimensional excess thermal energy can be expressed as $P(t,x)=C_{v}\,\Delta T(t,x)$. Under the Lévy flight paradigm the Fourier transform of $P$ is written as

[TABLE]

for a given Lévy exponent (also called special heat propagator) $\psi$. For the relativistic case under study here, this spatial heat propagator is simply a multiple of the function $\phi_{M}$ introduced in (2). Namely $\psi$ reads

[TABLE]

The prefactor $D_{\alpha}$ with unit mα/s denormalises the characteristic function for dimensionless space and time variables defined by Eq. (2) to its physical counterpart, and denotes the fractional diffusivity of the alpha-stable regime as we shall see shortly.

For thermal modelling purposes it is furthermore convenient to reformulate the process mass $M$, which has an exponent-dependent unit 1/mα, in terms of an associated characteristic length scale $x_{0}$ around which the transition from alpha-stable (Lévy) to Brownian dynamics takes place . Notice that according to our previous analysis leading up to Eq. (8) and (9) this length should be given by:

[TABLE]

This means that expression (14) can be recast as

[TABLE]

where $\tilde{\xi}\equiv|\xi|/|\xi_{0}|$. With those values of $M$ and $\xi_{0}$ in hand, we can translate (9) into an asymptotic transport limit as follows:

[TABLE]

The former corresponds to Lévy superdiffusion with characteristic exponent $\alpha$ and fractional diffusivity $D_{\alpha}$; the latter should recover to nominal diffusive transport $\psi(|\xi|)\equiv D|\xi|^{2}$. In order to make the last relation compatible with (17) we must set

[TABLE]

Finally, plugging (18) into (3.1) the heat propagator reads

[TABLE]

This formulation contains 3 material dependent parameters, each with an intuitive physical meaning: the characteristic exponent $\alpha$ of the alpha-stable regime; the nominal diffusivity $D$ of the Brownian regime; and the characteristic length scale $x_{0}$ around which the transition between those two asymptotic limits occurs (Fig. 1). In the sections that follow, we determine these parameter values for In0.53Ga0.47As by fitting a thermal model built upon the propagator (19) to time-domain thermoreflectance (TDTR) measurement signals.

3.2. Modelling of TDTR measurement signals

The central principle in TDTR is to heat up the sample with ultrashort pump laser pulses, and then monitor the thermal transient decay using a probe beam. Pulses from the laser are split into a pump beam and probe beam. The pump pulses pass through an electro-optic modulator (EOM) before being focused onto the sample surface through a microscope objective. A thin (50-100 nm) aluminium film is deposited onto the sample to act as measurement transducer: the metal efficiently absorbs the pump light and converts it to heat, and translates temperature variations to changes in surface reflectivity which can be captured by the probe. Lock-in detection at the pump modulation frequency $f_{\text{mod}}$ resolves the thermally induced reflectivity changes captured by the probe beam. A mechanical delay stage allows to vary the relative arrival time of the pump and probe pulses at the sample with picosecond resolution. To minimise the impact of random fluctuations in laser power and the variation of the pump beam induced by the delay stage, thermal characterisation is performed not on the raw lock-in signal itself but rather the ratio $-V_{\text{in}}/V_{\text{out}}$ of the in-phase and out-of-phase components as a function of the pump-probe delay.

Theoretical ratio curves $-V_{\text{in}}/V_{\text{out}}$ can be computed semi-analytically through mathematical manipulation of the semiconductor single-pulse response ((13)), as described in detail in Refs. [2, 21, 23]. Briefly, we first obtain the surface temperature response of a semi-infinite semiconductor to a cylindrically symmetric energy input via Fourier inversion of ((13)) with respect to the cross-plane coordinate. Next, a matrix formalism that accounts for heat flow in the metal transducer and across the intrinsic thermal resistivity $r_{\text{ms}}$ (in K-m2/W) of the metal-semiconductor interface provides the temperature response, weighted by the Gaussian probe beam, of the transducer top surface induced by a Gaussian pump pulse. Finally, harmonic assembly of this response at frequencies $n\cdot f_{\text{rep}}\pm f_{\text{mod}}\quad(\text{for }n=0,1,\ldots)$ accounting for the laser repetition rate $f_{\text{rep}}$, pump modulation frequency $f_{\text{mod}}$, and phase factors induced by the pump-probe delay $\tau$ yields the theoretical lock-in ratio signal $-V_{\text{in}}/V_{\text{out}}(\tau)$.

4. Experimental analysis

We have applied our model to TDTR measurements taken on a $\simeq$2 micron thick film of In0.53Ga0.47As ($C_{v}\simeq 1.55\,$MJ/m3-K) that was MBE-grown on a lattice-matched InP substrate. We note that although the semiconductor alloy under study (the InGaAs layer) is a geometrically thin film, thermally speaking it can still be considered, as is assumed by the thermal model, as a semi-infinite layer with good approximation. This is because the effective thermal penetration length $\ell=\sqrt{D/(\pi f_{mod})}$ stays firmly within the film over the experimentally probed modulation range $0.8\,\text{MHz}\lesssim f_{mod}\lesssim 18\,\text{MHz}$. The aluminium transducer deposited onto the sample measured 64 nm in thickness as determined by picoseconds accoustics. We used pump and probe beams with $1/e^{2}$ radii at the focal plane of 6.5 and 9 microns respectively, with respective powers of 17 and 8 mW at the sample surface.

In the thermal model with relativistic stable heat propagator (19), we fixed the heat capacity at the aforementioned 1.55 MJ/m3-K. Theoretical ratio curves were then collectively fitted through nonlinear least-square optimisation to signals measured at 7 different modulation frequencies to identify the 4 key thermal parameters: the characteristic exponent $\alpha$ of the Lévy superdiffusion regime; the quasiballistic-diffusive transition length scale $x_{0}$ associated to the mass $M$; the nominal thermal conductivity $\kappa=C_{v}D$ of the diffusive regime; and the thermal resistivity $r_{ms}$ of the transducer/semiconductor interface. The resulting best-fitting values $\alpha$ = 1.695 , $x_{0}$ = 0.86$\,\mu$m , $\kappa$ = 5.82 W/m-K , $r_{ms}$ = 4.28 nK-m2/W yield an excellent agreement with the measured signals (Fig. 2). Theoretical curves with parameter values deviating from the best fitting ones (Fig. 3) furthermore visually reveal the sensitivity to each of the parameters and illustrate the good quality of the best fit.

5. Conclusions and Outlook

Quasi-ballistic heat propagation in materials can be studied using atomic parameters through a multitude of techniques. First principle calculations and multi spectral phonon Boltzmann transport equations are very powerful in this regard. However, their use in the study of heat propagation in multi-layer/anisotropic materials and materials with complex geometries is limited. The Feynman-Kac representation of solutions to partial differential equations with non-local parameters can potentially provide alternative approaches to explain experimental thermal data. In this article we have replaced the traditional heat equation by a different PDE, whose solution has a Feynman-Kac representation driven by the so-called relativistic stable Lévy process. The transition characteristics of this process is in harmony to the heat propagation behaviour exhibited by TDTR data. In general, numerical approximations of the PDE solution can also be achieved through Monte Carlo simulations of the driving stochastic process in the Feynman-Kac formula. In particular, these numerical computations may provide substitute techniques to optimize materials or source geometry in order to reduce heating from nanoscale and/or ultrafast devices.

Our next challenge in this direction will be to model multidimensional transport in multilayer structures. To this aim, we shall investigate two methods: (i) Monte Carlo simulation according to our Feynman-Kac representation (1), taking into account jumps and change of media. (ii) Related PDEs involving the non local operator $\mathcal{L}^{M}=M-(-\Delta+M^{2/\alpha})^{\alpha/2}$, with boundary terms corresponding to the different layers. Both methods rely crucially on the relativistic Lévy representation advocated in this paper. They will be subject of future publications.