Abstract

The Fourier and the hyperbolic heat conduction equations were solved numerically to simulate a frequency-domain thermoreflectance (FDTR) experiment. Numerical solutions enable isolation of pump and probe laser spot size effects and use of realistic boundary conditions. The equations were solved in time domain and the phase lag between the temperature of the transducer (averaged over the probe laser spot) and the modulated pump laser signal was computed for a modulation frequency range of 200 kHz–200 MHz. Numerical calculations showed that extracted values of the thermal conductivity are sensitive to both the pump and probe laser spot sizes, while analytical solutions (based on Hankel transform) cannot isolate the two effects. However, for the same effective (combined) spot size, the two solutions are found to be in excellent agreement. If the substrate (computational domain) is sufficiently large, the far-field boundary conditions were found to have no effect on the computed phase lag. The interface conductance between the transducer and the substrate was found to have some effect on the extracted thermal conductivity. The hyperbolic heat conduction equation yielded almost the same results as the Fourier heat conduction equation for the particular case studied. The numerically extracted thermal conductivity value (best fit) for the silicon substrate considered in this study was found to be about 82–108 W/m/K, depending on the pump and probe laser spot sizes used.

1 Introduction

The ability to control and manipulate heat transport at the nanoscale is important for the advancement of thermal management strategies in electronic and optoelectronic devices. In the past decade, noncontact optical pump probe techniques based on thermoreflectance have been used extensively for the study of heat transport at very small time and length scales. The two most commonly used techniques are the time domain thermoreflectance (TDTR) technique and the frequency domain thermoreflectance (FDTR) technique. In TDTR, the sample is heated using a pulsed laser that is modulated, and the surface is probed using a time-delayed laser that measures the change in the reflectivity of the surface caused by the change in temperature of the sample [13]. In FDTR, the sample is, instead, heated using a modulated continuous wave pump laser beam resulting in surface temperature (reflectivity) oscillations, which are then monitored using a probe laser. The lag in phase between the pump and probe laser signals is recorded. In either method, the surface of the target material (whose thermal properties are sought) is covered with an ultrathin metallic layer, often referred to as the transducer.

Extraction of the thermal conductivity of the substrate from measured thermoreflectance data requires use of a thermal transport model. Since both TDTR and FDTR experiments measure surface temperature (or some indicator of it), while thermal conductivity is a volumetric property, the two quantities can only be related through a thermal transport model. Furthermore, the presence of the transducer complicates matters since the heat must now transfer through multiple layers and the imperfect contact between the two layers. The most common model used for this purpose is based on the solution of the Fourier heat conduction equation in frequency domain. This was brought to the limelight by Cahill [1] and has been used by the vast majority of researchers since then. Multiple layers are treated using the well-known Feldman algorithm [4]. Interfaces between layers are treated as artificial layers whose thermal properties are adjusted to reproduce measured interface conductance values.

The thermal conductivity extracted from FDTR experiments using such a Fourier law-based model has been found to change when the modulation frequency of the pump laser is changed [13,5]. This behavior is attributed to the fact that when the laser modulation frequency is high, the thermal penetration depth, which is inversely proportional to the square root of the modulation frequency, is small and can often be smaller than the mean free path of some of the energy-carrying phonons. As a result, some phonons hardly scatter. This results in so-called ballistic-diffusive transport or quasi-ballistic transport. In this regime of transport, the effective thermal conductivity has been found to be smaller than the bulk value—a phenomenon known as thermal conductivity suppression [57].

In an effort to capture ballistic effects and develop a model that can predict thermal conductivity suppression across all modulation frequencies, researchers have proposed various enhancements to Fourier law-based models [818]. One class of these models makes use of the hyperbolic heat conduction equation, which accounts for finite group velocity of the phonons. The Cattaneo equation, the Cattaneo–Vernotte model, and the dual phase lag model fall in this category [8,9]. Even though the hyperbolic heat conduction equation accommodates finite phonon speed, all phonons, regardless of their type and frequency, are assumed to have the same speed. In an effort to remove this deficiency, Ramu and Bowers [10] proposed a two-band model in which a cutoff frequency is used to classify the phonons into ballistic and diffusive phonons. The ballistic phonons are then treated by adding a higher order correction term to the Fourier law that is derived from the Boltzmann transport equation (BTE) for phonons. In a similar model, Ma [11] treated nondiffusive effects by introducing an additional term in the Fourier heat conduction equation that involves the characteristic ballistic heat transport length as an additional parameter to the thermal conductivity of the substrate. This characteristic length is an additional tuning parameter in this model. In the ballistic-diffusive model proposed by Chen [12], and later expanded to complex three-dimensional geometry by Mittal and Mazumder [13,14], the phonon intensity is split into a diffusive component and a ballistic component. The diffusive component, by virtue of being more or less isotropic, can be treated using the method of spherical harmonics, while the ballistic component is treated using a surface-to-surface exchange formulation. More recently, a model that introduces a hydrodynamic term in the Fourier heat conduction equation—analogous to the advective term in the Navier–Stokes equation—has been proposed to capture ballistic effects [17,18]. In principle, the multidimensional phonon BTE encapsulates the necessary phonon physics, and should be used as the model to interpret the measured data, as demonstrated by Ali and Mazumder [19] for TDTR experiments. All of the aforementioned approximate models were proposed and exercised to avoid full-fledged (and very time-consuming) solution to the phonon BTE, and continue to be used for the extraction of thermal conductivity from both TDTR and FDTR experimental data.

