On the Determination of Thermal Conductivity From Frequency Domain Thermoreflectance Experiments

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 [1–3]. 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 [1–3,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 [5–7].

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 [8–18]. 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.

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.

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.

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

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

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.

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.

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.

rather than the Fourier law, i.e., Eq. (12) is used instead of $q=−k∇T$⁠. 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.

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.

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/m²/K, as suggested by Cahill [1], although for a different material pair. The pump and probe laser 1/e² 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.

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.

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.

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.

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/m²/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/m²/K to 250 W/m²/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].

One of the critical inputs to the model considered here is the interface conductance, $GC$⁠. For the baseline calculations, a value of 200 MW/m²/K was used, in accordance with Cahill [1]. Other researchers have suggested a value in the range 160–250 MW/m²/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/m²/K than for the other two values.

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.

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.

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%.

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.

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

$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/m²/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/m²/K)

$hs$ =

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

$ht$ =

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

$reff$ =

effective laser radius (m)

$rprobe$ =

probe laser 1/e² radius (m)

$rpump$ =

pump laser 1/e² radius (m)

$GC$ =

thermal contact conductance (W/m/K)

$Ln$ =

thickness of nth layer (m)

$T∞$ =

ambient temperature (K)

$zT$ =

thickness of transducer (m)

$q″L$ =

laser heat flux (W/m²/K)

$Δt$ =

time-step size (sec)

Greek Symbols

$α$ =

thermal diffusivity (m²/s)

$ρ$ =

density (kg/m³)

$τ$ =

relaxation time scale (s)

$ω$ =

modulation frequency (Hz)