Heat Transport in Nanoscale Heterosystems: A Numerical and Analytical Study

Abstract

The numerical integration of the heat diffusion equation applied to the Bi/Si-heterosystem is presented for times larger than the characteristic time of electron-phonon coupling. By comparing the numerical results to experimental data, it is shown that the thermal boundary resistance of the interface can be directly determined from the characteristic decay time of the observed surface cooling, and an elaborate simulation of the temporal surface temperature evolution can be omitted. Additionally, the numerical solution shows that the substrate temperature only negligibly varies with time and can be considered constant. In this case, an analytical solution can be found. A thorough examination of the analytical solution shows that the surface cooling behavior strongly depends on the initial temperature distribution which can be used to study energy transport properties at short delays after the excitation.

1. Introduction

The ultrafast dynamics in a solid subsequent to short pulse laser irradiation are of great interest yielding insights into the microscopic processes of the energy transfer between different degrees of freedom of a solid. After excitation, a number of relaxation processes lead to the establishment of a thermal equilibrium in the laser-irradiated solid. Electron-electron thermalization is a rather fast process taking place on a subpicosecond timescale [1, 2]. Compared to this, the heat exchange between the electron and phonon subsystems is a slow process [3–5]. The nonequilibrium between these two subsystems can be described by the two-temperature model [3, 4].

For timescales larger than the electron-phonon coupling time, which typically lies between 1 picoseconds and 10 picoseconds, the system can be described by the usual heat conduction equation:𝜌𝑐𝜕𝑇(^→𝑟,𝑡)𝜕𝑡=∇⋅[𝐾∇𝑇(^→𝑟,𝑡)]+𝐴(^→𝑟,𝑡) (1)

with 𝜌 and 𝑐 the specific mass and specific heat capacity, respectively. In (1), Fourier's law has been used relating the heat flux to the temperature-gradient. 𝐾 is the thermal conductivity which generally is a tensor. The source term 𝐴(^→𝑟,𝑡) represents the heat generation in the solid per unit time and unit volume subsuming all microscopic steps from the absorption of a photon to the formation of the thermal equilibrium between electrons and phonons. Equation (1) can be solved analytically for simple forms of 𝐴(^→𝑟,𝑡) and assuming isotropic media with constant thermal conductivity 𝐾 [6].

For heterostructures, the assumption of an isotropic medium is no longer valid and (1) has to be solved independently for the different materials applying appropriate boundary conditions. It is well known that the boundary between two materials acts as a barrier to thermal heat diffusion [7–9]. The heat flow 𝑄̇ across an interface is given by𝑄̇=1𝑅_𝐾⋅Δ𝑇, (2)

where 𝑅_𝐾 denotes the thermal boundary resistance, which couples the heat flow across the interface to the temperature jump Δ𝑇 at the interface. If this heat transport is determined by phonons only, the temperature jump arises because phonons incident on the interface are partially reflected and only a fraction is transmitted across the interface. The simplest models that can be used to calculate the phonon transmission probabilities are the acoustic mismatch model (AMM) and the diffusive mismatch model (DMM) [7, 8]. If the wavelength of the phonons is larger than the interface roughness, the AMM can be applied for the calculation of the transmission coefficient [10]. The AMM treats phonons as acoustic waves that are reflected and refracted at the interface. The transmission probability is calculated by applying the acoustic analog of the Fresnel equations in optics [11]. In the framework of the DMM, which applies if the phonon wavelength is comparable to or smaller than the interface roughness, leading to a strong diffuse scattering at the interface, the transmission probability depends on the phonon density of states of the two adjacent media [7, 8, 10]. In the past, the above mentioned models have been extended and refined including lattice-dynamical calculations [12–15], heat transport in superlattices [16–18], ballistic transport [19], exact phonon dispersion curves instead of the usually applied Debye-approximation [20], and electronic contributions to the heat transport across the interface [21–23]. Recently, the Kapitza-effect, that is, the formation of a temperature jump at an interface, has also been observed in molecular dynamics simulations [24].