One advantage of using the Fourier law and the aforementioned transform method is the efficiency of the calculation process. Numerical solution of the Fourier heat conduction equation in time domain, although considerably simpler than solving the BTE, is still time-consuming. Such calculations have been conducted by Xing et al. [20] using the commercial finite-element solver comsol. The numerical calculations were motivated by the desire to address nonlinearity (due to temperature dependence of the thermal conductivity) and to assess the effect of external heat loss. It was found that heat loss plays an insignificant role unless the modulation frequency is very low. Braun and Hopkins [21] have conducted numerical calculations using the Fourier heat conduction equation in multilayered materials under continuous and pulsed laser energy input with the goal of investigating the thermal penetration depth and its dependence on various parameters. In this study, we reassess the efficacy of numerical solutions by solving the governing heat conduction equation in time domain. This allows us to assess the validity of some of the assumptions invoked in transform methods. Specifically, the objective is to answer the following questions: (1) is an effective laser spot size, defined as reff=rpump2+rprobe2, adequate for extracting the thermal conductivity, or do the individual values of rpump and rprobe have any bearing on the results? (2) What is the effect of heat loss from the system (effect of boundary conditions)? (3) What is the effect of the finite thickness of the substrate? (4) What is the effect of the interface conductance between the substrate and the transducer? In addition to the Fourier heat conduction equation, numerical solution of the hyperbolic heat conduction equation is also explored in an effort to capture quasi-ballistic effects.

2 Theory and Solution Method

2.1 Analytical Solution: Fourier Heat Conduction.

The analytical solution to the Fourier heat conduction equation with periodic heat flux boundary conditions is obtained by first transforming the original equation in time domain to an equation in frequency domain, followed by invoking a Fourier transform. Details of this procedure are provided elsewhere [22]. In this particular case, the frequency domain solution for a single layer is obtained by taking the Hankel transform of the response due to the laser beam since the sample is assumed to have cylindrical symmetry. For a semi-infinite solid being heated by a periodic heat flux operating at an angular frequency ω, the response of the surface is given by [1,22]
ΔT=2πA0G(p)exp[π2p2(rpump2+rprobe2)/2]pdp
(1)
where p is the Hankel transform parameter, rpump and rprobe are the 1/e2 radius of the pump and the probe beam, respectively, and A is the amplitude factor of the heat flux (due to the pump laser). The quantity ΔT denotes the weighted (over the radius of the probe laser beam) average temperature and has both a real and an imaginary part. The quantity G(p) in Eq. (1) is given by
G(p)=1k4π2p2+q2;q2=iωα
(2)

where k is the thermal conductivity of the solid (or substrate) and α its thermal diffusivity. It is clear from Eq. (1) that the effect of individual variation of the values of rpump and rprobe cannot be delineated in the analytically calculated phase lag; rather the effective radius, defined as reff=rpump2+rprobe2, is the only quantity that matters.

When the solution is extended to a layered structure using the Feldman algorithm [4], as proposed in [1], the layers are numbered starting from n =1 being the one closest to the surface being heated by the laser. The thickness, thermal conductivity, and thermal diffusivity of each layer are denoted by Ln, kn, and αn, respectively. The expression for G(k), shown in Eq. (2), is then replaced by
G(p)=1γ1(B1+B1B1++B1);γn=knun;un=4π2p2+qn2
(3)
The quantities B+ and B are the growth and decay exponents along the thickness. For an N-layered structure, the furthest layer only has decay (BN=1) and has no growth (BN+=0). The growth and decay exponents for the other layers are determined using the following recursive relation starting from the Nth layer
(Bn+Bn)=12γn[exp(unLn)00exp(unLn)][γn+γn+1γnγn+1γnγn+1γn+γn+1](Bn+1+Bn+1)
(4)

In this method, the interface between layers is treated as an artificial layer for which the thermal conductivity, thermal diffusivity, and thickness are chosen such that a known thermal conductance value (=kn/Ln) is matched, and the thermal mass is small. However, lateral (radial) conduction within the interface, which is unphysical, cannot be eliminated.

2.2 Numerical Solution: Fourier Heat Conduction.

The starting point of the numerical solution is the transient Fourier heat conduction equation. Under the axisymmetric assumption, and constant thermophysical properties, this equation may be written in cylindrical coordinates as [22]
ρcTt=k2T=k2Tz2+krr(rTr)
(5)
where z is in the axial or through-plane direction, and r is in the radial or in-plane direction, as shown in Fig. 1. The boundary conditions in the section where the laser is shining (rrL) is written as
kTz|(t,0,rrL)=qL(r)[1+sinωt]
(6a)
where qL(r) is the radially varying laser heat flux profile, and is assumed to be Gaussian in shape [1]. As noted in Eq. (6a), the laser flux also has a temporally varying component with modulation frequency ω. At the axis of symmetry, the boundary condition is written as
Tr|(t,z,0)=0
(6b)
Fig. 1
Schematic of the computational domain and boundary conditions
Fig. 1
Schematic of the computational domain and boundary conditions
Close modal
On the top surface outside of the laser spot and the external (side and bottom) surfaces of the substrate, a Newton cooling boundary condition was used, namely
Topsurface:kTz|(t,0,r>rL)=ht[T(t,0,r>rL)T]
(6c)
Sidesurface:kTr|(t,z,rS)=hs[T(t,z,rS)T]
(6d)
Bottomsurface:kTz|(t,zT+zS,r)=hb[T(t,zT+zS,r)T]
(6e)

where ht, hs, and hb are the heat transfer coefficients on the top, side, and bottom surfaces, respectively, and T is the ambient temperature.

