Upper limit to the thermal penetration depth during modulated heating of multilayer thin films with pulsed and continuous wave lasers: A numerical study

Pulsed and continuous wave (CW) laser applications are ubiquitous in both industry and scientific research. Many modern applications rely on laser heating to manipulate structural, mechanical, and electrical properties of materials. Pulsed laser heating, having versatile windows of energy duration ranging from femtoseconds to nanoseconds, is routinely used in materials and microstructure fabrication,^1–3 materials processing,⁴ laser cutting and micro-machining,⁵ and property measurements through pump-probe techniques.^6–9 Additionally, with picosecond and sub-picosecond pulse durations, pulsed lasers can be used to provide high energy densities locally, providing a means for surface modification.^10,11 Continuous wave laser heating, likewise, has major applications in welding, cutting, surface treatment¹² such as near-surface melting to modify wear and corrosion,¹³ and surface patterning.² In all such applications, it is necessary to understand the spatially varying temperature profile that results from laser heating, since this temperature dictates the material response.

In this study, we derive a method for predicting this spatially dependent temperature profile based on the framework of the radially symmetric heat diffusion equation. From this methodology, we examine spatially dependent temperature profiles to study the thermal penetration depth (TPD) in bulk and thin film systems. Our results provide answers to the questions of how the TPD is dependent on time and length scales of the laser heat source and sample geometry. In the case of optical pump-probe techniques that rely on locking into a frequency modulated surface heating event, namely, time- and frequency-domain thermoreflectance (TDTR and FDTR, respectively),^7,14–16 the TPD is of great importance as it relates a measurement length scale to the length scale of physical phenomena of interest. For example, interpretation of TDTR/FDTR extracted thermal properties' relation to physical processes relies on the idea that the TPD resulting from a modulated heating event can be orders of magnitude smaller than the pump radius, making it the defining length scale of the thermal response within the measurement volume. In fact, many of the recent advances in mean free path spectroscopy¹⁷ stem from relating phonon mean free paths to the TPD imposed by a modulated heating event.^18–23 Therefore, it is of utmost importance to properly characterize this length scale. As we will show in the discussion of our results, this length scale is highly dependent on a multitude of variables that make accurate analytical predictions difficult for realistic sample geometries and conditions, suggesting a need for the rigorous calculations presented in this study. Such variables include time, modulation frequency, material and interface properties, and heater radius.

The 1/e temperature decay TPD used to interpret TDTR and FDTR results is often estimated using the analytical solution to a one-dimensional semi-infinite sample subjected to a CW heating source modulated in the high frequency limit (⁠$ω0≫D/r02$⁠, where ω₀ is modulation angular frequency, r₀ is pump radius, and D is thermal diffusivity) taking the form

However, this equation breaks down when experimental conditions deviate from those on which this expression was derived. Considering commonly studied materials like silicon (Si), amorphous silicon dioxide (a-SiO₂), and aluminum (Al) on Si, a-SiO₂ reveals that such deviations can occur at modulation frequencies and pump radii common in TDTR and FDTR. To comply with Eq. (1) and the definition typically used in TDTR/FDTR, we refer to the TPD as the distance away from the spatial coordinate of maximum temperature at which the temperature profile has decayed to 1/e of this maximum temperature; we observe the effects of altering this definition as well.

In order to establish the upper limit to the obtainable TPD from a modulated heating event, we first discuss the TPD associated with an unmodulated heating event. Note that since we are considering both modulated and unmodulated temperature profiles in this study, we generalize the definition of TPD so that we can quantify the temperature decay both radially and in depth. We do this because, when considering the modulated temperature profile simulating the heating event to which lock-in detection occurs during TDTR and FDTR, the temperature profile is typically one-dimensional; that is, the temperature decay extends in depth much less than radially. However, as the modulation frequency approaches zero, the TPD becomes large, on the order of pump radius, making the radial decay length important. Before discussing these concepts in more detail, we first establish the framework on which all calculations are based.

To begin, the heat diffusion equation in cylindrical coordinates with radial symmetry and a surface heating source is given by

where κ_r and κ_z are in-plane and cross-plane thermal conductivity, respectively, T is the temperature (relative to some initial temperature) at a given location and time, r denotes radius (in-plane), z denotes depth (cross-plane) and is orthogonal to r, C_v is the volumetric heat capacity, and t is the time. For a surface adjacent to an insulating medium (such as air), the boundary conditions are given by

where T_top(r, t) and Q_top(r, t) are the temperature and heat flux at the top surface (z = 0) of a given material. The initial condition is given by T(z, r, 0) = 0, where 0 is defined as the reference temperature to which Eq. (2) describes changes in temperature thereof. The heating source term, Q_top(r, t), is defined by