The most simple heterosystem which can be studied is a layered structure of two materials, that is, a thin film on a substrate. The thermal boundary resistance of such heterosystems for a large number of different material combinations and even the heat transport properties of nanolaminates has been determined in the past using the time-domain thermoreflectance technique [8, 25–29]. In a previous study, the temporal surface temperature evolution of a thin Bismuth-film on a Silicon substrate has been investigated by means of ultrafast electron diffraction [30, 31]. In this technique, a short electron pulse is diffracted at the surface with variable delays from an initial fs-laser pulse excitation [32–36]. The surface temperature is extracted from the diffraction spot intensity, which is affected by the Debye-Waller-effect. For the Bi/Si-system, it was observed that the initial temperature increase is followed by an exponential surface temperature decay. From the decay time constant, the thermal boundary resistance has been extracted, which was found to be in good agreement with values calculated from the AMM and DMM [31]. A direct determination of the thermal boundary resistance from the exponential decay time constant, however, is only possible if the substrate temperature is constant during the experiment [8]. Generally, (1) has to be solved separately for the two materials with a thermal boundary resistance as an input parameter yielding the best agreement between experiment and simulation [8, 25, 26].

In this paper, we will present numerical simulations for the temporal evolution of the surface temperature of Bi films of various thicknesses on a Silicon substrate by applying (1). It will be shown that for the Bi/Si-system, the thermal boundary resistance can be directly extracted from the experimentally observed decay constant without an elaborate comparison of simulations and experimental data. In addition, an analytical solution for the heat transport in the Bi/Si-system will be presented and discussed. It will be shown that the analytical approach well describes the results of the numerical simulation and observed surface cooling behavior. Although the discussion presented here is carried out for the Bi/Si-system, the results are applicable to other material combinations. That the temporal substrate temperature changes need to be small compared to the temperature changes in the film is the only restriction. If this condition is fulfilled, even information on the heat transfer mechanisms at times shorter than the electron-phonon coupling time can even be obtained. This constitutes a new approach for studying such processes.

2. Numerical Study

The numerical integration of (1) is accomplished in a standard, explicit FTCS-scheme (forward time centered space) in one dimension, namely, the cross-plane, which in the following will be referred to as the “𝑧-direction" [37]. The spatial and temporal discretization are both taken at finite and constant values obeying the stability criterion for the numerical integration [37]. In the above mentioned ultrafast electron diffraction experiment, the probed area (300 𝜇m × 4 mm) was an order of magnitude smaller than the excited area (spot size ≈4 mm) resulting in a laterally homogeneous heating of the probed area. The lateral heat-diffusion is, therefore, neglected in the following. Additionally, the tensor properties of the heat diffusion constant 𝐾 are reduced to a single, but 𝑧 -dependent constant in this one-dimensional treatment.

At the surface the Neumann-type boundary condition, ∇𝑇=0 has been applied. This condition is valid if no particles, that is, electrons, atoms, or clusters are leaving the surface and radiation losses can be neglected. In an experiment, both conditions can be fulfilled using low laser excitation energies. In the studied Bi/Si-systems, the maximum surface temperature was below 300 K, which means that radiative losses are negligible. Additionally, no ablation has been observed [30, 31].

At the second boundary, the back side of the simulated volume, a constant temperature has been applied (Dirichlet-condition), which can be considered as a heat sink (thermostat). The effect of this boundary condition on the temporal surface temperature evolution has been minimized by choosing a sufficiently large simulation slab.

The source term in (1) for one dimension can be written as 𝐴(𝑧,𝑡)=𝐼₀(1−𝑟)𝛼exp(−𝛼𝑧)𝑞(𝑡), (3)

where 𝛼 is the inverse absorption length and 𝑟 the reflection coefficient for light at a given wavelength 𝜆 [4, 6]. 𝑞(𝑡) is the temporal profile for the heat generation, that is, the temporal profile of the laser pulses. For the simulation, 𝑞(𝑡) has been taken constant for a period of 45 femtoseconds according to the laser pulse duration used in the experiment. 𝐼₀ is the integrated laser-pulse intensity which was set to 8.5×10¹⁴ W/m² in the simulation depicting the experimental value. Note, the above treatment of the source term is only correct for strong electron-phonon coupling resulting in a thermal equilibrium of electrons and phonons on a short timescale. Generally, the two-temperature model has to be applied in order to obtain the time-dependent temperature distributions of the electron and phonon subsystems [3, 4, 26]. In section 3, the effect of strong and weak electron-phonon coupling on the temporal temperature evolution will be discussed in more detail. In addition to the application of the two-temperature model, multiple reflection of light in the thin film and internal reflection, that is, the dynamic change of the optical properties during the excitation process, have to be considered for the initial heating dynamics [38, 39]. As the main concern of this work is the cooling behavior of thin Bi-films, we assume the source term (3) to be valid and use bulk values for the optical constants.

According to (3), the excitation of the Bi/Si-system with light of wavelength 800 nm results in a temperature jump at the interface because of the largely different absorption coefficients of the two materials (cf. Table 1).

Table 1 