For numerical solution of Eq. (5) subject to the boundary conditions in Eq. (6), the solution domain is first split into two regions: the transducer and the substrate. The transducer is very thin, and consequently, may be assumed to have negligible temperature gradient in the axial or z-direction. Thus, it may be split into a series of radial control volumes, as shown in Fig. 2. Applying the finite volume procedure [23] to these control volumes, we obtain
[kTzT(2rj+Δrj)Δrj+Δrj+1+kTzT(2rjΔrj)Δrj+Δrj1+ρTcTzT2rjΔrjΔt]TT,j[kTzT(2rj+Δrj)Δrj+Δrj+1]TT,j+1[kTzT(2rjΔrj)Δrj+Δrj1]TT,j1=ρTcTzT2rjΔrjΔtTT,jold+2qt,jrjΔrj2qC,jrjΔrj
(7)
where qt,j is the heat flux on the top surface of the jth control volume (or cell) of the transducer, and is given by either Eq. (6a) or Eq. (6c) depending on the location of the cell. The density, thermal conductivity, and specific heat capacity of the transducer are denoted by ρT, kT and cT, respectively, while zT denotes its thickness. The radius of the jth cell's center is denoted by rj, while the radial span (grid size) is denoted by Δrj. Equation (7) is derived using the backward Euler time advancement method [23], wherein TT,jold denotes the temperature of the jth cell of the transducer at the previous time-step. The heat flux through the bottom surface of the transducer or the interfacial contact (see Fig. 2) is denoted by qC,j, and is related to the contact conductance by the following relation
qC,j=GC(TT,jTS,j,top)
(8)
where GC is the contact conductance (in W/m2/K), and TS,j,top are the temperatures on the top surface of the substrate. Equation (8), in fact, serves to connect the transducer to the substrate, details of which will be discussed shortly. Substitution of Eq. (8) into Eq. (7), followed by rearrangement yields
[kTzT(2rj+Δrj)Δrj+Δrj+1+kTzT(2rjΔrj)Δrj+Δrj1+ρTcTzT2rjΔrjΔt+2GCrjΔrj]TT,j[kTzT(2rj+Δrj)Δrj+Δrj+1]TT,j+1[kTzT(2rjΔrj)Δrj+Δrj1]TT,j1=ρTcTzT2rjΔrjΔtTT,jold+2qt,jrjΔrj+2GCTS,j,toprjΔrj
(9)
Fig. 2
Schematic representation of the stretched mesh used for computations in the two different regions (transducer on top and substrate at the bottom) and their coupling
Fig. 2
Schematic representation of the stretched mesh used for computations in the two different regions (transducer on top and substrate at the bottom) and their coupling
Close modal

Equation (9) represents a set of N tri-diagonal equations that can be solved readily using the tridiagonal matrix algorithm [23].

In order to derive similar equations for the cells in the substrate, the finite volume procedure is applied again, but now to the full two-dimensional equation, to yield
[kSΔzk(2rj+Δrj)Δrj+Δrj+1+kSΔzk(2rjΔrj)Δrj+Δrj1+4kSrjΔrjΔzk+Δzk+1+4kSrjΔrjΔzk+Δzk1+ρScSΔzk2rjΔrjΔt]TS,j,k[kSΔzk(2rj+Δrj)Δrj+Δrj+1]TS,j+1,k[kSΔzk(2rjΔrj)Δrj+Δrj1]TS,j1,k[4kSrjΔrjΔzk+Δzk+1]TS,j,k+1[4kSrjΔrjΔzk+Δzk1]TS,j,k1=ρScSΔzk2rjΔrjΔtTS,j,kold
(10)

where thermal properties with subscript “S” denote those of the substrate. Similar equations may be derived for cells adjacent to the boundaries. These equations are not presented here for the sake of brevity. Equation (10), along with similar equations for the boundary cells, represent a system of five-banded linear algebraic equations that may be solved using any iterative solver tailored for banded systems. In this particular case, the Stone's strongly implicit method [23] is used.

Equations (9) and (10) are coupled, and are solved using the following iterative procedure:

  1. The temperature on the top surface of the substrate, TS,j,top, is first guessed. A reasonable guess may be the temperature of the ambient, i.e., T.

  2. Equation (9) is then solved using a tridiagonal matrix algorithm. This yields the temperature of all cells of the transducer, namely TT,j .

  3. Equation (8) is next used to compute the heat flux through the interface or contact.

  4. The computed value of the flux serves as a boundary condition for the top surface of the substrate. With this boundary condition, Eq. (10), along with similar equations for the boundary cells in the substrate, are solved using the Stones's method.

  5. Steps 1–4 are repeated until convergence.

  6. Once convergence has been attained for the time step in question, this solution replaces the initial condition, and the procedure is repeated for the next time step.

2.3 Hyperbolic Heat Conduction Equation.

The Fourier law can be derived from the BTE in the limit where the number of scattering events of the heat carriers, i.e., phonons, is infinite. This limit may also be manifested by assuming an infinite group velocity of the phonons. The Fourier heat conduction equation does not contain any information pertaining to either the group velocity or the relaxation time scale of the phonons. It is incapable of capturing any physical phenomena where the finite speed of the phonons may be of importance.

To overcome the aforementioned limitations of infinite wave speed as predicted by the Fourier heat conduction equation, a modification to Eq. (5) was proposed by Catteneo [8,9], in which, an additional transient term was introduced
ρcτ2Tt2+ρcTt=k2T=k2Tz2+krr(rTr)
(11)

Equation (11) is a damped wave equation with wave speed k/ρcτ, and is referred to as the hyperbolic heat conduction equation. This equation can also be derived from the BTE by taking its first moment [8,9]. It assumes that the time scales of interest are of the same order of magnitude as the relaxation time, whereas the length scales are much larger than the characteristic length scale for local thermodynamic equilibrium. This makes the hyperbolic heat conduction equation nonlocal in time but local in space. Therefore, the ballistic behavior of phonons is only partially captured by this model. In addition, the frequency dependent behavior of phonons cannot be modeled using this equation.

As far as discretization and the numerical solution of Eq. (11) is concerned, the procedure is very similar to the one described in Sec. 2.2. The right-hand side is discretized and spatial boundary conditions are applied in exactly the same manner as discussed earlier. The second derivative in time is discretized using the backward Euler method, as described elsewhere [23]. The numerical solution also follows the same algorithm outlined for the Fourier heat conduction equation. Once the temperature distribution has been computed, the heat fluxes are computed using
τqt+q=kT
(12)