such that Q_top is taken to be the absorbed heat flux supplied by the laser, Q_laser. This term is decoupled into a radial distribution and a time varying function, G(t). The right hand side of Eq. (5) assumes that the radial distribution takes on a Gaussian profile with r₀ describing the 1/e² radius of the pump beam. As the boundary conditions suggest, in this study, we focus on the situation in which the laser energy is absorbed near the surface, which is the typical condition used in TDTR and FDTR where a metal thin-film transducer absorbs the laser energy. In the situation when this is not the case, modifications to the energy absorption profile must be implemented, as described elsewhere.^24,25

We will consider two cases for G(t), a pulsed source and a CW source. In the pulsed case, the time dependence of the beam is assumed Gaussian such that G(t) can be mathematically represented as

where the absorbed laser energy at the surface is given by the absorption of the sample, α, multiplied by energy per pulse, E, defined as A/f_s, where A is the time-averaged laser power, and f_s is the laser repetition rate. Finally, σ is the pulse width standard deviation and relates to the experimentally obtainable full-width at half maximum (FWHM) time of the pulse width, through $FWHM=22ln(2)σ$⁠. In this study, we consider only sub-picosecond pulses, such that the approximation made on the right-hand side of Eq. (6) can be assumed without significant loss of accuracy. The sum is taken over all pulses to account for any accumulation of heat over time. The continuous wave laser heating event is simple by comparison, given by

such that G(t) is constant with time. In both cases, G(t) is described in units of [W] such that Q_top(r, t) obtains the desired units of [W/m²].

Typically in optical pump-probe experiments such as TDTR and FDTR, a modulation at angular frequency ω₀ [rad/s] is applied to the pump beam. To keep this analysis broad, we derive all equations using a modulated heat source, noting that the solution is equivalent to that of an unmodulated source when ω₀ = 0. Thus, given the same absorbed average power of αA, the fully extinguished sinusoidally modulated source terms for pulsed and CW sources, respectively, can be expressed as

These power functions are displayed in the inset to Fig. 1(a). Given the radial symmetry and the periodic nature of the time varying heat source, Eq. (2) is conveniently solved in Hankel and frequency space. Taking the Hankel and Fourier transforms, the heat diffusion equation and boundary condition, after simplification, become

where k denotes the Hankel transform variable and ω denotes the frequency space variable. The notation used here, consistent with a previous work,²⁶ is such that for a function $Y(z,r,t), Ỹ≡Ỹ(z,r,ω)$ denotes the Fourier transform of Y, while $Ỹ̃≡Ỹ̃(z,k,ω)$ denotes the Hankel transform of $Ỹ$⁠. For notational clarity, we drop the k and ω dependence of transformed functions. The solution to this ordinary differential equation in z, after applying boundary conditions and noting that the heat flux component in depth is $Q̃̃(z)=−κz∂T̃̃∂z$⁠, is given in matrix form as

This solution is generalized to layered systems by noting that Eq. (13) is valid for any layer in which $T̃̃top$ and $Q̃̃top$ for that layer are known. While the heat flux is prescribed at layer 1's surface, the temperature at this surface is unknown. To obtain this quantity, we invoke a semi-infinite boundary condition at the back side of the final layer of some system having n layers (refer to Fig. 1(a)), such that the following condition holds:

Applying this condition to the nth layer allows the determination of $T̃̃top,n$ and $Q̃̃top,n$⁠, defined to be the temperature and heat flux component in depth at the surface of layer n, i.e., at $z=∑j=1n–1dj$⁠. These quantities can then be related to the bottom temperature and heat flux of the (n – 1)th layer through the thermal boundary conductance (TBC), h_(n–1)/n, mathematically represented in matrix notation by

so that the bottom temperature and heat flux of the (n – 1)th layer are known, allowing for the solution to the temperature and heat flux at the top of this layer as well. This process can be repeated for an arbitrary number of layers until ultimately the temperature at the top of layer 1 is determined. This matrix approach^14,27–29 can then be generalized so that the temperature and heat flux within the nth layer can be determined from

Here, z_n is defined as $zn=z−∑jn–1dj$ such that z_n = 0 corresponds to the top surface of the nth layer. Note that this set of equations can be used to solve the temperature profile in any given layer once T_top is obtained, since n can represent any layer of interest. To proceed, though, without loss of generality we assume the nth layer denotes the final, semi-infinite layer of this geometry. Then, $T̃̃top$ is determined by invoking the aforementioned semi-infinite boundary condition described in Eq. (14), such that