Literature values for the material constants used in the numerical simulation. 𝑟 is the reflectivity for light of wavelength 800 nm, 𝛼 is the absorption coefficient, 𝐾 is the heat conductivity, 𝜌 and 𝑐 are the mass density and specific heat, respectively [40].

At the interface the energy flux is given by (2) which is taken into account by applying𝐾_𝑓𝜕𝑇_𝑓|𝜕𝑧_𝑧=𝑑=𝐾_𝑠𝜕𝑇_𝑠|𝜕𝑧_𝑧=𝑑1=−𝑅_𝐾(𝑇_𝑓−𝑇_𝑠), (4) as Neumann-type boundary conditions [25, 26, 37]. In (4), 𝑇_𝑓 and 𝑇_𝑠 are the film- and substrate temperatures, respectively. 𝐾_𝑓 and 𝐾_𝑠 are the heat diffusion constants in the film and the substrate. The boundary condition (4) guarantees that no heat is accumulated at the interface: the heat flow toward the interface equals the heat flow from the interface into the substrate. Due to the existence of a Schottky-barrier between Bi and Si and the small density of states in the interval up to 1.55 eV (energy of the photons with wavelength 800 nm) above the valence band maximum in Silicon, the energy transport across the interface by electrons is neglected [40, 41].

The numerical integration uses literature values for the material parameters that are tabulated in Table 1. It should be noted that the material parameters can depend on the dimension of the material, for example, the film-thickness. As their thickness dependencies are unknown, we use bulk values for Bismuth. The thermal boundary resistance is set to 𝑅_𝐾=9.76×10⁻⁸ (K m²)/W as experimentally determined [30, 31]. For all numerical integrations, a starting temperature of 80 K is used.

In Figure 1, the results of the numerical integration are shown. Figure 1(a) displays the spatial and temporal temperature distributions of a 10 nm thin Bi-film. The total simulation slab dimension in spatial direction was 110 nm, but for a better representation only the first 30 nm are shown. The boundary between Bi and Si is clearly visible for all delays at 𝑧=10 nm. The maximum temperature of the Bi-film at short delays is 240 K. For larger delays, the Bi-film cools down due to heat transport across the interface. On the contrary, the Si-temperature is below 81 K for all times. As discussed, above the small heat generation in Silicon for short delays is explained by the three orders of magnitude smaller absorption coefficient for light of wavelength 800 nm compared to Bi. In addition, at larger delays heat transmitted through the interface is efficiently dissipated in the Si-substrate which has a thermal conductivity that is two orders of magnitude higher than Bi (cf.Table 1). The Silicon substrate acts as a thermostat and the energy loss in the Bi-film is given by (2). It has to be noted that the thickness of the silicon substrate in the simulation 𝑑_Si=100 nm is smaller than the phonon mean free path in silicon (𝜆^Si_Ph≈250 nm) and the validity of Fourier's law in (1) is not fulfilled. However, as shown above, the temperature gradient in the Si-substrate is very small. We found that larger substrate thicknesses has no effect on the temperature evolution of the Bi-film reversely justifying the application of Fourier's law.

(a)

(b)

Figure 1 

Result of the numerical integration of the heat diffusion equation using literature values and $R_K=9.76\times {10}^{-8}$ (K?m2)/W. (a) shows a temporal and spatial temperature of a 10?nm thin Bi-film on a Si-Substrate. (b) displays the temporal surface temperature (distance from surface $z=0$?nm) evolution for different Bi-film-thicknesses. For all numerical integrations, the thickness of the Si-slab was 100?nm, but for a better representation only the first 30?nm of the total integration slab are shown in (a). Note, the unusual, nonlinear color scale in (a) is for descriptive reasons.

For a comparison with the experiment, the temporal evolution of the temperature is obtained by taking slices of the integration slab at the surface layer (𝑧=0 nm). A compilation of the time-dependent surface temperature for different film-thicknesses is displayed in Figure 1(b). The surface temperature evolution of Bi-films with thicknesses below ∼100 nm can be divided into two contributions: a fast temperature decay at short delays followed by a slow decay for long delays. These different decay behaviors can be explained by studying the spatial temperature profiles for different delays.