rather than the Fourier law, i.e., Eq. (12) is used instead of q=kT. One of the critical additional inputs to the hyperbolic heat conduction equation (Eq. (11)) is the effective relaxation time-scale of the phonons, τ. As to how it is determined is discussed in Sec. 3.

3 Results and Discussion

For the purposes of this study, the experimental data reported by Regner et al. [2] have been used for extraction of the thermal conductivity. The substrate in this experiment is a silicon block that is 525 μm thick, i.e., zS= 525 μm. The radial extent of substrate, rS, is not known, and was assumed to be also equal to 525 μm for the numerical calculations. The transducer is a bilayer transducer with 55 nm of gold and 5 nm of chromium, resulting in zT = 60 nm. The thermophysical properties of silicon, gold, and chromium that were used for calculations are shown in Table 1. Based on the thicknesses of the gold and chromium layers, effective values of the properties of the transducer were estimated and used. These are also shown in Table 1.

Table 1

Thermophysical properties of the various materials used in the calculations

Transducer
SiliconGoldChromiumEffective
Density (kg/m3)232919,320714018,290.8
Specific heat capacity (J/kg/K)702129450155.7
Thermal conductivity (W/m/K)Calculated31093.9266.6
Transducer
SiliconGoldChromiumEffective
Density (kg/m3)232919,320714018,290.8
Specific heat capacity (J/kg/K)702129450155.7
Thermal conductivity (W/m/K)Calculated31093.9266.6

For numerical calculations, a 200 × 200 nonuniform mesh with a stretching factor not exceeding 1.2 was used, as shown schematically in Fig. 2. This mesh was found to yield grid independent solutions. The modulation frequency, which is as an input parameter, was varied between 200 kHz and 200 MHz as in the experiments, and calculations were conducted for 22 different frequencies in this range. Each sinusoidal cycle of the laser was split into 5000 time steps. This implies that for high frequencies, a very small time-step size (∼1 ps) was used to ascertain accurate temporal solutions. The nominal value of the interfacial (between the substrate and the transducer) contact conductance, GC, was taken to be 200 MW/m2/K, as suggested by Cahill [1], although for a different material pair. The pump and probe laser 1/e2 radii, denoted by rpump and rprobe, respectively, are inputs in the model. A Gaussian laser flux profile (in r) was used for all calculations. The laser power was adjusted to attain a temperature rise of approximately 5 deg at the center of the laser spot, as reported in the experimental description [2]. For each modulation frequency, the calculations were advanced in time for several tens of cycles until the system exhibited quasi-periodic behavior. Three cases were considered:

  • Case 1 (baseline case): rpump = 4.1 μm, rprobe = 2.8 μm, resulting in reff = 4.95 μm

  • Case 2: rpump = 3.5 μm, rprobe = 3.5 μm, resulting in reff = 4.95 μm

  • Case 3: rpump = 2.8 μm, rprobe = 2.8 μm, resulting in reff = 3.95 μm

The radii for the baseline case were selected based on data presented in the supplementary material document of Regner et al. [2]. Cases 1 and 2 constitute a parametric variation in which the pump and probe radii are altered while keeping the effective radius unchanged. Cases 1 and 3, on the other hand, constitute a parametric variation in which the pump radius is varied while keeping the probe radius unchanged.

Figure 3 shows a plot of the temperature–time history of the center of the transducer (origin in Fig. 1) for four different frequencies for the baseline case. At low frequency, the transducer receives the laser flux for a longer period of time and, consequently, heats up more. This time-domain signal is postprocessed to calculate the phase shift or lag between the pump laser and the computed temperature at the center of the transducer. The same calculation was repeated for several different thermal conductivity (of substrate) values. For each thermal conductivity value, a phase lag versus frequency plot was generated. The error norms—both L1 norm and L2 norm—between the computed phase lag and the experimentally measured phase lag were then computed. The extracted (or desirable) thermal conductivity was chosen to be the one that minimized the error norm. In other words, it is a value that best fits the data over the entire range of frequency considered in this study. It was found that depending on which norm (L1 or L2) is used for best fit, the extracted value of thermal conductivity is slightly different, as discussed shortly. To improve the quality of fit, previous researchers have suggested that, perhaps, one should consider a frequency-dependent thermal conductivity [5], while others have suggested using an anistropic thermal conductivity [24]. Yet others have tried to adjust the interface thermal conductance to obtain a better fit [11]. Wilson et al. [25], and later Regner et al. [26], proposed and used a more sophisticated nonequilibrium two-temperature model in which phonons and electrons have different temperature, and exchange of energy at the interface between the metal transducer and the semiconductor substrate is driven by the difference in these temperatures. While such a model better represents the energy transport at the interface, it requires solution of two transient heat equations, which is beyond the scope of this study, and also has additional unknown (uncertain) parameters.

Fig. 3
Computed (Fourier equation) temperature history at the center of the top surface of the transducer for different frequencies for the baseline case
Fig. 3
Computed (Fourier equation) temperature history at the center of the top surface of the transducer for different frequencies for the baseline case
Close modal