Since $Q̃̃top$ is a well-defined function given by Eq. (12), $T̃̃top$ can be obtained. To obtain this solution in radial space and time, we take the inverse Hankel and Fourier transforms on $T̃̃top$⁠. Taking first the inverse Hankel transform yields

where J₀ is the zeroth-oder Bessel function of the first kind and $L̃(r,ω)$ is defined to simplify notation. This Hankel transform can be solved numerically; the upper bound of the integral can be set to 2/r₀ without significant loss of accuracy.¹⁴ Next, taking the inverse Fourier transform yields the desired temperature solution in radial and time coordinates

Thus, T_top is completely characterized, in principle. In practice, there are several simplifying assumptions that can be made based on the specific form of $G̃(ω)$⁠.

Beginning with the modulated CW source described by Eq. (9), we write $G̃(ω)$ as

whereas the modulated pulsed source given by Eq. (8) becomes

where ω_s is defined as 2πf_s and $F$ denotes the Fourier transform operation. Interestingly, in the frequency domain the only difference between the pulsed and CW source terms is the additional sum over n pulses for the pulsed case. The delta function representation of these expressions makes the transformation from frequency space to time relatively simple. In the CW case, the expression for T_top becomes

whereas in the pulsed case, the expression for T_top takes the form

It is clear that in both the CW and pulsed modulated heating cases, there are two contributions to the temperature rise at the surface: a steady state, “DC” temperature rise governed by the average power deposited into the sample, and a modulated, “AC” temperature rise associated with frequency modulation. Decoupling the terms and recognizing that T_top = T_top,DC + T_top,AC, the temperature rises for the CW case are given by

while those of the pulsed case are given by

These equations completely describe the temperature rise at the top surface of layer 1. With T_top known and Q_top prescribed, T(z, r, t) can be explicitly calculated in any layer within a sample having layered structures by invoking Eq. (16) and transforming back to real space and time.

Before determining the spatially dependent temperature profile, we first examine the modulated surface temperature profile associated with TDTR and FDTR to show the effects of pulse accumulation and modulation frequency on this temperature response. We note that Schmidt et al.²⁹ examined this with a simple exponential thermal response for demonstrative purposes; here, we show the realistic temperature response of a 100 nm Al film on a Si substrate, a typical system studied in TDTR and FDTR. Although the temperature distribution contains modulated and unmodulated components, to correspond with the lock-in amplification process used in TDTR/FDTR, we consider only the modulated component. Additionally, the photodetector measures a reflectance signal proportional to a gaussian-averaged temperature distribution about the probe radius, r₁, represented by

This modulated temperature rise, along with the full radial temperature profile in time given by Eq. (31), is determined for a 100 nm Al film on a Si substrate assuming r₀ = r₁= 15 μm and αE/f_s = αA = 5 mW when a 10 MHz modulation frequency is applied. Figure 2 depicts this scenario for a pulsed source having an 80 MHz repetition rate. Due to the description for modulation invoked, there exists a real and imaginary component to the surface temperature distribution (Fig. 2(a)), while the magnitude describes the absolute value of this temperature (Fig. 2(b)). During TDTR/FDTR experiments, the amplitude of voltage and phase lag are detected to extract an in-phase and out-of-phase voltage, which correspond to the proportional reflectivity change resulting from these real and imaginary temperature rises, respectively. That is, while the $eiω0t$ description for modulation frequency is used for mathematical convenience, the surface temperature in real time takes the form of a sinusoidal wave with some phase offset (relative to the heating event) oscillating about the steady-state temperature rise. However, because the amplitude of this oscillation dictates the depth of temperature decay, the magnitude is typically used to discuss the temperature profile associated with the TPD, as we will show later.

Note that in Fig. 2, only the AC surface temperature rise is shown, so that 0 K corresponds to the reference temperature resulting from the steady-state (DC) heating associated with the average power of the heat source. The peak temperature corresponds to a time of 100 ps to ensure validity of the thermal diffusion equation invoked here. Taking the probe-averaged temperature (that observed in TDTR) described in Eq. (32) leads to the profile shown in Fig. 2(c), while subtracting out the modulation component leads to the profile shown in Fig. 2(d). Typically, this latter temperature is related to the single-pulse (impulse) response through the real component and the pulse accumulation response through the imaginary component when a pulsed heating event is considered. For comparison, Figs. 2(c) and 2(d) also depict the solution for a CW source of equivalent average power.