Figure 2 shows the spatial temperature profiles at 1  picosecond and 6 picoseconds delay of a 10 nm thin Bi-film. At 1 picosecond, delay the temperature distribution across the film is inhomogeneous. Heat diffusion in the film drives the system into an equilibrated state with a homogeneous temperature distribution across the film. For a 10 nm thin Bi-film this state is reached after ∼6 picoseconds (Figure 2 dashed line). This equilibration-time for heat diffusion in the film is smaller than the characteristic time for diffusion obtained from (21) and from [25]. We attribute this discrepancy to the application of Fourier's law in (1). As the phonon mean free path in Bi is 𝜆^Bi_Ph≈13 nm, the thermal energy transport is overestimated for film-thicknesses below the phonon mean free path resulting in a faster surface temperature decay. For such thin films the Boltzmann transport equations have to be solved which is beyond the scope of this study. For comparison with the experiment, we assume that the heat diffusion can be treated using Fourier's law because the surface temperature dynamics for these short timescales is not accessible due to the limited experimental resolution (cf. Section 4).

Figure 2 

Spatial temperature profiles at 1?picosecond (solid line) and 6?picoseconds (dashed line) delay of a 10?nm thin Bi-film. The thermal boundary resistance in the simulation is $R_K=9.76\times {10}^{-8}$?(K?m2)/W. The total simulation slab thickness is 100?nm but for a better representation only the first 30?nm are shown. At 1?picosecond, a temperature gradient in the Bi-film is still present. After $~$6?picoseconds the temperature gradient has vanished and results in a homogeneous spatial temperature distribution.

With increasing film-thickness, the time required to establish a homogeneous temperature distribution increases. Concomitantly, the temperature level of the homogeneous state decreases as the amount of absorbed energy is distributed over a larger volume. For larger film-thicknesses, this results in temporal surface temperature evolutions which are similar to the surface temperature evolution of a Bismuth single crystal. In fact, the surface temperature evolution of the 500 nm thick Bi-film is the same as obtained for a Bi-single crystal in the displayed range. We conclude that the fast temperature decay at short delays is governed by heat diffusion in the film itself.

The evolution at delays after a homogeneous temperature distribution is formed, follows an exponential behavior. This is well understood in terms of the thermal boundary resistance. Rewriting (2) results in𝜌𝑐𝑑𝜕𝑇_𝑓1𝜕𝑡=−𝑅_𝐾(𝑇_𝑓−𝑇_𝑠) (5) with 𝑑 denoting the film-thickness. 𝜌 and 𝑐 are the mass density and heat capacity of the thin film, respectively. For a constant substrate temperature 𝑇_𝑠, the above equation results in an exponential decay of the film temperature 𝑇_𝑓 with a time constant:𝜏_𝐾=𝑐𝜌𝑅_𝐾𝑑. (6)

The validity of (6) will be further discussed insection 3 By determining the decay constant, the thermal boundary resistance 𝑅_𝐾 can be extracted from the exponential surface temperature decay which can also be used to crosscheck the accuracy of the numerical integration. It turns out that the value of 𝑅_𝐾 determined from the decay constant of the numerical simulation critically depends on the size of the spatial discretization Δ𝑧.Figure 3(a) shows a linear relationship between the extracted thermal boundary resistance and the spatial discretization. For finite step widths, the difference between the input value and the extracted value for the thermal boundary resistance can be up to 25% (for Δ𝑧=5 nm, cf. Figure 3(a)). Only for Δ𝑧→0 the extracted thermal boundary resistance is equal to the input value. In contrast, the temporal discretization has a relative influence on the order of only 10⁻⁴ for the extracted thermal boundary resistance as shown inFigure 3(b). For the simulation, we used a spatial and temporal discretization of Δ𝑧=0.25 nm and Δ𝑡=0.01 femtosecond, respectively, which results in a thermal boundary resistance that is 2% larger than the input value.

(a)

(b)

Figure 3 

Dependence of the thermal boundary resistance $R_K$ extracted from the decay of the surface temperature obtained from the numerical integration of (1) on (a) the spatial discretization $\mathrm{?}z$ and (b) temporal discretization $\mathrm{?}t$. The temporal discretization in (a) is $\mathrm{?}t=0.01$ femtosecond. As evident from the figure, the extracted $R_K$ linearly depends on the spatial discretization $\mathrm{?}z$ and reaches the input value only for $\mathrm{?}z?0$. In (b), the numerical integration is performed with a spatial step-width of $\mathrm{?}z=2$?nm. The numerical integration uses the literature values fromTable 1 and $R_K=9.76\times {10}^{-8}$?(K?m2)/W. The Bi-film-thickness is 10?nm.

3. Analytical Solution of the Heat Diffusion Equation after Excitation

