Heat Transport in Nanoscale Heterosystems : A Numerical and Analytical Study

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.


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.Electronelectron 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][4][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 picosecond and 10 picoseconds, the system can be described by the usual heat conduction equation: with ρ and c the specific mass and specific heat capacity, respectively.In (1), Fourier's law has been used relating the heat flux to the temperature-gradient. K is the thermal conductivity which generally is a tensor.The source term A( r, t) 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 A( r, t) and assuming isotropic media with constant thermal conductivity K [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][8][9].The heat flow Q across an interface is given by where R K denotes the thermal boundary resistance, which couples the heat flow across the interface to the temperature jump ΔT 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][13][14][15], heat transport in superlattices [16][17][18], ballistic transport [19], exact phonon dispersion curves instead of the usually applied Debyeapproximation [20], and electronic contributions to the heat transport across the interface [21][22][23].Recently, the Kapitzaeffect, 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][26][27][28][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 fslaser pulse excitation [32][33][34][35][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 electronphonon coupling time can even be obtained.This constitutes a new approach for studying such processes.

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 "z-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 K are reduced to a single, but z -dependent constant in this one-dimensional treatment.
At the surface the Neumann-type boundary condition, ∇T = 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 where α is the inverse absorption length and r the reflection coefficient for light at a given wavelength λ [4,6].q(t) is the temporal profile for the heat generation, that is, the temporal profile of the laser pulses.For the simulation, q(t) has been taken constant for a period of 45 femtoseconds according to the laser pulse duration used in the experiment.I 0 is the integrated laser-pulse intensity which was set to 8.5 × 10 14 (a) at 300 K.
W/m 2 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).
At the interface the energy flux is given by (2) which is taken into account by applying as Neumann-type boundary conditions [25,26,37].In (4), T f and T s are the film-and substrate temperatures, respectively.K f and K s 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 R K = 9.76 × 10 −8 (K m 2 )/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 z = 10 nm.The maximum temperature of the Bifilm 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 Sisubstrate 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 Bifilm is given by (2).It has to be noted that the thickness of the silicon substrate in the simulation d 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.
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 (z = 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 Bifilm.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 filmthicknesses 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).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 Bisingle 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 with d denoting the film-thickness.ρ and c are the mass density and heat capacity of the thin film, respectively.For a constant substrate temperature T s , the above equation results in an exponential decay of the film temperature T f with a time constant: The validity of (6) will be further discussed in section 3.By determining the decay constant, the thermal boundary  resistance R K 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 R K determined from the decay constant of the 0.9 As evident from the figure, the extracted R K linearly depends on the spatial discretization Δz and reaches the input value only for Δz→0.In (b), the numerical integration is performed with a spatial step-width of Δz = 2 nm.The numerical integration uses the literature values from Table 1 and R K = 9.76 × 10 −8 (K m 2 )/W.The Bi-film-thickness is 10 nm.
numerical simulation critically depends on the size of the spatial discretization Δz. 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 Δz = 5 nm, cf. Figure 3(a)).Only for Δz→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 −4 for the extracted thermal boundary resistance as shown in Figure 3(b).For the simulation, we used a spatial and temporal discretization of Δz = 0.25 nm and Δt = 0.01 femtosecond, respectively, which results in a thermal boundary resistance that is 2% larger than the input value.

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 onedimensional heat diffusion equation after excitation: with κ = K/(ρc) the diffusivity.Introducing the dimensionless variables: where T i is the initial surface temperature (z = 0), T s the substrate temperature, and d the film-thickness.For the derivation of the analytical solution, we apply the same boundary conditions used in the numerical simulation which are ∂θ ξ, t ∂ξ ξ=0 = 0, ∂θ ξ, t ∂ξ ξ=1 = − σ K θ ξ, t .(10) Standard procedure for solving (7) yields the series expansion: where the Eigenvalues λ n are obtained by evaluating the transcendental equation: The coefficients E n in (11) are determined from the initial condition θ(ξ, t = 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 electronphonon 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  With this initial condition the coefficients E n are given by [43] E n = 2 sin λ n λ n + sin λ n cos λ n .(14) On the other hand, if electron-phonon coupling occurs on a time scale which is faster than energy dissipation by hot electrons across the film, the heat generation can be described by ( 3) and the initial condition will be an exponential temperature distribution given by with α = αd and α the absorption coefficient used before.With this temperature distribution, the coefficients E n are given by 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, t) of the analytical solution (11) including the first 200 expansion terms are shown in Figure 4(a) for σ K = 0.01297, which corresponds to a 10 nm thin Bi-film on an Si-substrate.For t > 1, the surface cooling follows an exponential behavior with the same decay constant regardless of the temperature distribution at t = 0.At small t, the temperature evolution strongly depends on the initial condition as can be seen in the close up view for 0 < t < 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 ( t < 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 The similarity of the exponential surface temperature decay for t > 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 t > 0.3 as evident from Figure 4(b).For an initially constant temperature distribution, θ λ1 even describes the surface cooling behavior for t < 0.3.
From (12), it is seen that the Eigenvalues λ n are only determined by σ K .The dependence of the Eigenvalues λ 1 to λ 4 are shown in Figure 4(c), which explains the dominant contribution of θ λ1 to the series expansion for t > 0.3.For σ K < 1, λ 1 is more than an order of magnitude smaller, and for σ K > 1, λ 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)): Figure 4(c) shows that within 5%, the Eigenvalue λ 1 is equal to σ K for σ K < 1 (thick solid line).Using ( 9) and ( 17), the decay constant of the surface cooling τ K is given by which is the same result as stated previously by ( 6).This means that for a given set of constants ρ, c, and R K , the decay constant linearly depends on the film-thickness as long as σ K < 1 is fulfilled.The upper limit for the film-thickness is with the Kapitza-length l K 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 l K = 770 nm using literature value for the thermal conductance K and the experimentally determined thermal boundary resistance R K .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 σ K > 1, the Eigenvalue λ 1 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 R K : The decay constant obtained by ( 20) is associated with the heat diffusion in the film itself.The condition σ K < 1 implies a lower limit of the decay constant τ K for the temperature decay which is determined by the heat transport across the interface (cf. ( 18)): 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 l K , 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.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 τ exp = (1205 ± 70) picoseconds.Convolving the numerical result with a rectangular function of 70 picoseconds width has no influence on the time constant which is τ sim = 1187 picoseconds for the 10.4 nm thin Bi-film.
As an example, for a 10-nm-thin Bi-film τ K 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 t > 0.3 corresponds to 5 picoseconds for a 10 nm thin Bi-film.

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 R K = (9.7 ± 0.6) × 10 −8 (K m 2 )/W using bulk values for c 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].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 upper limit 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 compares the temperature evolution of 50 nm thick Bi-film, σ K = 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 t > 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 in section 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 filmthicknesses are in progress.

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 filmthickness 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 Kapitzalength, the Bi/Si-heterosystem is an ideal candidate for such studies.