Solving for the temperature and heat flux profiles for a solid in contact with a thermally conducting layer opposite the plane of irradiation follows a similar approach to that outlined in Section II A, but the boundary conditions change enough to warrant a separate discussion here. This type of geometry is typically implemented in TDTR/FDTR experiments to study solid/liquid systems,³⁷ such as that depicted in Fig. 1(b), so that we derive all equations for this example. The liquid is placed on a metal film and the laser heating source is applied to the top of the metal layer adjacent to a transparent substrate, such that the following boundary conditions apply:

where m and s subscripts denote metal and (transparent) substrate, respectively. The heat flux boundary condition implies that the laser heating is distributed to both the metal/liquid (m/l) side and the substrate (s) side, depicted in Fig. 1(b). Q_laser takes the same form of Eq. (5) as in Section II A. To be clear, in this section, the “top” refers the surface of each layer closest to z = 0 where the laser heat flux is absorbed. Thus, T_top,m is the temperature rise of the metal surface adjacent to the transparent substrate, and T_top,s is the temperature rise of the substrate surface adjacent to the metal film, both located at z = 0. Likewise, Q_top,m and Q_top,s refer to the heat fluxes at the metal and substrate surfaces, respectively. To determine T_top,m, from which we can obtain the spatially dependent temperature profile on both sides, we follow the same procedure as in the insulated case. The difference here is that two semi-infinite conditions will be invoked, one in the substrate layer and the other in the liquid layer. Proceeding first with the substrate side, the matrix solution governing the thermal transport is

Note that the negative sign in the argument of the hyperbolic terms results from defining the negative z-direction as the direction within the substrate away from the laser heating event. Applying the semi-infinite condition ensures $Q̃̃s(z)→0$ as z → –∞, so that the following relation holds:

On the metal/liquid side, the matrix solution becomes

where the l subscript denotes the liquid layer and z_l = z – d_m. Applying the semi-infinite condition ensures $Q̃̃l(z)→0$ as z → ∞, so that the following relation holds:

While expressions for T_top,m are obtained for either side, an expression in terms of a prescribed quantity, namely, Q_laser, is desired. Equations (36) and (39), together with the condition given in Eq. (34), can be solved to yield the following solution:

From this equation, we can follow the same procedure as in the insulated solution to explicitly determine the temperature profile as a function of z, r, and t for this bidirectional heat flow geometry.

We first consider an unmodulated heating event to understand how the steady-state temperature profile and TPD compare when a source is constant (CW laser) vs. time-dependent (pulsed laser). We begin with this case because it fundamentally establishes the upper limit to the attainable TPD in the modulated case. The surface temperature for CW and pulsed sources is given by Eqs. (28) and (30); we use these solutions to solve for the entire depth dependent temperature profile utilizing Eq. (13) and, when considering a layered system, Eq. (16).

For a bulk system (homogeneous, semi-infinite material) irradiated with an unmodulated CW source, the TPD is determined to be independent of both power and thermal properties of the system; only radius of the heating event is directly correlated with the TPD. Shown in Fig. 3 is the temperature profile for both bulk a-SiO₂ and bulk Si subjected to the same CW source having a heater radius r₀ = 15 μm and average absorbed power of 1 mW. In both cases, the temperature decay from the maximum temperature is equivalent when normalized. However, it is clear that a-SiO₂, given the same heat flux, has a significantly higher temperature distribution. Nonetheless, the temperature decay with depth is proportionally the same in both cases. While this may seem counterintuitive in the context of optical thermometry techniques, as the thermal perpetration depth is traditionally assumed to increase with increasing thermal conductivity, these results can be understood in the context of a simple Fourier analysis, given that the source is a constant heat flux. For example, in a one-dimensional example, given a constant heat flux boundary condition, Fourier's law implies Q_top = –κ(dT/dz), so that the temperature is described by a linear function in z, T(z) = –(Q_top /κ)z + T(0). The 1/e temperature decay, then, is simply $T(0)/e=−(Qtop/κ)zTPD+ T(0)$⁠. Solving for TPD, $zTPD=T(0)(1−1/e)(κ/Qtop)$⁠, but in accordance with Eq. (21), for a semi-infinite medium $T(0)∝(Qtop/κ)$⁠, so that the dependency of both heat flux and thermal conductivity is nullified. As a result, for simple isotropic homogeneous bulk materials, the TPD proves to be independent of the thermal properties and magnitude of heat flux. However, the heater radius becomes important in determining the TPD when extending this concept to consider the effects of radial temperature decay.