One result of the above numerical simulation is the negligible temperature variation of the substrate temperature in the Bi/Si-heterosystem. For this case, that is, thin film on a substrate with an approximately constant temperature, an analytical solution can be derived. Consider the one-dimensional heat diffusion equation after excitation:𝜕𝜕𝜕𝑡𝑇(𝑧,𝑡)=𝜅²𝜕𝑧²𝑇(𝑧,𝑡) (7) with 𝜅=𝐾/(𝜌𝑐) the diffusivity. Introducing the dimensionless variables:𝜃(𝜉,𝑡^{0𝑥002𝑑𝑐})=𝑇(𝑧,𝑡)−𝑇_𝑠𝑇_𝑖−𝑇_𝑠, (8)𝑧𝜉=𝑑,𝑡^{0𝑥002𝑑𝑐}=𝜅𝑡𝑑²,𝜎^{0𝑥002𝑑𝑐}_𝐾=𝑑𝑅_𝐾𝐾, (9) where 𝑇_𝑖 is the initial surface temperature (𝑧=0), 𝑇_𝑠 the substrate temperature, and 𝑑 the film-thickness. For the derivation of the analytical solution, we apply the same boundary conditions used in the numerical simulation which are𝜕𝜃(𝜉,𝑡^{0𝑥002𝑑𝑐})|𝜕𝜉_𝜉=0=0,𝜕𝜃(𝜉,𝑡^{0𝑥002𝑑𝑐})|𝜕𝜉_𝜉=1=−𝜎^{0𝑥002𝑑𝑐}_𝐾𝜃(𝜉,𝑡^{0𝑥002𝑑𝑐}). (10)

Standard procedure for solving (7) yields the series expansion:𝜃(𝜉,𝑡^{0𝑥002𝑑𝑐})=^∞_𝑛=1𝐸_𝑛𝑒^{−𝜆2𝑛0𝑥002𝑑𝑐𝑡}cos(𝜆_𝑛𝜉), (11) where the Eigenvalues 𝜆_𝑛 are obtained by evaluating the transcendental equation:𝜆_𝑛tan(𝜆_𝑛)=𝜎^{0𝑥002𝑑𝑐}_𝐾. (12) The coefficients 𝐸_𝑛 in (11) are determined from the initial condition 𝜃(𝜉,𝑡^{0𝑥002𝑑𝑐}=0). Within the scope of this paper, the initial condition is the state after thermal equilibration between the electron- and phonon-subsystems. Depending on the timescale on which the electron-phonon thermalization occurs, two limiting cases can be distinguished yielding different initial conditions. For weak electron-phonon coupling, the energy transport at short timescales is mediated by ballistic electrons resulting in a fast dissipation of the excitation energy in the film [42]. Provided that these hot electrons are reflected at the interface and remain in the film, this results in a homogeneous distribution of the absorbed energy and the phonon temperature will be constant across the film

𝜃(𝜉,𝑡^{0𝑥002𝑑𝑐}=0)=1,for0<𝜉<1. (13) With this initial condition the coefficients, 𝐸_𝑛 are given by [43]𝐸_𝑛=2sin(𝜆_𝑛)𝜆_𝑛+sin(𝜆_𝑛)cos(𝜆_𝑛). (14) On the other hand, if electron-phonon coupling occurs on a time scale which is faster than energy dissipation by hot electron across the film, the heat generation can be described by (3) and the initial condition will be an exponential temperature distribution given by𝜃(𝜉,𝑡^{0𝑥002𝑑𝑐}=0)=𝑒^{−0𝑥002𝑑𝑐𝛼𝜉} (15) with ^0𝑥0𝛼^02𝑑𝑐=𝛼𝑑 and 𝛼 the absorption coefficient used before. With this temperature distribution, the coefficients 𝐸_𝑛 are given by𝐸_𝑛=2𝜆_𝑛(𝛼^{0𝑥002𝑑𝑐}(1−𝑒^{−0𝑥002𝑑𝑐𝛼}cos(𝜆_𝑛))+𝜆_𝑛𝑒^{−0𝑥002𝑑𝑐𝛼}sin(𝜆_𝑛))(𝛼^{0𝑥002𝑑𝑐2}+𝜆²_𝑛)(𝜆_𝑛+sin(𝜆_𝑛)cos(𝜆_𝑛)). (16) This condition leads to an analytical solution of the heat diffusion equation that is comparable to the numerical simulation shown in Figures 1 and 2.