The sensitivity of the extracted thermal conductivity to the pump and probe radii was first investigated, since in many experiments, these quantities are often uncertain. Cases 1 and 3 constitute a parametric variation in which the pump radius is varied, while the probe radius is unchanged. Furthermore, in both cases, the output temperature signals were processed in two ways: the center spot of the transducer (corresponds to the limiting case of rprobe=0), and an average over rprobe=2.8μm. Figure 4(a) shows the phase lag computed for the two different pump radii and rprobe=2.8μm for k = 98 W/m/K. As seen, the computed phase lag is sensitive to the pump laser spot size, especially at low and intermediate frequencies when the thermal wave gets an opportunity to penetrate more, and the radial transport of heat becomes more pronounced with a smaller laser spot. Figure 4(b) shows the phase lag computed with two different probe/pump radii combinations, while maintaining the same effective radius of 3.5 μm. Again, k = 98 W/m/K is used to generate this phase plot. It is clear from Fig. 4(b) that the numerically calculated phase is affected by individual values of the pump and probe radii, as opposed to just the effective radius in the analytical model. The extracted thermal conductivity values that result in best fit to the experimental data are shown in Table 2. For the baseline case, the analytically and numerically extracted thermal conductivity values are in excellent agreement: 99 W/m/K versus 98 W/m/K. The fact that these two results start deviating from each other in the other cases is simply due to the fact that the experimental data are for the baseline case and using it for extracting the thermal conductivity under other conditions is not valid, strictly speaking. The analytical value also agrees with the reported value [2] of k = 99±6 W/m/K, and the phase plots for these cases are shown in Fig. 5. The numerical model curve represents the best fit to the experimental data using the L2 norm minimization criterion. For the numerical model, the best fit is obtained by a thermal conductivity value of k = 98 W/m/K. If instead, the L1 norm minimization criterion is used, the best fit is obtained with a thermal conductivity value of 99 W/m/K. From Fig. 4 and Table 2, it is also clear that both pump and probe radii have noticeable effects on the predicted phase lag and the resulting thermal conductivity. Most importantly, the numerical results show that effective radius alone is not adequate to extract the thermal conductivity and variation of individual pump and probe radii do have an impact on the extracted thermal conductivity. Any uncertainty in the probe and pump radii can have a significant impact on the results.

Fig. 4
Numerically computed phase lags for the different cases considered: (a) Cases 1 and 3 showing effect of variation of pump radius and (b) Cases 1 and 2 showing effect of variation of both pump and probe radii but same effective radius. In each case, k = 98 W/m/K was used.
Fig. 4
Numerically computed phase lags for the different cases considered: (a) Cases 1 and 3 showing effect of variation of pump radius and (b) Cases 1 and 2 showing effect of variation of both pump and probe radii but same effective radius. In each case, k = 98 W/m/K was used.
Close modal
Fig. 5
Phase lag calculated using the analytical and numerical (Fourier) models: the numerical solution is for k = 98 W/m/K, while the analytical solutions are with k = 99 ± 6 W/m/K, with the small dashed lines showing the ±6 W/m/K envelop
Fig. 5
Phase lag calculated using the analytical and numerical (Fourier) models: the numerical solution is for k = 98 W/m/K, while the analytical solutions are with k = 99 ± 6 W/m/K, with the small dashed lines showing the ±6 W/m/K envelop
Close modal
Table 2

Effect of pump and probe laser radius on the extracted thermal conductivity (in W/m/K)

rpump (μm)
2.83.54.1
rprobe (μm)AnalyticalNumericalAnalyticalNumericalAnalyticalNumerical
0476465707873
2.87682889998
3.599108109
rpump (μm)
2.83.54.1
rprobe (μm)AnalyticalNumericalAnalyticalNumericalAnalyticalNumerical
0476465707873
2.87682889998
3.599108109

The extracted value of thermal conductivity, namely, k = 98 W/m/K, is significantly smaller than the bulk thermal conductivity of silicon at 300 K, which is about 148 W/m/K. This reduction in thermal conductivity represents the so-called thermal conductivity suppression reported in the literature [6,7]. The phase lag versus frequency behavior with the bulk value of thermal conductivity is shown in Fig. 6. The phase lag is underpredicted at all frequencies if the bulk value of 148 W/m/K for thermal conductivity is used.

Fig. 6
Phase lag computed using the numerical model (Fourier equation) with two different thermal conductivity values
Fig. 6
Phase lag computed using the numerical model (Fourier equation) with two different thermal conductivity values
Close modal

The spatiotemporal evolution of the nondimensional temperature distribution and the heat wave is illustrated in Fig. 7 for a modulation frequency of 200 kHz. For clarity of visualization, only a small portion of the substrate is shown. Since the penetration depth of the heat wave is inversely proportional to the inverse square root of the modulation frequency, and 200 kHz is the smallest frequency considered in this study, the distributions shown in Fig. 7 represent the highest penetration of the heat wave among all cases considered. After 30 cycles, the system has almost reached quasi-steady-state. Since the substrate (and computational domain) is of size 525 μm in both directions, and the heat wave has not even penetrated 200 μm in the worst-case scenario, it is fair to conclude that the computational domain is large enough for the far-field boundary conditions at the top and side to have minimal effect on the results. This is corroborated using additional calculations described in the next paragraph. The computed 1/e penetration depths were compared against those reported by Braun and Hopkins [21], although the present system includes a transducer layer, while the calculations reported in [21] does not. At pump laser frequencies of 200 kHz and 10 MHz, the penetration depths reported in [21] are 7 μm and 1.8 μm, respectively. In contrast, because of the presence of a transducer (and the contact resistance between the transducer and the substrate), the penetration depths in the present study was found to be expectedly lower: 2.4 μm and 1.1 μm, respectively, for the same two frequencies.

Fig. 7
Nondimensional temperature distributions computed using the Fourier heat conduction equation in the substrate after various instances of time with a modulation frequency of 200 kHz for the baseline case. Dimensions shown are in micrometer.
Fig. 7
Nondimensional temperature distributions computed using the Fourier heat conduction equation in the substrate after various instances of time with a modulation frequency of 200 kHz for the baseline case. Dimensions shown are in micrometer.
Close modal

For baseline numerical calculations with rpump=4.1μm and averaging over rprobe=2.8μm, the top surface of the transducer outside of the pump laser spot was assumed to have heat loss with ht = 10 W/m2/K. On the other hand, the far-field boundaries on the side and bottom were assumed to be isothermal at the ambient temperature (T= 300 K), which essentially corresponds to a very high heat transfer coefficient, i.e., hs,hb. To investigate the sensitivity of this boundary condition on the computed results, the solution was recomputed with an adiabatic boundary condition instead of an isothermal boundary condition. An adiabatic boundary condition represents the other extreme wherein the heat transfer coefficients are zero. Figure 8 shows the difference in the temperature along the centerline (axis) of the substrate. Clearly, there is no meaningful difference implying that the boundary condition has no impact on the results. This may be attributed to the fact that the computational domain is much larger than the penetration depth. Likewise, the heat transfer coefficient on the top of the transducer was changed from 10 W/m2/K to 250 W/m2/K, and the maximum change in temperature anywhere on the transducer was found to be 5.5 μK. The general trend of these results appears to agree with the findings of Xing et al. [20].