To this end, we next consider the effect of varying the pump radius. Plotted in Fig. 4(a) is the TPD of a semi-infinite bulk material as a function of pump radius for a CW gaussian source. We note that while the most common definition for TPD is the 1/e decay length, the definition invoked can make a drastic difference on the calculated thermal penetration depth. For example, changing the thermal penetration depth definition to the 1/e² or 1/e³ temperature decay length can increase this length by an order of magnitude, depending on the pump radius. No matter what the definition for TPD, however, a linear relationship is found between the pump radius and TPD.

Additionally, we can quantify the exact percent temperature decay as a function of depth for any given radius. Shown in Fig. 4(b) is this decay percentage as a function of depth for two pump radii, 1 and 100 μm. Comparing the two cases reveals that the smaller radius achieves the same decay percentage at much shallower depths than the larger radius. However, when normalizing depth to the pump radius, as defined by z/r₀, the decay percentages are exactly the same. This allows us to generalize these results to predict the decay percentage at any given depth along the r = 0 contour for any pump radius, r₀. One conclusion from these results is that, while the definition of 1/e decay length suggests a comparable TPD to 1/e² pump radius, invoking different definitions for the TPD (e.g., 99% decay length) can result in the TPD being orders of magnitude greater than this pump radius. Although the examples here are given for an unmodulated case, this idea extends similarly to the modulated TPD. Thus, interpretations of this length scale dictating physical phenomena observed in TDTR/FDTR become complicated simply by the definition of the TPD invoked.

Next, we consider the differences between pulsed and CW temperature profiles and explore the relationship such profiles hold with the repetition rate. Our calculations lead to the following general conclusions. First, the steady-state temperature profile, while independent of time in the CW case, displays a strong time dependence when considering a pulsed source. This results from the heat flux being supplied by impulses of energy such that the steady-state temperature distribution after pulse accumulation has saturated still varies with time between each pulse. Second, the solution of the steady-state temperature profile for a pulsed heating event relaxes in depth to that for a CW source of equivalent average power; the depth it takes this co-incidence of solutions to occur is dependent on the repetition rate of the system. In the limit as the repetition rate becomes slow enough that the temperature decays to zero between pulses, this co-incidence occurs only as z → 0. In this case, the solution is equivalent to that of the single pulse (impulse) temperature response. At the other extreme, as the repetition rate approaches infinity, the heating source becomes a CW source so that the thermal response relaxes to that of a CW solution very close to the surface.

To further explore the effect of the pulsed laser repetition rate on the temperature profile and quantify these concepts, we calculate the radial- and depth-dependent temperature profiles for silicon at repetition rates spanning roughly four orders of magnitude (1 MHz to 800 MHz); the results are shown in Fig. 5 for the temperature decay in depth along the r = 0 contour. For higher repetition rate pulsed sources (e.g., 800 MHz), the solution is almost indistinguishable from the CW solution at any given time. 80 MHz is a commonly used repetition rate for TDTR as it is a typical output of most commercial Ti:Sapphire oscillators. For this repetition rate, the temperature decay has a time dependence, but beyond 1 to 2 μm the depth-dependent temperature decays to be equivalent to that of the CW source solution. The temperature profile resulting from a 10 MHz repetition rate source displays similar characteristics, but takes longer in time and depth to relax to the CW solution. Still, this similarity to the 80 MHz case suggests that pulse accumulation is still present; that is, the thermal decay does not reach zero before the next pulse arrives. Finally, the 1 MHz case displays a unique trend in that the solution does not approach that of the CW source until both decay to zero. This suggests that there is little, if any, pulse accumulation. Furthermore, the t = 100 ps curve shows that before the temperature equilibrates, the magnitude of temperature is relatively low and shows a flat profile in depth beyond the near-surface region, confirming this lack of significant pulse accumulation.