The temporal surface temperature evolution 𝜃(𝜉=0,𝑡^{0𝑥002𝑑𝑐}) of the analytical solution (11) including the first 200 expansion terms are shown in Figure 4(a) for ^0𝑥0𝜎^02𝑑𝑐_𝐾=0.01297, which corresponds to a 10 nm thin Bi-film on an Si-substrate. For ^0𝑥0𝑡^02𝑑𝑐>1, the surface cooling follows an exponential behavior with the same decay constant regardless of the temperature distribution at ^0𝑥0𝑡^02𝑑𝑐=0. At small 𝑡,^{0𝑥002𝑑𝑐} the temperature evolution strongly depends on the initial condition as can be seen in the close up view for 0<𝑡^{0𝑥002𝑑𝑐}<1 shown in Figure 4(b). If the initial temperature distribution is exponential, the result from the analytical model is equal to the numerical simulation: a fast temperature decrease (^0𝑥0𝑡^02𝑑𝑐<0.3) is followed by a slow surface cooling. Using the same arguments as above, the fast temperature decrease is driven by heat diffusion in the film itself until a homogeneous temperature distribution across the film is established. Subsequently, the surface cooling is determined by the heat transport across the interface. Consequently, an initially already homogeneous temperature distribution directly results in the slow surface cooling (Figures 4(a) and 4(b) dotted line).

(a)

(b)

(c)

Figure 4 

(a) Surface temperature evolution $?(?=0,\tilde{t}=?t/d^2)$ including the first 200 terms of the analytical solution of the heat diffusion equation for ${\tilde{s}}_K=d/(R_{K}K)=0.01297$, corresponding to a 10?nm thin Bi-film on an Si-substrate with an initially exponential (solid line) and constant (dashed line) temperature distribution. (b) close up view of (a) for $0<\tilde{t}<1$. In (b), the first expansion terms ${?}_{{?}_1}$ corresponding to ${?}_1$ are separately displayed. (c) Double-logarithmic plot of Eigenvalues ${?}_1$ to ${?}_4$ versus ${\tilde{s}}_K$. The solid and dotted line represent $\sqrt{{\tilde{s}}_K}$ and $p/2,$ respectively.

The similarity of the exponential surface temperature decay for ^0𝑥0𝑡^02𝑑𝑐>1 regardless of the initial conditions is explained by the dominant contribution of first term 𝜃_𝜆1 of the series expansion (11). Independent of the initial condition 𝜃_𝜆1 describes the temporal evolution for ^0𝑥0𝑡^02𝑑𝑐>0.3 as evident from Figure 4(b). For an initially constant temperature distribution, 𝜃_𝜆1 even describes the surface cooling behavior for ^0𝑥0𝑡^02𝑑𝑐<0.3.

From (12), it is seen that the Eigenvalues 𝜆_𝑛 are only determined by ^0𝑥0𝜎^02𝑑𝑐_𝐾. The dependence of the Eigenvalues 𝜆₁ to 𝜆₄ are shown in Figure 4(c), which explains the dominant contribution of 𝜃_𝜆1 to the series expansion for ^0𝑥0𝑡^02𝑑𝑐>0.3. For ^0𝑥0𝜎^02𝑑𝑐_𝐾<1𝜆₁ is more than an order of magnitude smaller, and for ^0𝑥0𝜎^02𝑑𝑐_𝐾>1, 𝜆₁ is at least a factor of 3 smaller than the other Eigenvalues. On basis of this observation, the transient surface temperature evolution can be described with a single exponential decay characterized by a decay constant (cf. (11)):^{0𝑥002𝑑𝑐}1𝜏=𝜆²₁. (17) Figure 4(c) shows that within 5%, the Eigenvalue 𝜆₁ is equal to ^0𝑥0𝜎^02𝑑𝑐_𝐾 for ^0𝑥0𝜎^02𝑑𝑐_𝐾<1 (thick solid line). Using (9) and (17), the decay constant of the surface cooling 𝜏_𝐾 is given by𝜏_𝐾=^{0𝑥002𝑑𝑐}𝑑𝜏⋅²𝜅=_0𝑥001_2𝑑𝑐𝜎_𝐾⋅𝑑²𝐾/𝑐𝜌=𝑐𝜌𝑅_𝐾𝑑, (18) which is the same result as stated previously by (6). This means that for a given set of constants 𝜌, 𝑐, and 𝑅_𝐾, the decay constant linearly depends on the film-thickness as long as ^0𝑥0𝜎^02𝑑𝑐_𝐾<1 is fulfilled. The upper limit for the film-thickness is1>𝜎^{0𝑥002𝑑𝑐}_𝐾=𝑑𝑅_𝐾𝐾→𝑑<𝑅_𝐾𝐾=𝑙_𝐾 (19) with the Kapitza-length 𝑙_𝐾 defined previously [9, 44]. The Kapitza-length is a measure for the thermal resistance of the interface in terms of the thermal resistance of a perfect crystal. In the case of a thin Bi-film on an Si-substrate 𝑙_𝐾=770 nm using literature value for the thermal conductance 𝐾 and the experimentally determined thermal boundary resistance 𝑅_𝐾. This means that the temperature difference on the two sides of the interface is the same as the temperature difference of a stationary temperature profile between two sites in bulk-Bi that are 770 nm apart.