Fig. 8
Computed temperature difference (Fourier equation) along the centerline of the substrate between adiabatic and isothermal boundary condition at the bottom surface at quasi-steady-state
Fig. 8
Computed temperature difference (Fourier equation) along the centerline of the substrate between adiabatic and isothermal boundary condition at the bottom surface at quasi-steady-state
Close modal

One of the critical inputs to the model considered here is the interface conductance, GC. For the baseline calculations, a value of 200 MW/m2/K was used, in accordance with Cahill [1]. Other researchers have suggested a value in the range 160–250 MW/m2/K [11]. To understand the sensitivity of the extracted thermal conductivity to this unknown parameter, computations were performed for various values of GC. Figure 9 shows the predicted phase lag for various values of GC and a thermal conductivity of 98 W/m/K. There is clearly some difference in the results at intermediate and high modulation frequencies. The quality of the fit is best for GC= 200 MW/m2/K than for the other two values.

Fig. 9
Effect of interface conductance on phase lag computed using the Fourier heat conduction equation for the baseline case
Fig. 9
Effect of interface conductance on phase lag computed using the Fourier heat conduction equation for the baseline case
Close modal
As a final step in the analysis, all calculations were repeated using the hyperbolic heat conduction equation. The effective (or average) relaxation time-scale, τ, which appears as an input in the hyperbolic heat conduction equation, was estimated using the following relationship [27]
τ=pωmin,pωmax,pcω,p|υω,p|2τω,pdωpωmin,pωmax,pcω,p|υω,p|2dω
(13)

where cω,p is the spectral specific heat capacity of silicon, and υω,p is the spectral group velocity of the phonons. Both of these quantities can be computed directly for any polarization p and frequency ω using dispersion relationships of silicon, as explained elsewhere [8,28]. The value of the effective relaxation time-scale using Eq. (13) and frequency and polarization dependent time-scale expressions proposed by Holland [29,30] was found to be about 0.0101 ns. The error (difference) in temperature as predicted by the two models (Fourier versus hyperbolic) is depicted in Fig. 10. The error in temperature shown in the y-axis is the difference normalized by the peak-to-peak temperature of each case, and expressed as a percentage. The nondimensional distance shown in the x-axis is the distance normalized by the theoretical penetration depth for the lowest frequency (200 kHz) case. At high frequency, quasi-ballistic effects are expected to be the strongest and consequently, some difference is observed between the two models. At low frequency, the errors are practically nonexistent. At 200 MHz, the cycle time is 5 ns, which is still significantly larger than the estimated effective relaxation time scale. This is the reason why even for the highest frequency case, the largest error is only about 0.12%. The phase lag calculated using the two models was found to be almost identical, and the resulting thermal conductivity computed using the hyperbolic heat conduction model was found to be 97 W/m/K (as opposed to 98 W/m/K with the Fourier heat conduction model). Again, this closeness in the extracted value of thermal conductivity is to be expected since the relaxation time-scale used in the calculations is significantly smaller than the cycle time for any modulation frequency considered.

Fig. 10
Computed temperature difference (percentage) along the centerline of the substrate between the Fourier heat conduction and hyperbolic heat conduction models for various modulation frequencies after 30 cycles for the baseline case
Fig. 10
Computed temperature difference (percentage) along the centerline of the substrate between the Fourier heat conduction and hyperbolic heat conduction models for various modulation frequencies after 30 cycles for the baseline case
Close modal

The choice of the effective relaxation time-scale was further investigated. Figure 11 shows a plot of the relaxation time scales of the two acoustic phonon types along with the cumulative energy distribution function computed at 300 K. It shows that a significant amount of energy is carried by low-frequency phonons, which have a much larger relaxation time-scale compared to the average value. Based on this observation, computations were performed with τ = 0.05 ns instead of 0.01 ns. With this relaxation time scale, the thermal conductivity that exhibited best fit to the phase plot was still found to be 97 W/m/K, although the quality of the fit was not as good as the lower average time scale case. These results indicate that the conditions under which this particular set of experiments were conducted, are such that ballistic effects, if considered on an average, are almost negligible. Perhaps, phonons of extremely low frequency exhibit this effect, and this can only be captured by a higher fidelity model such as the BTE.

Fig. 11
Relaxation time scales of acoustic phonons computed using [26,27] and the cumulative energy distribution function computed at 300 K for silicon
Fig. 11
Relaxation time scales of acoustic phonons computed using [26,27] and the cumulative energy distribution function computed at 300 K for silicon
Close modal

4 Summary and Conclusions

Calculations were conducted both analytically and numerically to extract the thermal conductivity of a silicon substrate from FDTR experimental data. The analytical calculations used Hankel transforms of the Fourier heat conduction equation in frequency domain along with the Feldman algorithm to treat multiple finite-sized layers. The numerical calculations were conducted in time domain using the finite volume method and a backward Euler time advancement scheme. Both the Fourier heat conduction equation and the hyperbolic heat conduction equation were explored for numerical solutions. The numerical solutions enable use of more realistic boundary conditions, namely, convective heat loss at the surfaces of the substrate and transducer, and can also isolate the effect of variation of the spot sizes of the pump and probe lasers. The calculations were conducted in the range of modulation frequency going from 200 kHz to 200 MHz for 22 different frequencies. The phase lag between the modulated pump laser signal and the temperature of the transducer averaged over the probed spot was calculated in each case and a plot of phase lag versus frequency was generated. This plot was compared to the same plot measured experimentally, and the thermal conductivity was adjusted to obtain the best fit. The numerically extracted thermal conductivity was found to be in the range 82–108 W/m/K, depending on the pump and probe laser spot sizes used, choice of boundary conditions, and the interface conductance between the transducer and the substrate used in the calculations. For the baseline case, the extracted value agrees almost perfectly with the value of 99±6 W/m/K reported earlier [2] for the same experimental dataset. Parametric studies revealed that the extracted thermal conductivity is sensitive to both the pump and probe laser spot sizes, and the effect cannot be captured by an effective radius alone. It is quite insensitive to the boundary conditions at the surfaces of the substrate/transducer. For this particular case, the hyperbolic heat conduction equation yielded almost the same results as the Fourier heat conduction equation, with errors in predicted temperature below 0.12%.