We next discuss the spatially dependent temperature profiles when a modulated heat flux is applied to the surface. Since the laser heating magnitude is modulated about some average power, we note that there will always be a steady-state (DC) component to temperature rise as discussed in Section II A. In this section, however, we focus on only the modulated (AC) component to this temperature rise to comply with temperature profiles of interest to TDTR and FDTR. Equations (29) and (31) describe the surface temperature solution to such a heat flux condition. The expression given by Eq. (1) is commonly used in both TDTR and FDTR studies to correlate TPD with physical length scales associated with heat carriers.^18–23 However, this equation, given the assumptions invoked in its derivation, is limited in its practical applicability and clearly breaks down when true conditions deviate from such idealities. For example, in agreement with the conclusion drawn from the unmodulated case, the TPD depends strongly on time after pulse absorption for the pulsed case. Figure 6 demonstrates this, showing calculations of temperature rise as a function of depth for a (a) CW and (b) pulsed source. For a CW source, the temperature profile with depth at r = 0 (representing FDTR when r₁ ≪ r₀) shows that when considering the magnitude of temperature decay, Eq. (1) exactly predicts the calculated TPD for bulk a-SiO₂. However, when a pulsed response of equivalent power is observed at 100 ps, the TPD is clearly not accurately predicted by Eq. (1). Thus, the TPD varies with time in the pulsed case, making its use as a metric of measurement volume unreliable in many situations. Of course, it can be argued that since temperature deviations in pulsed vs. CW solutions occur only near-surface for high repetition rate pulsed sources (note, as per our previous discussion, this is not the case with low repetition rate sources), the experimental volume as a whole can still be reasonably predicted using the CW solution. Nonetheless, it is clear that Eq. (1) has limitations, such that its use can result in inaccurate predictions of TPD for certain experimental conditions. As it turns out, such conditions are well within the norm of those typically used in TDTR and FDTR experiments. Two such examples include: (i) the modulation frequency is on the order of the radially averaged thermal diffusivity (⁠$ω0 ≲ D/r02$⁠) and (2) layers and interfaces become significant to thermal transport. We consider both cases below.

The unmodulated temperature response suggests a physical limit to the thermal penetration depth obtainable by a modulated source; that is, the thermal penetration depth resulting from a modulated heat flux can never be greater than that resulting from an unmodulated heat flux with all other parameters staying equal. Thus, Eq. (1) breaks down at low modulation frequencies. Figure 7(a) displays the thermal penetration depth vs. modulation frequency as predicted via Eq. (1) and that calculated using the full solution discussed in the analysis above for both Si and a-SiO₂ (homogeneous semi-infinite substrates) using a pump beam radius of 15 μm. For bulk a-SiO₂, the thermal penetration depth predicted is reasonably accurate for the frequencies greater than ∼100 kHz, proving to be an effective predictor of the TPD in the range of typical modulation frequencies used in TDTR and FDTR. However, for Si, the prediction of the TPD via Eq. (1) is not accurate until ∼10 MHz using the thermal properties listed in Table I. We note that manipulating the pump radius can change the accuracy of Eq. (1) relative to the calculated results. Smaller pump radii lead to a decrease in the upper limit to the attainable TPD as calculated in the unmodulated case. As a result, the onset of disagreement between prediction and calculation occurs at higher modulation frequencies for smaller pump radii. To get a better idea of where this discrepancy begins, we normalize the TPD in a similar fashion to that of the unmodulated case (see Fig. 4) to determine the ratio of TPD/r₀ as a function of modulation frequency. Doing so reveals that the breakdown of calculation and theory, as implied by Si and a-SiO₂ examples, occurs when this ratio reaches approximately 0.2.

More importantly, however, we note that in a typical TDTR or FDTR experiment, a metal transducer (typically on the order of 100 nm) is deposited onto a bulk substrate. The addition of this thin film (and more importantly for high-thermal diffusivity materials, the interface between the thin film and substrate) serves to completely change the temperature decay profile such that for Si, the predicted TPD is never in agreement with the calculated TPD and can be off by orders of magnitude at modulation frequency extremes on both the high and low end. This is observed in Fig. 7(b), where the calculated TPD represents that within the Si or a-SiO₂ substrate after subtracting the 100 nm Al thickness. a-SiO₂, on the other hand, shows similar agreement to the predicted model as in the bulk case since the thermal resistance from the Al and Al/SiO₂ interface is negligible compared to that of the a-SiO₂.

From these observations, two generalizations can be made. First, low-thermal diffusivity materials are less prone to the breakdown of Eq. (1) within typical experimental conditions of TDTR and FDTR. Second, for high-thermal diffusivity materials, interfaces and layered structures completely change the thermal decay profile such that simple analytical expressions can no longer be used to accurately predict the TPD. Thus, to obtain accurate estimations, the full solution to the temperature decay profile with depth, radius, and time must be solved.

We now return to the exemplary case of the transparent substrate/metal film/liquid geometry supporting bidirectional heat flow (see Fig. 1(b)). This specific geometry is of interest for this study as it provides two complexities that may affect the TPDs in the semi-infinite media (i.e., liquid and substrate): (i) bidirectional heat flow and (ii) multiple layers and interfaces which can affect temperature rise profiles. We consider both the unmodulated and modulated CW heating below, noting that all the generalities and differences associated with pulsed source mentioned above remain true for this geometry.