For ^0𝑥0𝜎^02𝑑𝑐_𝐾>1, the Eigenvalue 𝜆₁ asymptotically reaches the value 𝜋/2 (cf. Figure 4 thick dotted line). In this limit, the decay-constant quadratically depends on the film-thickness, but is independent of the thermal boundary resistance 𝑅_𝐾: 𝜏_Diff=4𝜋²⋅𝑑²𝜅=4𝜋²𝑐𝜌𝐾𝑑². (20) The decay constant obtained by (20) is associated with the heat diffusion in the film itself.

The condition ^0𝑥0𝜎^02𝑑𝑐_𝐾<1 implies a lower limit of the decay constant 𝜏_𝐾 for the temperature decay which is determined by the heat transport across the interface (cf. (18)):𝜏_𝐾>𝑐𝜌𝐾𝑑²>𝜏_Diff, (21)

which is the same result as previously stated [25]: the thermal boundary resistance can be extracted from the temperature decay if the decay constant associated with the heat transport across the interface is larger than the time constant of the heat diffusion in the film itself. By comparing (19) and (21), this means that the film-thickness has to be smaller than the Kapitza-length 𝑙_𝐾, otherwise the surface cooling behavior is determined by the heat diffusion in the film.

The transition from a surface cooling that is determined by the heat transport across the interface to a surface cooling which is driven by heat diffusion in the film itself can also be seen in the numerical results (Figure 1) at a film-thickness on the order of the Kapitza-length.

As an example, for a 10 nm thin Bi-film 𝜏_𝐾 must be larger than 15 picoseconds which is well below the time constant of the numerical simulation of 𝜏_sim=1187 picoseconds and the experimentally observed value 𝜏_exp=1205 picoseconds (cf.Section 4). The condition that the thermal boundary resistance can only be extracted from the surface temperature decay for delays with ^0𝑥0𝑡^02𝑑𝑐>0.3 corresponds to 5 picoseconds for a 10 nm thin Bi-film.

4. Comparison of Numerical Results and Experiment

In Figure 5, an example of the temporal surface temperature evolution of a 10.4 nm thin Bi-film deposited on an Si(001) substrate at 300 K is shown. The transient surface temperature evolution is obtained from the (01)-diffraction spot by means of ultrafast electron diffraction [30, 31]. An exponential fit to the data (dashed line in Figure 5) yields a time constant for the decay of 𝜏_exp=(1205±70) picoseconds. As discussed above, the thermal boundary resistance can be determined from the decay constant which yields 𝑅_𝐾=(9.7±0.6)×10⁻⁸ (K m²)/W using bulk values for 𝑐 and 𝜌 (cf.Table 1). Within the error this value is the same as previously determined from the surface temperature decay of a 5.5 nm thin Bi-film [30, 31].

Figure 5 

Comparison of the simulation result and experimental data for a 10.4?nm thin Bi-film deposited on an Si(001) substrate. The experimental data are obtained from the temporal evolution of the (01)-spot intensity that is converted to a temperature using the Debye-Waller-effect [30, 31]. An exponential fit to the experimental data yields the decay constant $t_{\text{exp}}=(1205\pm 70)$?picoseconds. Convolving the numerical result with a rectangular function of 70?picoseconds width has no influence on the time constant which is $t_{\text{sim}}=1187$?picoseconds for the 10.4?nm thin Bi-film.

The result of the numerical integration is also shown in Figure 5 (dash-dotted curve). In order to account for the finite temporal resolution of the experiment, the result of the numerical integration has been convoluted with a boxed shaped function. The temporal resolution in the experiment is limited by the velocity mismatch between the probing electrons at grazing incidence and the pumping laser pulses at normal incidence [30]. During the travel time of the electrons across the surface, the resulting measured temperature is an average over the travel time which is 70 picoseconds for the above shown experiment (electron energy: 7 keV, sample width: 3 mm, incident angle: 5^°). Convolution with a rectangular function assumes that any spot of the sample is probed with an equal number of electrons. Note, this procedure gives an upperlimit for the temporal resolution. If the electron distribution is inhomogeneous across the sample width, the temporal resolution is increased because the major part of the electrons is diffracted from a smaller area of the sample. The transient temperature evolution obtained from the convolution of the simulation with a box-shaped function is also displayed in Figure 5 (solid line). Apart from the region of the initial temperature increase at short delays, both sets are similar and the surface temperature decays with a time constant of 𝜏_sim=1187 picoseconds which agrees with the experimentally determined value 𝜏_exp.