Figure 1 :
Figure 1: Result of the numerical integration of the heat diffusion equation using literature values and R K = 9.76 × 10 −8 (K m 2 )/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.

Figure 2 :
Figure2: 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 × 10 −8 (K m 2 )/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.

Figure 3 :
Figure3: 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 Δz and (b) temporal discretization Δt.The temporal discretization in (a) is Δt = 0.01 femtosecond.As evident from the figure, the extracted R K linearly depends on the spatial discretization Δz and reaches the input value only for Δz→0.In (b), the numerical integration is performed with a spatial step-width of Δz = 2 nm.The numerical integration uses the literature values from Table1and R K = 9.76 × 10 −8 (K m 2 )/W.The Bi-film-thickness is 10 nm.

Figure 4 :
Figure 4: (a) Surface temperature evolution θ(ξ = 0, t = κt/d 2 ) including the first 200 terms of the analytical solution of the heat diffusion equation for σ 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 < t < 1.In (b), the first expansion terms θ λ1 corresponding to λ 1 are separately displayed.(c) Doublelogarithmic plot of Eigenvalues λ 1 to λ 4 versus σ K .The solid and dotted line represent σ K and π/2, respectively.

Figure 5 :
Figure 5: Comparison of the simulation result and experimental data for a 10.4 nm thin Bi-film deposited on a 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 τ exp = (1205 ± 70) picoseconds.Convolving the numerical result with a rectangular function of 70 picoseconds width has no influence on the time constant which is τ sim = 1187 picoseconds for the 10.4 nm thin Bi-film.

Figure 6 :
Figure6: Surface temperature evolution θ(0, t) for a 50 nm thick Bi-film ( σ 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 τ = 5823 picoseconds.For delays larger than 200 picoseconds, both models yield the same exponential decay.

Table 1 :
[40]rature values for the material constants used in the numerical simulation.r is the reflectivity for light of wavelength 800 nm, α is the absorption coefficient, K is the heat conductivity, ρ and c are the mass density and specific heat, respectively[40].