First, the unmodulated CW source solution is depicted in Fig. 8(a), where the steady-state heating profile is shown as a function of radius and depth for both the metal/liquid and substrate sides. In this case, a-SiO₂ was used as the transparent substrate, water was used as the liquid, and a 100 nm Al film was assumed for the metal transducer. All thermal properties used in the calculation are listed in Table I. Counter-intuitively, the thermal penetration depths in both the substrate layer and metal/liquid layer are equivalent. Moreover, neglecting interfaces (i.e., assuming an infinite TBC), the temperature distribution is the same in both sides, such that there exists symmetry about z = 0.

This can be explained as follows: a constant heat flux at z = 0 results in a surface temperature rise based on the thermal properties of each side. The resulting heat flux in each side will be proportional to the thermal conductivity of that side, such that temperature decay can be understood based on a parallel thermal resistor model. As such, the heat flux is divided non-symmetrically to allow for the temperature drop to be symmetric about z = 0. To confirm this theory, in Fig. 8(b), we show the heat flux component in the z-direction as a function of depth and radius. It is clear that an asymmetry exists about z = 0 such that a higher magnitude of heat flux exists in the substrate side to compensate for its higher thermal conductivity than the liquid on the opposite side. We further discuss the role of interfaces on this temperature profile by changing the liquid used from water to air. Doing so (not shown) reveals a similar, symmetric temperature profile about z = 0 when the Al/air thermal boundary conductance is ignored and set to ∞. However, when set to a realistic value on the order of typical forced convection coefficients (∼1000 W m⁻² K⁻¹), the temperature insulation of air is preserved to reveal minimal temperature rises in the air with virtually all heat flux moving towards the substrate side. Thus, we find that the insulative nature of air stems not from its low thermal effusivity but rather its low thermal boundary conductance with most solids. In other words, heat transfer at an air-insulated boundary is limited by the convection coefficient.

Next, we consider the modulated temperature decay in solid/liquid geometries. Considering the same a-SiO₂/100 nm Al/water geometry from the unmodulated case, we find that Eq. (1) accurately predicts the TPD, within reason, for both the a-SiO₂ and water (subtracting the length of the Al layer). The only difference between calculated and predicted TPD comes from temperature drop associated with the Al/a-SiO₂ and Al/water interfaces. The temperature decays on both the substrate side and metal/liquid side are depicted in Fig. 9(a) for a 10 MHz modulated CW source with a pump radius of 15 μm and average absorbed power of 1 mW to reveal the accuracy of the TPD prediction. Because both a-SiO₂ and water have low thermal conductivities, the temperature drops at interfaces resulting from thermal boundary resistances are negligible compared to the temperature drops in the materials themselves (using the values given in Table I). As such, Eq. (1), which neglects the role of layers and interfaces, is valid within reason. However, when interfacial thermal resistances become high compared to the resistance in these materials, a clear deviation is evident. This can happen in two ways. First, the thermal conductivity of one side can become relatively large; for example, consider using sapphire (Al₂O₃), with a thermal conductivity of 35 W/m K and heat capacity of 3.06 J/cm³ K, instead of a-SiO₂ as a substrate. The thermal resistance of the Al/Al₂O₃ interface, although set to the same value as the a-SiO₂ case, is significant in this case, as shown in Fig. 9(b). Second, the thermal boundary conductance at an interface can become relatively low (for example, adding hydrophobic self assembled monolayers to the metal can significantly reduce metal/water TBC^38–40) compared to the thermal conductivity of the adjacent material to that interface. Referring again to Fig. 9(b), the Al/water TBC was reduced from 100 MW m⁻² K⁻¹ to 10 MW m⁻² K⁻¹ to show that in this case, Eq. (1) is not an accurate predictor of the TPD, despite water having such low thermal diffusivity. Thus, it is clear that relative thermal resistances determine the accuracy of Eq. (1), proving the need for rigorous modeling for accurate determination of the TPD in many cases.

Taken together, the unmodulated and modulated cases reveal the following generalizations regarding the TPD of bidirectional heat flow geometries. Much like the layered solid case with an insulative boundary condition on one side, interfaces can play a major role in determining the TPD. In the unmodulated case, while the magnitude of temperature rise is determined by all thermal parameters, the proportional temperature distribution and the TPD on both sides are independent of thermal parameters of the material, but are highly dependent on the TBC at interfaces; that is, neglecting the role of interfaces, a symmetric temperature distribution is realized on both sides no matter the difference in thermal conductivities of each side. In the modulated case, determination of the TPD depends on relative thermal resistances of interfaces and adjacent materials. When the interfacial thermal resistance is low (i.e., thermal boundary conductance is high) compared to material thermal resistances, however, we find Eq. (1) to be a valid predictor of the TPD in both the substrate and liquid, despite bidirectional heat flow. This is due to the fact that this bidirectional heat flux changes the magnitude of temperature at the sample surface (z = 0) but not the relative decay in temperature distribution within each side. Moreover, a clear difference is observed when comparing the modulated solution to the unmodulated solution. In the modulated case, the temperature distribution is not symmetric about z = 0 when neglecting interfaces. This is due to the modulation of heat flux, such that material parameters sensitive to changes in flux over the period of modulation play a role in determination of the temperature profile.