Acknowledgment

The Ohio Supercomputer Center (OSC) and the Department of Mechanical and Aerospace Engineering at the Ohio State University are gratefully acknowledged for providing computational resources.

Funding Data

  • National Science Foundation Directorate for Computer and Information Science and Engineering (Award No. 2003747) (Funder ID: 10.13039/100000083).

Nomenclature

c =

specific heat capacity (J/kg/K)

h =

heat transfer coefficient (W/m/K)

k =

thermal conductivity (W/m/K)

n =

time index

p =

Hankel transform parameter

r =

radial coordinate (m)

z =

axial coordinate (m)

q =

heat flux vector (W/m2/K)

A =

amplitude factor (see Eq. (1))

N =

number of layers

T =

temperature (K)

hb =

external heat transfer coefficient at the bottom surface of the tank (W/m2/K)

hs =

external heat transfer coefficient at the side surface of the tank (W/m2/K)

ht =

external heat transfer coefficient at the top surface of the tank (W/m2/K)

reff =

effective laser radius (m)

rprobe =

probe laser 1/e2 radius (m)

rpump =

pump laser 1/e2 radius (m)

GC =

thermal contact conductance (W/m/K)

Ln =

thickness of nth layer (m)

T =

ambient temperature (K)

zT =

thickness of transducer (m)

qL =

laser heat flux (W/m2/K)

Δt =

time-step size (sec)

Greek Symbols

Greek Symbols
α =

thermal diffusivity (m2/s)

ρ =

density (kg/m3)

τ =

relaxation time scale (s)

ω =

modulation frequency (Hz)

References