The differences between the convoluted and original temperature profiles at small delays are two-fold. One is the linear increase of the surface temperature of the convoluted profile compared to the step-like increase of the unconvoluted profile. The other is, the fast decrease of the surface temperature for small delays, driven by heat diffusion in the film itself, is leveled out by the convolution. For thicker films, however, the fast temperature decay for small delays is still present even after convolution with a 70 picoseconds wide rectangular function which is demonstrated in Figure 6.

Figure 6 

Surface temperature evolution $?(0,t)$ for a 50?nm thick Bi-film (${\tilde{s}}_K=0.065$). The curves are obtained by convolving the results from the analytical solution with a rectangular function of 70?picoseconds width in order to account for the limited temporal resolution of the experiment (see text for detailed discussion). Compared to the fast decay at short times resulting from the model with an initially exponential temperature distribution, the surface of a homogeneously heated film cools down exponentially with a time-constant $t=5823$?picoseconds. For delays larger than 200?picoseconds, both models yield the same exponential decay.

Figure 6 compares the temperature evolution of 50 nm thick Bi-film, ^0𝑥0𝜎^02𝑑𝑐_𝐾=0.065, obtained from the analytical solution of the heat diffusion equation with the two different initial conditions, that is, exponential and constant initial temperature distribution. To account for the finite experimental temporal resolution, the two results have been convolved with a rectangular function of 70 picoseconds width. For delays larger than 200 picoseconds, corresponding to ^0𝑥0𝑡^02𝑑𝑐>0.53, both evolutions follow an exponential decay with time-constant 𝜏=5823 picoseconds. Evidently, the cooling behavior for delays below 200 picoseconds drastically depends on the initial condition. An initially constant temperature distribution results in the slow exponential decay similar to delays larger than 200 picoseconds. Compared to this slow cooling, the surface temperature drops much faster if the initial temperature distribution is exponential since heat diffusion in the film itself dominates the cooling behavior at these short delays.

From the discussion of the initial condition insection 3, the two different surface cooling behaviors can be used to gain information on the energy transport at short times after the excitation. If the electron and phonon subsystems are weakly coupled, the main transport mechanism is via hot electrons resulting in homogeneously heating of the phonon system. In the limit of an instantaneous electron to phonon energy transfer, the phonon subsystem is expected to be an exponential temperature distribution across the film. To study this property, further experiments with varying film-thicknesses are in progress.

5. Summary

In conclusion, we have applied the one-dimensional heat diffusion model to the heterosystem of a thin Bi-film on an Si-substrate. The numerical integration was carried out for timescales larger than the typical electron-phonon coupling times when the two subsystems are in thermal equilibrium. The Bi-surface temperature has been determined for different film-thicknesses. For film-thicknesses below 100 nm, the surface cooling is characterized by two different timescales. For short delays after the initial heat pulse, the surface cooling is dominated by thermal diffusion in the film resulting in a homogeneous temperature distribution across the film. Subsequently, the surface temperature evolution is determined by the heat transport across the interface resulting in a slow exponential temperature decay. For these films, the decay constant linearly depends on the film-thickness. For film-thicknesses larger than 100 nm, the surface cooling is virtually determined by the heat diffusion in the film itself.

Because of its thermal properties, the Silicon substrate temperature stays constant for all delays, which allows the derivation of an analytical solution for the heat diffusion equation. The examination of the analytical solution shows that the linear dependence of the decay-constant is valid up to film-thicknesses on the order of the Kapitza-length. For thicknesses larger than the Kapitza-length, it is shown that the decay-constant quadratically depends on the film-thickness and the thermal boundary resistance cannot be determined from surface temperature decay. Additionally, the analytical study yields insight into the initial dynamics of energy dissipation. Two limiting cases were studied: strong versus weak electron-phonon coupling. These two cases result in different initial temperature distributions across the film which have a strong effect on the surface cooling behavior. We have shown that for certain sets of parameters, the experimentally obtained temporal resolution is sufficient to resolve the initial surface cooling behavior and we motivate further experimental studies on the ultrashort dynamics of energy dissipation. Due to its large Kapitza-length, the Bi/Si-heterosystem is an ideal candidate for such studies.

Acknowledgments

The authors gratefully acknowledge valuable discussion with S. Thomae and E. Pehlke. financial support by the Deutsche Forschungsgemeinschaft through SFB 616 “Energy Dissipation at Surfaces" is gratefully acknowledged.