It is important to note the limitations associated with the heat diffusion framework used in this analysis. First, solutions assume negligible thermal property deviations with temperature. For small temperature changes, such as those observed in the modulated heating event given CW and high-repetition rate pulsed sources or when the heater spot size is very large, this condition is generally true. However, high temperature excursions can result in significant inaccuracies when thermal parameters are highly temperature dependent, such as when cryogenic temperatures are used. Iterative algorithms can be used whereby the solution is repeated in time and space until a convergence criterion is met to establish a temperature profile, but this proves to be difficult in practice and computationally expensive. Thus, numerical simulations such as finite element methods are likely more appropriate for such cases.

Another limitation we note is that the solution relies on a surface heating event. This is, for a wide variety of applications, a reasonable assumption. For example, in TDTR and FDTR, a metal transducer is deposited onto the sample of interest so that the optical penetration depth in most cases is on the order of ∼10 nm, much less than the thermal penetration depth of most materials. However, for bulk semiconductors and insulators, this optical penetration depth can be much larger, such that surface heating may not accurately describe the heat flux provided by laser interaction. Modification to the heat transfer analysis in this case has been discussed previously.^24,25

The final limitation of this model is that it does not consider the physical processes associated with non-equilibrium interactions that occur during laser absorption before thermal diffusion is valid; such interactions include electron-phonon, magnon-phonon, and exciton-phonon coupling.^41–48 These interactions typically occur before thermal diffusion can be described or are limited in their contribution to thermal transport, which is why in this paper we restrict all time scales to a minimum of 100 ps to ensure a thermal diffusion time regime for all calculations. Still, because of ballistic electron length scales and other possible carrier interactions having longer time-scales, the heating event described may not be valid in certain materials. Inclusion of these more complex heat sources is a future direction of study.

In conclusion, in this study we have developed a method to solve for the depth- and radially dependent temperature decay for bulk, homogeneous materials and layered thin-films subjected to pulsed and CW laser heating, both modulated and unmodulated. We also include analyses of geometries with conductive boundary conditions on either side of a deposited heat flux (e.g., substrate/metal film/liquid geometries). We find that in the unmodulated CW case, the thermal penetration depth, defined as the depth at which the surface temperature decays to its 1/e value, is dependent only on the radius of the Gaussian heating source, while an additional temporal dependence is observed when considering a pulsed laser source. In particular, as the pulsed laser repetition rate approaches infinity, the temperature rise profile approaches that of the CW profile. At the other extreme, as the repetition rate is reduced to a regime of negligible pulse accumulation, the thermal penetration depth becomes ill-defined without some equilibrium temperature distribution to which the decay can settle in time. The solution to the modulated heating event reveals that for a CW source, the typically used expression for predicting TPD (Eq. (1)) works reasonably well for low-thermal diffusivity materials (e.g., a-SiO₂) but breaks down for high-thermal diffusivity materials (e.g., Si). Of course, the exact point of failure depends on specific experimental parameters. Yet, we find that such parameters are the norm for most TDTR and FDTR experiments, leading to overestimations of TPD by, in some cases, orders of magnitude. For the pulsed case, these conclusions remain the same and are even exacerbated by the time dependence of TPD for a pulsed source. In the case of bidirectional heat flow, we find that the steady-state temperature distribution resulting from an unmodulated source is symmetric about the plane of laser absorption no matter the thermal properties intrinsic to the materials on each side due to an asymmetry of heat flux about this plane. However, interfacial thermal resistance is responsible for breaking this symmetry in temperature; for example, when air is used as the liquid layer, the thermal insulation is entirely due to the convection coefficient between metal and air rather than the low thermal effusivity of air. This study is therefore important to applications relying on dissipation of heat, both modulated and unmodulated, as well as optical pump-probe experiments for relating length scales of physical phenomena to thermal penetration length scales observed experimentally.

This work was funded through the Office of Naval Research, Grant No. N00014-15-12769. This research was conducted with Government support under and awarded by DoD, Air Force Office of Scientific Research, National Defense Science and Engineering Graduate (NDSEG) Fellowship, 32 CFR 168a.