1.
Cahill
,
D. G.
,
2004
, “
Analysis of Heat Flow in Layered Structures for Time Domain Thermo-Reflectance
,”
Rev. Sci. Instrum.
,
75
, p.
5119
.10.1063/1.1819431
2.
Regner
,
K. T.
,
Sellan
,
D. P.
,
Su
,
Z.
,
Amon
,
C. H.
,
McGaughey
,
A. J. H.
, and
Malen
,
J. A.
,
2013
, “
Broadband Phonon Mean Free Path Contributions to Thermal Conductivity to Thermal Conductivity Measured Using Frequency Domain Thermoreflectance
,”
Nat. Commun.
,
4
, p.
1640
.10.1038/ncomms2630
3.
Minnich
,
A. J.
,
Johnson
,
J. A.
,
Schmidt
,
A. J.
,
Esfarjani
,
K.
,
Dresselhaus
,
M. S.
,
Nelson
,
K. A.
, and
Chen
,
G.
,
2011
, “
Thermal Conductivity Spectroscopy Technique to Measure Phonon Mean Free Paths
,”
Phys. Rev. Lett.
,
107
, p.
095901
.10.1103/PhysRevLett.107.095901
4.
Feldman
,
A.
,
1999
, “
Algorithm for Solutions of the Thermal Diffusion Equation in a Stratified Medium With a Modulated Heating Source
,”
High Temp. High Press.
,
31
, pp.
293
298
.10.1068/htrt171
5.
Koh
,
Y. K.
, and
Cahill
,
D. G.
,
2007
, “
Frequency Dependence of the Thermal Conductivity of Semiconductor Alloys
,”
Phys. Rev. B
,
76
, p.
075207
.10.1103/PhysRevB.76.075207
6.
Minnich
,
A. J.
,
2012
, “
Determining Phonon Mean Free Paths From Observations of Quasiballistic Thermal Transport
,”
Phys. Rev. Lett.
,
109
, p.
205901
.10.1103/PhysRevLett.109.205901
7.
Hua
,
C.
, and
Minnich
,
A. J.
,
2014
, “
Transport Regimes in Quasiballistic Heat Conduction
,”
Phys. Rev. B
,
89
, p.
094302
.10.1103/PhysRevB.89.094302
8.
Tien
,
C. L.
,
Majumdar
,
A.
, and
Gerner
,
F. M.
, eds.,
1998
,
Microscale Energy Transport
,
Taylor and Francis
,
Washington, DC
.
9.
Zhang
,
Z.
,
2007
,
Nano/Microscale Heat Transfer
,
McGraw-Hill
,
New York
.
10.
Ramu
,
A. T.
, and
Bowers
,
J. E.
,
2015
, “
A Compact Heat Transfer Model Based on an Enhanced Fourier Law for Analysis of Frequency Domain Thermoreflectance Experiments
,”
Appl. Phys. Lett.
,
106
, p.
263102
.10.1063/1.4923310
11.
Ma
,
Y.
,
2014
, “
A Two-Parameter Nondiffusive Heat Conduction Model for Data Analysis in Pump–Probe Experiments
,”
J. Appl. Phys.
,
116
, p.
243505
.10.1063/1.4904355
12.
Chen
,
G.
,
2001
, “
Ballistic-Diffusive Heat-Conduction Equations
,”
Phys. Rev. Lett.
,
86
(
11
), pp.
2297
2300
.10.1103/PhysRevLett.86.2297
13.
Mittal
,
A.
, and
Mazumder
,
S.
,
2011
, “
Generalized Ballistic-Diffusive Formulation and Hybrid SN-PN Solution of the Boltzmann Transport Equation for Phonons for Non-Equilibrium Heat Conduction
,”
ASME J Heat Transfer-Trans. ASME
,
133
(
9
), p.
092402
.10.1115/1.4003961
14.
Mittal
,
A.
, and
Mazumder
,
S.
,
2011
, “
Hybrid Discrete Ordinates—Spherical Harmonics Solution to the Boltzmann Transport Equation for Phonons for Non-Equilibrium Heat Conduction
,”
J. Comput. Phys.
,
230
(
18
), pp.
6977
7001
.10.1016/j.jcp.2011.05.024
15.
Hua
,
C.
,
Lindsay
,
L.
,
Chen
,
X.
, and
Minnich
,
A. J.
,
2019
, “
Generalized Fourier's Law for Nondiffusive Thermal Transport: Theory and Experiment
,”
Phys. Rev. B
,
100
, p.
085203
.10.1103/PhysRevB.100.085203
16.
Romano
,
G.
,
Kolpak
,
A. M.
,
Carrete
,
J.
, and
Broido
,
D.
,
2019
, “
Parameter-Free Model to Estimate Thermal Conductivity in Nanostructured Materials
,”
Phys. Rev. B
,
100
, p.
045310
.10.1103/PhysRevB.100.045310
17.
Torres
,
P.
,
Ziabari
,
A.
,
Torelló
,
A.
,
Bafaluy
,
J.
,
Camacho
,
J.
,
Cartoixà
,
X.
,
Shakouri
,
A.
, and
Alvarez
,
F. X.
,
2018
, “
Emergence of Hydrodynamic Heat Transport in Semiconductors at the Nanoscale
,”
Phys. Rev. Mater.
,
2
, p.
076001
.10.1103/PhysRevMaterials.2.076001
18.
Beardo
,
A.
,
Hennessy
,
M. G.
,
Sendra
,
L.
,
Camacho
,
J.
,
Myers
,
T. G.
,
Bafaluy
,
J.
, and
Alvarez
,
F. X.
,
2020
, “
Phonon Hydrodynamics in Frequency-Domain Thermoreflectance Experiments
,”
Phys. Rev. B
,
101
, p.
075303
.10.1103/PhysRevB.101.075303
19.
Ali
,
S. A.
, and
Mazumder
,
S.
,
2017
, “
Phonon Boltzmann Transport Equation Based Modeling of Time Domain Thermo-Reflectance Experiments
,”
Int. J. Heat Mass Transfer
,
107
, pp.
607
621
.10.1016/j.ijheatmasstransfer.2016.11.077
20.
Xing
,
C.
,
Jensen
,
C.
,
Hua
,
Z.
,
Ban
,
H.
,
Hurley
,
D. H.
,
Khafizov
,
M.
, and
Kennedy
,
J. R.
,
2012
, “
Parametric Study of the Frequency-Domain Thermoreflectance Technique
,”
J. Appl. Phys.
,
112
, p.
103105
.10.1063/1.4761977
21.
Braun
,
J. L.
, and
Hopkins
,
P. E.
,
2017
, “
Upper Limit to the Thermal Penetration Depth During Modulated Heating of Multilayer Thin Films With Pulsed and Continuous Wave Lasers: A Numerical Study
,”
J. Appl. Phys.
,
121
, p.
175107
.10.1063/1.4982915
22.
Carslaw
,
H. S.
, and
Jaeger
,
J. C.
,
1959
,
Conduction of Heat in Solids
, 2nd ed.,
Oxford University Press
,
London
, UK.
23.
Mazumder
,
S.
,
2016
,
Numerical Methods for Partial Differential Equations: Finite Difference and Finite Volume Methods
, 1st ed.,
Academic Press
,
New York
.
24.
Feser
,
J. P.
, and
Cahill
,
D. G.
,
2012
, “
Probing Anisotropic Heat Transport Using Time-Domain Thermoreflectance
,”
Rev. Sci. Instrum.
,
83
, p.
104901
.10.1063/1.4757863
25.
Wilson
,
R. B.
,
Fesser
,
J. P.
,
Hohensee
,
G. T.
, and
Cahill
,
D. G.
,
2013
, “
Two-Channel Model for Nonequilibrium Thermal Transport in Pump-Probe Experiments
,”
Phys. Rev. B
,
88
, p.
144305
.10.1103/PhysRevB.88.144305
26.
Regner
,
K. T.
,
Wei
,
L. C.
, and
Malen
,
J. A.
,
2015
, “
Interpretation of Thermoreflectance Measurements With a Two-Temperature Model Including Non-Surface Heat Deposition
,”
J. Appl. Phys.
,
118
, p.
235101
.10.1063/1.4937995
27.
Allu
,
P.
, and
Mazumder
,
S.
,
2018
, “
Comparative Assessment of Deterministic Approaches to Modeling Quasiballistic Phonon Heat Conduction in Multi-Dimensional Geometry
,”
Int. J. Therm. Sci.
,
127
, pp.
181
193
.10.1016/j.ijthermalsci.2018.01.024
28.
Mazumder
,
S.
, and
Majumdar
,
A.
,
2001
, “
Monte Carlo Study of Phonon Transport in Solid Thin Films Including Dispersion and Polarization
,”
ASME J Heat Transfer-Trans. ASME
,
123
, pp.
749
759
.10.1115/1.1377018
29.
Holland
,
M. G.
,
1963
, “
Analysis of Lattice Thermal Conductivity
,”
Phys. Rev.
,
132
(
6
), pp.
2461
2471
.10.1103/PhysRev.132.2461
30.
Holland
,
M. G.
,
1964
, “
Phonon Scattering in Semiconductors From Thermal Conductivity Studies
,”
Phys. Rev.
,
134
(
2A
), pp.
A471
A480
.10.1103/PhysRev.134.A471