We previously demonstrated an extension of time-domain thermoreflectance (TDTR) which utilizes offset pump and probe laser locations to measure in-plane thermal transport properties of multilayers. However, the technique was limited to systems of transversely isotropic materials studied using axisymmetric laser intensities. Here, we extend the mathematics so that data reduction can be performed on non-transversely isotropic systems. An analytic solution of the diffusion equation for an N-layer system is given, where each layer has a homogenous but otherwise arbitrary thermal conductivity tensor and the illuminating spots have arbitrary intensity profiles. As a demonstration, we use both TDTR and time-resolved magneto-optic Kerr effect measurements to obtain thermal conductivity tensor elements of <110> α-SiO2. We show that the out-of-phase beam offset sweep has full-width half-maxima that contains nearly independent sensitivity to the in-plane thermal conductivity corresponding to the scanning direction. Also, we demonstrate a Nb-V alloy as a low thermal conductivity TDTR transducer layer that helps improve the accuracy of in-plane measurements.

Time domain thermoreflectance (TDTR) is a pump-probe technique that is used to measure the thermal properties of materials and interfaces. TDTR utilizes surface heating of a material via a train of laser pulses (pump beam) and monitors the resulting temperature change through the reflectance of the surface, sampled using a time-delayed sensing beam. The train of pump pulses is modulated, which is used for lock-in detection of the thermoreflectance signal and generates useful pulse accumulation effects.1,2 The physical implementation of our two-tint TDTR system and the data reduction for layered materials has been described previously.2,3 TDTR has been applied to measure the thermal transport properties of a wide range of thin-films,4 multilayers,5 bulk materials, and their interfaces.6 Recently, we have shown that when TDTR is performed using spatially offset pumping/sensing beams, the resulting signal gives information about the in-plane heat spreading.7 

To extract thermal properties from TDTR experiments, it is necessary to compare the data to an accurate thermal simulation. However, previous descriptions of heat flow in TDTR have been limited to transversely isotropic materials studied with axisymmetric laser intensities (with the notable exception of spatially offset Gaussian beams used to measure radially symmetric materials1,7). There are several common situations of importance which are not captured by this model: (1) materials that are transversely isotropic, such as hexagonal or tetragonal crystals, but are not cut perpendicular to their c-axis (2) illumination with non-axisymmetric beams; such illumination includes elliptical spots (typical of aberrations in lasers) and complex patterns arising from gratings. (3) Crystals and metamaterials with lower symmetry also exist, and require the use of a more complete model. Common examples include orthorhombic materials such as gallium.

In the current work, we derive an analytic model for the TDTR signals obtained from an N-layer multilayer system, where each layer has a homogeneous but otherwise arbitrary thermal conductivity tensor and the illuminating spots have arbitrary intensity profiles. We show that by performing offset-beam measurements along multiple axes, the in-plane components of the thermal conductivity tensor can be determined independently to a good approximation. We demonstrate the method by performing TDTR and TR-MOKE (time-resolved magneto-optic Kerr effect) experiments on ⟨110⟩ α-SiO2 (“X-Cut quartz”).

When temperature fluctuations are small such that the heat transport equations are linear and time-invariant, the thermoreflectance signal measured during a TDTR experiment can be written in terms of the frequency-domain response of the thermal system:2 

\begin{eqnarray}{\mathop{\rm Re}\nolimits} [\Delta R_M (t_d)] &=& \frac{{dR}}{{dT}}\sum\limits_{m = - \infty }^\infty ( \Delta T(mf_{{\rm rep}} + f_{{\rm mod}})\nonumber\\[-2pt]&& + \Delta T(mf_{{\rm rep}} - f_{\bmod }))\exp (im2\pi f_{{\rm rep}} t_d),\nonumber\\\end{eqnarray}
Re [ΔRM(td)]=dRdTm=(ΔT(mf rep +f mod )+ΔT(mf rep fmod))exp(im2πf rep td),
(1)
\begin{eqnarray}{\mathop{\rm Im}\nolimits} [\Delta R_M (t_d)] &=& - i\frac{{dR}}{{dT}}\sum\limits_{m = - \infty }^\infty ( \Delta T(mf_{{\rm rep}} + f_{{\rm mod}})\nonumber\\[-2pt]&& - \Delta T(mf_{{\rm rep}} - f_{{\rm mod}}))\exp (im2\pi f_{{\rm rep}} t_d),\nonumber\\\end{eqnarray}
Im [ΔRM(td)]=idRdTm=(ΔT(mf rep +f mod )ΔT(mf rep f mod ))exp(im2πf rep td),
(2)

where frep, fmod, and td are the laser repetition rate (in Hz), modulation frequency, and time delay, respectively. The frequency-domain response ΔT(f) is defined as the average temperature fluctuation sampled by the sensing beam, due to a continuous pump beam that is modulated sinusoidally in time at frequency f; the signal measured during a frequency domain thermoreflectance experiment8 is proportional to ΔT(f). In the case of TDTR systems where the mechanical stage advances the arrival time of the pump beam rather than delays the probe, there is an additional phase factor ΔR(td) = ΔRM(td)exp (ifmodtd). The real and imaginary parts of the reflectivity fluctuations are recorded as the “in-phase” and “out-of-phase” signals, Vin and Vout, respectively. In conventional TDTR, “the ratio,” −Vin/Vout is fit to a thermal simulation used to extract the thermal properties from the signal.2 

Equations (1) and (2) hold for any similar transient optical experiment where the temperature is measured by a train of sensing pulses that produces a signal linearly proportional to the surface temperature fluctuations. For example, in time-resolved magneto-optic Kerr effect (TR-MOKE) measurements, the signal is the fluctuation of the Kerr rotation (Δθ ∝ ΔM) instead of the fluctuation in reflectance (ΔR), and the sensing coefficient (dR/dT) in Eqs. (1) and (2) is replaced by the thermo-magneto-optic Kerr rotation coefficient, (dθ / dT = (dθ / dM)(dM / dT)).

If the pumping and sensing beam intensities are given by p(x, y) and s(x, y), then the intensity-weighted temperature fluctuation detected by the sensing beam is

\begin{equation}\Delta T(f) = \frac{1}{{A_S }}\int\limits_\xi\!\!\! {\int\limits_\eta {\hat G(f,\xi,\eta)\hat P(\xi,\eta)\hat S( - \xi, - \eta)d\xi d\eta } },\end{equation}
ΔT(f)=1ASξηĜ(f,ξ,η)P̂(ξ,η)Ŝ(ξ,η)dξdη,
(3)

where |$\hat P(\xi,\eta)$|P̂(ξ,η) and |$\hat S(\xi,\eta)$|Ŝ(ξ,η) are the Fourier transforms of the pump and sensing intensity, AS is the total power of the sensing beam, and |$\hat G(f,\xi,\eta)$|Ĝ(f,ξ,η) is the Green's function solution for the temperature response due to a sinusoidally oscillating point source of unit strength.2 We obtained the solution, |$\hat G(f,\xi,\eta)$|Ĝ(f,ξ,η) using the Feldman transfer matrix technique.9 The derivation is given as the Appendix. We allow the thermal conductivity to be a homogenous 2nd rank tensor written in Cartesian coordinates with local heat flux, q, given, using Einstein index notation, by

\begin{equation}q_i = - \kappa _{ij} \frac{{\partial T}}{{\partial x_j }}.\end{equation}
qi=κijTxj.
(4)

The second law of thermodynamics requires that the symmetric part of κij to be positive semi-definite. Although there is no rigorous proof that the anti-symmetric part of κij must vanish, two consequence of a non-zero anti-symmetric part are that (1) heat flow can occur without entropy production (this does not violate the 2nd law of thermodynamics) and (2) temperature disturbances can be maintained indefinitely without the need for an external heat flux.10 Since these two points are counter to observation, we consider only a symmetric matrix with up to 6 independent components defined as

\begin{equation}\kappa _{ij} = \left[ {\begin{array}{l@{\quad}c@{\quad}c}{\kappa _{xx} } & {\kappa _{xy} } & {\kappa _{xz} } \\{\kappa _{xy} } & {\kappa _{yy} } & {\kappa _{yz} } \\{\kappa _{xz} } & {\kappa _{yz} } & {\kappa _{zz} } \\\end{array}} \right],\end{equation}
κij=κxxκxyκxzκxyκyyκyzκxzκyzκzz,
(5)

where z corresponds to the through-plane direction. Physically, the interpretation is that there may be up to 3 positive, independent, and orthogonal principal thermal conductivities and 3 angles that describe their orientation with respect to the laboratory coordinates.

The intensity profile of the pump and sensing beam can be arbitrarily specified by providing their Fourier transforms, |$\hat P(\xi,\eta)$|P̂(ξ,η) and |$\hat S(\xi,\eta)$|Ŝ(ξ,η), for use in Eq. (3), but in the current work we will consider only elliptical Gaussian beams with major and minor axis aligned with Cartesian coordinate axes x and y, which are the typical output of lasers, where the pump is centered on the origin and the sensing beam may be offset in the x, y direction by x0 and y0, respectively. Then, the corresponding transforms are

\begin{equation}\hat P(\xi,\eta) = \frac{{A_P }}{{4\pi }}\exp \big[ { - \big( {\xi ^2 w_{x,0}^2 + \eta ^2 w_{y,0}^2 }\big)/8}\big],\end{equation}
P̂(ξ,η)=AP4πexpξ2wx,02+η2wy,02/8,
(6)
\begin{eqnarray}\hat S(\xi,\eta) &=& \frac{{A_S }}{{4\pi }}\exp \big[ - \big( {\xi ^2 w_{x,0}^2 + \eta ^2 w_{y,0}^2 }\big)/8 \big]\nonumber\\&&\times\exp [ {i(\xi x_0 + \eta y_0)}].\end{eqnarray}
Ŝ(ξ,η)=AS4πexpξ2wx,02+η2wy,02/8×exp[i(ξx0+ηy0)].
(7)

Using Eqs. (6), (7), and (A23) the frequency domain response (Eq. (3)) can be calculated, from which the experimental signal can be calculated using Eqs. (1) and (2). Note that in this work, due to the use of Eqs. (6) and (7), we only consider elliptical spots whose major and minor axis are aligned with the beam offset direction, since this closely corresponds to our experimental situation (a common aberration in Ti:Sapphire ultrafast lasers). For other illumination scenarios, including rotated ellipses and grating profiles, the only change that would be needed would be to use the relevant Fourier transforms, |$\hat P(\xi,\eta)$|P̂(ξ,η) and |$\hat S(\xi,\eta)$|Ŝ(ξ,η).

In the case of transversely isotropic material properties with concentric spots, we recover the solution given by Cahill,2 and in the case of transverse isotropy with offset beams we have verified numerically that the results are identical to those of Schmidt1 and Feser.7 The solution in the current work has the practical disadvantage that simulations require about 2 orders of magnitude more computational time than Ref. 7 due to the multidimensional integral in Eq. (3), which suggests that its use is not desirable except for systems which lack transverse isotropy. Using our code implementation, we find that to generate a simulation of ΔR(td) required ∼10–1000 s on a personal computer (6-core Intel Xeon processor at 3.2 GHz using MATLAB and the Parallel Computing Toolbox), depending on the degree of heat spreading and the beam offset distance. Generally, larger beam offset distances require larger simulation times due to the fast oscillations that arise from S(ξ, η) in Eq. (7) that must be evaluated to obtain the integral in Eq. (3).

In the high frequency limit corresponding to 1D through-plane transport, the only transport property affecting the frequency domain solution is the through-plane thermal conductivity, κzz. The primary utility of the general-tensor frequency-domain formulation of this paper is the ability to analyze experiments where in-plane heat spreading in non-negligible, and either the anisotropy of the material or illumination is not axisymmetric (or reducible to an equivalent axisymmetric system by the method given in Ref. 7).

We have recently shown that for experiments with transversely isotropic materials, offsetting the pump and probe beam in TDTR provides useful information about in-plane heat spreading.7 In particular, the in-phase and out-of-phase signals scale differently with the offset distance such that the out-of-phase signal decays more slowly with offset, due to the fact that the physical origin of the out-of-phase signal is pulse accumulation whereas the in-phase signal is approximately interpreted as instantaneous pulse absorption/decay (this approximation is valid at high modulation frequency, ffrep/8 and for short, positive time-delay). Thus, the in-phase signal nearly traces the cross-correlation between the pumping/sensing beam intensities, and the out-of-phase signal traces the temperature field associated with the modulation frequency, which has broader width if significant heat spreading exists. By using the in-phase signal to precisely measure the beam sizes, and by fitting the ratio between the in-phase and out-of-phase signals, we showed that in-plane properties of thin films and bulk materials can be measured using TDTR.

Using the analytical model for TDTR presented in Section II, the beam-offset technique can be extended to measure the in-plane thermal conductivity along multiple in-plane axes, independently. This is accomplished by sweeping the offset between the pumping and sensing beams using the apparatus described in Ref. 7, which can control the offset with a precision of ∼30 nm when using a 1 μm spot size (the precision scales inversely with the objective magnification). In this paper, we also demonstrate that the technique is compatible with TR-MOKE detection. The beam-steering optics for TR-MOKE is identical, but the detection scheme is altered as follows11: using an initially linearly polarized probe beam, we detect fluctuations in the polarization angle that results from reflectance of the probe from a perpendicularly magnetized transducer multilayer of CoPt (detailed structure is described later). To detect the polarization rotation, a non-polarizing cube beam splitter is inserted between the steering polarizing beam splitter and the microscope objective lens which diverts the reflected pump and probe beams toward a detection path. In the detection path, a sharp optical filter (Semrock, Razor Edge, short-pass, 129 cm−1 transition) removes the pump beam, and the polarization fluctuations are monitored by splitting the probe into orthogonal polarization states with a Wollaston prism (Thorlabs WP10-B) and detecting the changes in the intensity using a balanced detector. A half-wave plate placed before the Wollaston prism is used to balance the average intensity on each photodiode.

Rather than analyzing the ratio as is done with conventional TDTR,2 we have recently found that when extracting in-plane thermal properties it is often better to measure and analyze the full-width half maximum (FWHM) of the out-of-phase signal with respect to offset distance at negative delay time in order to determine the in-plane thermal conductivities. There are several reasons why this is advantageous: (1) Similar to the ratio, the FWHM is independent of amplitude of the temperature fluctuations for small excursions. (2) The FWHM in any particular offset direction is more sensitive to the in-plane thermal conductivity along that direction compared to the ratio. (3) The FWHM is nearly independent of thermal conduction occurring in orthogonal directions. This is not usually the case if analyzing the ratio. We have shown previously that, even if the material is only transversely isotropic, the through-plane thermal conductivity must generally be well known in order to extract the in-plane thermal conductivity using the ratio.7 This is not the case if the FWHM is used as will be shown. (4) For cases where in-plane heat penetration depth is large compared to the spot size, the FWHM along different principal axes shows more anisotropy than the ratio.

The reason for the last point can be illustrated by appealing to the origin of the in-phase and out-of-phase signals. At positive time delay, the in-phase signal is approximately related to the cross-correlation of the pump and probe beam,7 and thus has form |$V_{in} \sim\exp ( - x_0^2 /w_0^2)$|Vinexp(x02/w02) which is independent of the offset direction for a Gaussian spot. In the high frequency limit, the out-of-phase signal has exactly the same shape, |$V_{out} \sim\exp ( - x_0^2 /w_0^2)$|Voutexp(x02/w02), and thus neither the shape of the ratio nor the FWHM of Vout contains any information about heat spreading. However, if heat spreading is comparable or larger than the size of the spot, |$V_{out} \sim\exp ( - x_0^2 /L^2)$|Voutexp(x02/L2), where L is a length scale that is larger than the spot size and depends on the in-plane thermal penetration depth. For |$L^2 \gg w_0^2 $|L2w02, the ratio approaches |$ - V_{in} /V_{out} \sim\exp ( - x_0^2 /w_0^2)$|Vin/Voutexp(x02/w02) and, so again, the shape of the ratio does not contain any information about the heat spreading (i.e., the shape is always isotropic with respect to beam offset). On the other hand, Vout does have information about the thermal penetration depth in the offset direction, and thus the FWHM of Vout is effectively a direct measurement of the thermal penetration depth in that direction if |$L^2 \gg w_0^2 $|L2w02. In the intermediate regime, both the ratio and the FWHM have information about heat spreading, but since the FWHM of Vout is more robust, we adopt that approach for determining in-plane thermal conductivities. By measuring the penetration depth along multiple offset directions, values and the orientation of the in-plane principal axes of thermal conductivity can be ascertained. To our knowledge, it is not possible to perform the experiment such that the values of κxz and κyz can be obtained independently, although κzz can be easily obtained by independently performing experiments using large spot size and high-modulation frequency (i.e., by probing 1D through-plane heating). The current model does have the ability to model the effect of non-zero κxz and κyz if their values can be inferred from material symmetry and knowledge of crystal orientation, as is often the case.

As an example of how this type of analysis can be used to measure in-plane properties of samples that lack in-plane symmetry, experiments were performed on ⟨110⟩ α-SiO2 purchased from MTI Corp (optical grade, X-cut quartz) coated with a 63 nm thick film of a Nb0.43V0.57 alloy for TDTR, or a 19 nm CoPt multilayer for TR-MOKE experiments, which served as transducers. Quartz is a trigonal crystal and has higher thermal conductivity parallel to its c-axis than in the plane perpendicular to it. In X-cut quartz, the c-axis lies in-plane (Fig. 1). To simplify the analysis, the sample was aligned using XRD such that the two orthogonal beam offset sweeps corresponded to: (1) an x-offset perpendicular to the c-axis of the crystal and (2) an y-offset parallel to the c-axis. A 50× magnification objective was used to produce beams with ∼1 μm spot size. The precise size of the spots was measured in each sweep direction by probing at high frequency (9.8 MHz) and small positive delay time (100 ps); the result is shown in Fig. 2, and the spot size is extracted by fitting Vin to a Gaussian, |$V_{in} \sim\exp ( - x_0^2 /w_{0,i}^2)$|Vinexp(x02/w0,i2) in each direction. Then, the beam offset sweep was performed again, but at lower frequency (1.11 MHz) and at short negative delay time (−100 or −20 ps). Performing the experiment at negative delay time reduces the amplitude of Vin; reducing Vin prevents phase jitter or small errors in the lockin-amplifier reference phase from altering Vout. Performing the experiment at low frequency increases the in-plane thermal penetration depth and consequently the sensitivity of the experiment. We use a heat capacity of α-SiO2 from literature12 of Cα-SiO2 = 2.02 MJ/m3-K.

FIG. 1.

Alignment of the <110> SiO2 principal thermal conductivity axes with the experimental coordinate system: κc > κa.

FIG. 1.

Alignment of the <110> SiO2 principal thermal conductivity axes with the experimental coordinate system: κc > κa.

Close modal
FIG. 2.

Experimental beam offset data for <110> α-SiO2. Vin at positive time delay and high frequency used to measure the spot size for (a) TDTR and (c) TR-MOKE methods. Vout at negative delay time and low frequency, used to extract FWHM for (b) TDTR and (d) TR-MOKE methods. Closed- and open-circles correspond to beam offset scans along the in-plane direction perpendicular to and parallel to the c-axis, respectively.

FIG. 2.

Experimental beam offset data for <110> α-SiO2. Vin at positive time delay and high frequency used to measure the spot size for (a) TDTR and (c) TR-MOKE methods. Vout at negative delay time and low frequency, used to extract FWHM for (b) TDTR and (d) TR-MOKE methods. Closed- and open-circles correspond to beam offset scans along the in-plane direction perpendicular to and parallel to the c-axis, respectively.

Close modal

The sensitivity values of the experiment are shown in Fig. 3, and help to illustrate the general characteristics and advantages of the FWHM measurement technique. The sensitivity here is defined as

\begin{equation}S_i^{FWHM} \equiv \frac{{\partial [\ln (FWHM)]}}{{\partial [\ln (\alpha _i)]}}\end{equation}
SiFWHM[ln(FWHM)][ln(αi)]
(8)

and represents the relative change in the FWHM for a relative change in parameter, α. In Fig. 3, we consider the thermal conductivities (κ), volumetric heat capacity (C), film thickness (h), and thermal interface conductance, (G), and x,y direction 1/e2 spot radii (w) as the model input parameters. Figure 3 shows that the sensitivity is small for all substrate κij with the exception of the thermal conductivity along the offset direction. Furthermore, Fig. 3 shows that the FWHM is more sensitive to spot size than any other parameter, and thus spot size must be known precisely in order to measure in-plane thermal conductivity; we do this by measuring the beam offset shape at high modulation frequency and at short positive time-delay. If the spot shape is elliptical, to a good approximation the extracted thermal conductivity only depends on the axis of the ellipse aligned along the offset direction. Note that the sensitivity has equal but oppositely signed values for the in-plane thermal conductivities and heat capacity; this indicates that the measurement is actually extracting the thermal diffusivity, κ/C, along the offset direction. This is in contrast to conventional TDTR which extracts the thermal effusivity, |$\sqrt {\kappa C} $|κC, and is not surprising since the current measurement is essentially a direct measurement of the in-plane thermal penetration depth, |$L\sim\sqrt {(\kappa /(Cf)} $|L(κ/(Cf). It is not possible to define a meaningful sensitivity for the cross-component of the thermal conductivity tensor, κxz and κyz, since they are zero in the current experiment; however, the derivatives ∂(FWHMx)/∂κxz, ∂(FWHMy)/∂κxz, ∂(FWHMx)/∂κyz, and ∂(FWHMy)/∂κyz are still well-defined and can be shown to be identically zero for κxz = κyz = 0 as in the current experiment which indicates that, at least for small values, they cannot be measured via the FWHM method.

FIG. 3.

Sensitivity values for the (a) TDTR and (b) TR-MOKE α-SiO2 beam offset experiment described in the text. Sensitivities are of the FWHM perpendicular (solid) and parallel to the c-axis (open). The first four parameters are properties of the substrate; the subscript M stands for “metal,” i.e., the transducer layer properties.

FIG. 3.

Sensitivity values for the (a) TDTR and (b) TR-MOKE α-SiO2 beam offset experiment described in the text. Sensitivities are of the FWHM perpendicular (solid) and parallel to the c-axis (open). The first four parameters are properties of the substrate; the subscript M stands for “metal,” i.e., the transducer layer properties.

Close modal

In order to measure the low thermal conductivity of α-SiO2 and avoid loss in sensitivity due to in-plane heat spreading within the metal transducer, we found it necessary to develop a new, low thermal conductivity TDTR transducer, NbV alloy (Nb0.43V0.57 as measured by Rutherford backscattering spectrometry). The thermal conductivity of NbV is between 15 and 22 W/m-K depending on thickness (15 W/m-K @ 32 nm, 18 W/m-K @ 65 nm as measured by the beam offset technique described in Ref. 7 and 22 W/m-K @ 500 nm as measured through-plane by conventional TDTR), which is low enough to not dominate the in-plane thermal conductance of the current experiment (Fig. 3). Note that it is possible that the through-plane thermal conductivity is different than the in-plane thermal conductivity due to columnar grain growth, and may not truly represent thickness dependence. The use of a thin, low thermal conductivity transducer is important because the sensitivity of the measurement to the conductance of the transducer increases nearly linear with transducer conductance (κMhM). Therefore, if a 70 nm aluminum (κ ∼ 200 W/m-K) or gold (κ ∼ 300 W/m-K) transducer was used, the measurement would become ∼20–30 times more sensitive to the properties of the transducer relative to the thermal conductivity of the substrate. Consequently, the measurement sensitivity would become dominated by the transducer, and an accurate measurement would cease to be possible.

NbV was found to have a moderate thermoreflectance coefficient (dR/dT ∼ 2 × 10−5 at λ = 785 nm), high absorption coefficient (∼0.37), and a clear picosecond acoustic response. We measured a speed-of-sound, v[110] = 5.40 nm/ps, using picosecond acoustic reflections with thickness obtained independently using X-ray reflectivity. Sputtered NbV alloy was found to deposit with strong [110] texturing as determined by X-ray diffraction, which is typical for bcc metals. We experimentally obtained the heat capacity, C = 2.65 ± 0.13 MJ/m3-K at room temperature by performing TDTR on a sample with known thermal conductivity (amorphous thermally grown SiO2) and known transducer thickness obtained from XRR. The Dulong-Petit limit of the lattice heat capacity can be computed using the lattice constant as Clattice, DP = 2.55 MJ/m3-K, which compares favorably with the measured value but is slightly lower. Note the Nb and V elements have relatively high predicted electronic heat capacities;13 Cel = γT yields Cel, Nb = 0.22 MJ/m3-K and Cel, V = 0.35 MJ/m3-K at 300 K, which could explain why the measured heat capacity is larger than the Dulong-Petit limit.

The TR-MOKE transducer used was CoPt multilayer with a total thickness of 19 nm as measured by X-ray reflectivity; the nominal target structure, from top to bottom, was Pt(1 nm)/[Co(0.5 nm) Pt(1 nm)] × 6/Pt(10 nm). The estimated heat capacity is 2.93 J/cm3-K and the in-plane thermal conductivity, estimated from the electronic part of electrical conductivity using the Weidemann-Franz law with degenerate Lorentz number, was 26.7 W/m-K. Note that measurements using the MOKE transducer are less affected by the properties of the transducer (Fig. 3(b)), because of its small thickness and low thermal conductivity. In addition, the temperature dependence of polarization rotation to transducer magnetization is larger than effects arising from the temperature-dependent optical indices of the non-magnetic sublayers; thus, semitransparency of the transducer layer is more tolerable in TR-MOKE compared to TDTR, where the thermoreflectance may be strongly influenced by light interaction with sub-layers if the transducer is semi-transparent. In that sense, the TR-MOKE technique is well suited for measurements of in-plane properties.

To extract the thermal conductivity from the beam offset experiment, simulations of Vout vs. beam offset, (x0, y0) were performed, and the FWHMs of the resulting profiles were tabulated as function of the thermal conductivity. By matching the measured FWHM to the simulated FWHM, the corresponding thermal conductivity was found (Fig. 4). For TDTR, the extracted thermal conductivity in the case of <110> α-SiO2 is κa = 6.6 W/m-K (Fig. 4(a)) and κc = 11.0 W/m-K (Fig. 4(b)) along the directions perpendicular to and parallel to the c-axis, respectively. For the TR-MOKE experiment, the corresponding values are κa = 6.9 W/m-K (Fig. 4(c)) and κc = 11.0 W/m-K (Fig. 4(d)). These are in reasonable agreement with the recommended values for α-SiO2 (κc = 10.4 W/m-K and κa = 6.2 W/m-K).14 Note that the extracted a-direction thermal conductivity is obtained nearly independently of the through-plane thermal conductivity in this approach (see Fig. 3). We have found that the experimental setup has a FWHM repeatability of <50 nm for experiments with a FWHM of ∼3 μm (or <1%–2% uncertainty). This repeatability is only slightly larger than the minimum incremental motion of our apparatus for the given microscope objective, ∼34 nm. Thus, although the sensitivity to in-plane conductivity ∂[ln (FWHM)]/∂[ln (kxx)] ∼ 0.2 (see Fig. 3), that is large enough that the error due to uncertainty in the FWHM is only 5%–10% which is comparable to other errors in the experiment, such those due to uncertainty in the spot size. We have visually indicated this in Fig. 4 using dashed lines to represent the range of FWHM uncertainty.

FIG. 4.

Comparison of the experimental FWHM (thick solid lines) with simulated FWHM as function of the in-plane thermal conductivities, used to extract the thermal properties. The dashed lines are the uncertainty in the FWHM estimated from the repeatability of the measurement (∼0.05 μm). Simulations used nominal κz = 5.5 W/m-K (fixed), κc = 10.4 W/m-K (left graph), and κa = 6.2 W/m-K (right graph).

FIG. 4.

Comparison of the experimental FWHM (thick solid lines) with simulated FWHM as function of the in-plane thermal conductivities, used to extract the thermal properties. The dashed lines are the uncertainty in the FWHM estimated from the repeatability of the measurement (∼0.05 μm). Simulations used nominal κz = 5.5 W/m-K (fixed), κc = 10.4 W/m-K (left graph), and κa = 6.2 W/m-K (right graph).

Close modal

By crystal symmetry, the through-plane thermal conductivity should be the same as the extracted κa; however, we have found that conventional TDTR using large spot sizes (w = 10 μm) and high frequency (f = 9.8 MHz) consistently gives through-plane thermal conductivity values that are slightly lower than the recommended value (κzz = κa = 5.5 ± 0.5 W/m-K). Previous work has shown that anomalously low thermal conductivities extracted from high-modulation rate or high-spatial frequency measurements may be attributed to ballistic-phonon effects.15–18 Given the relatively low bulk thermal conductivity (and thus mean-free-path) of α-SiO2, this would seem to be a somewhat surprising observation.

If the in-plane sample orientation is not known before the experiment, it is possible to obtain it. By mapping Vout(x0, y0), the in-plane principal axes can be identified as the major and minor axes of the resulting elliptical contours, and the FWHM of those axes can be used to determine the value of the in-plane principal thermal conductivities. An example map is shown for the same quartz sample using TDTR, but with the c-axis intentionally rotated ∼ 120° relative to the x-axis. Although the thermal anisotropy is relatively small, the orientation can be easily discerned, and the FWHM of the ellipse can be estimated using interpolating contour software (we used IgorPro to numerically generate the contour in Fig. 5).

FIG. 5.

Experimental beam offset map of Vout for a non-aligned α-SiO2 sample (c-axis rotated ∼120° from x-axis). The FWHM ellipsoid is shown as a thick black line.

FIG. 5.

Experimental beam offset map of Vout for a non-aligned α-SiO2 sample (c-axis rotated ∼120° from x-axis). The FWHM ellipsoid is shown as a thick black line.

Close modal

We have shown that beam-offset pump-probe measurements, coupled with a mathematical model that considers general anisotropy and surface intensity profiles, can be used to measure elements of the thermal conductivity tensor whenever the spot-sizes are comparable or smaller to the in-plane diffusion of heat. We have shown that measuring the FWHM of the beam-offset sweep at negative time-delays provides a convenient method to isolate thermal conductivity tensor components from one another. To enhance the sensitivity of TDTR to in-plane substrate thermal conductivities, we have demonstrated NbV alloy as a low thermal conductivity thermal transducer. However, since TR-MOKE experiments can utilize a much thinner (and low thermal conductivity) transducer, we have found that it also provides sufficient sensitivity to in-plane properties for the purposes of measuring anisotropic sample properties.

This work was supported by the Army Research Office Grant Nos.W911NF-11-10526 and W911NF-14-1-0016. TDTR data were acquired using the equipment in the Laser Facility of the Frederick Seitz Materials Research Laboratory (MRL) at the University of Illinois at Urbana-Champaign.

Assuming homogeneous, but anisotropic material properties, the general heat transport equation for each layer is

\begin{eqnarray}C\frac{{\partial T}}{{\partial t}} &=& k_{xx} \frac{{\partial ^2 T}}{{\partial x^2 }} + k_{yy} \frac{{\partial ^2 T}}{{\partial y^2 }} + k_{zz} \frac{{\partial ^2 T}}{{\partial z^2 }} + 2k_{xy} \frac{{\partial ^2 T}}{{\partial x\partial y}}\nonumber\\&& + 2k_{xz} \frac{{\partial ^2 T}}{{\partial x\partial z}} + 2k_{yz} \frac{{\partial ^2 T}}{{\partial y\partial z}}.\end{eqnarray}
CTt=kxx2Tx2+kyy2Ty2+kzz2Tz2+2kxy2Txy+2kxz2Txz+2kyz2Tyz.
(A1)

Taking the Fourier transform with respect to the in-plane coordinates and time, |$T(x,y,z,t) \leftrightarrow \hat T(\xi,\eta,z,\omega)$|T(x,y,z,t)T̂(ξ,η,z,ω)

\begin{eqnarray}\left( {\frac{{iC\omega }}{{k_{zz} }}} \right)\hat T &=& - \left[ {\frac{{k_{xx} }}{{k_{zz} }}\xi ^2 + 2\frac{{k_{xy} }}{{k_{zz} }}\xi \eta + \frac{{k_{yy} }}{{k_{zz} }}\eta ^2 } \right]\hat T\nonumber\\&& +\, 2i\left[ {\frac{{k_{xz} }}{{k_{zz} }}\xi + \frac{{k_{yz} }}{{k_{zz} }}\eta } \right]\frac{{\partial \hat T}}{{\partial z}} + \frac{{\partial ^2 \hat T}}{{\partial z^2 }},\end{eqnarray}
iCωkzzT̂=kxxkzzξ2+2kxykzzξη+kyykzzη2T̂+2ikxzkzzξ+kyzkzzηT̂z+2T̂z2,
(A2)

or more compactly,

\begin{equation}\frac{{\partial ^2 \hat T}}{{\partial z^2 }} + \lambda _2 \frac{{\partial \hat T}}{{\partial z}} - \lambda _1 \hat T = 0,\end{equation}
2T̂z2+λ2T̂zλ1T̂=0,
(A3)

where

\begin{equation}\lambda _1 \equiv \frac{{i\omega C}}{{k_{zz} }} + \left[ {\frac{{k_{xx} }}{{k_{zz} }}\xi ^2 + 2\frac{{k_{xy} }}{{k_{zz} }}\xi \eta + \frac{{k_{yy} }}{{k_{zz} }}\eta ^2 } \right],\end{equation}
λ1iωCkzz+kxxkzzξ2+2kxykzzξη+kyykzzη2,
(A4)
\begin{equation}\lambda _2 \equiv 2i\left[ {\frac{{k_{xz} }}{{k_{zz} }}\xi + \frac{{k_{yz} }}{{k_{zz} }}\eta } \right].\end{equation}
λ22ikxzkzzξ+kyzkzzη.
(A5)

Equation (A3) has general solution

\begin{equation}\hat T = C^ + \exp (u^ + z) + C^ - \exp (u^ - z),\end{equation}
T̂=C+exp(u+z)+Cexp(uz),
(A6)

where

\begin{equation}u_n ^ \pm = - \left( {\frac{{\lambda _2 }}{2}} \right) \pm \sqrt {\left( {\frac{{\lambda _2 }}{2}} \right)^2 + \lambda _1 }.\end{equation}
un±=λ22±λ222+λ1.
(A7)

Here n is an index for the layer number. To obtain the solution for an N layer system, we apply the multilayer formalism described by Feldman9 as follows. We can define a temperature vector based on equation (A6),

\begin{equation}\left[ {\begin{array}{*{20}c}{B_n^ + } \\[5pt]{B_n^ - }\end{array}} \right] \equiv\left[ {\begin{array}{*{20}c}{C_n^ + \exp (u_n^ + z)} \\[5pt]{C_n^ - \exp (u_n^ - z)} \\\end{array}} \right].\end{equation}
Bn+BnCn+exp(un+z)Cnexp(unz).
(A8)

Then to translate a distance hn from the bottom of layer n to the top, we can write

\begin{eqnarray}\left[ {\begin{array}{c}{C_n^ + \exp (u_n^ + z_t)} \\[5pt]{C_n^ - \exp (u_n^ - z_t)} \\\end{array}} \right] &=& \left[ {\begin{array}{c}{C_n^ + \exp (u_n^ + (z_b - h))} \\[5pt]{C_n^ - \exp (u_n^ - (z_b - h))} \\\end{array}} \right] \nonumber\\&=& \left[ {\begin{array}{c@{\quad}c}{\exp ( - u_n^ + h_n)} & 0 \\[5pt]0 & {\exp ( - u_n^ - h_n)} \\\end{array}} \right]\nonumber\\&&\times\left[ {\begin{array}{c}{C_n^ + \exp (u_n^ + z_b)} \\[5pt]{C_n^ - \exp (u_n^ - z_b)} \\\end{array}} \right],\end{eqnarray}
Cn+exp(un+zt)Cnexp(unzt)=Cn+exp(un+(zbh))Cnexp(un(zbh))=exp(un+hn)00exp(unhn)×Cn+exp(un+zb)Cnexp(unzb),
(A9)

or

\begin{equation}\left[ {\begin{array}{*{20}c}{B_n^ + } \\[5pt]{B_n^ - } \\\end{array}} \right]_t = \Pi (h_n)\left[ {\begin{array}{*{20}c}{B_n^ + } \\[5pt]{B_n^ - } \\\end{array}} \right]_b,\end{equation}
Bn+Bnt=Π(hn)Bn+Bnb,
(A10)

where

\begin{equation}\Pi (h_n) \equiv \left[ {\begin{array}{c@{\quad }c}{\exp ( - u_n^ + h_n)} & 0 \\[5pt]0 & {\exp ( - u_n^ - h_n)} \\\end{array}} \right],\end{equation}
Π(hn)exp(un+hn)00exp(unhn),
(A11)

which defines a transfer matrix for translation, Π(hn). We can obtain another transfer matrix for propagating the solution to the other side of an interface between materials. From layer n+1 to layer n, applying the condition that the temperatures and heat fluxes are equal on both sides of the interface, temperature continuity requires

\begin{equation}B_n^ + + B_n^ - = B_{n + 1}^ + + B_{n + 1}^ -.\end{equation}
Bn++Bn=Bn+1++Bn+1.
(A12)

The heat flux in the z direction in real-space coordinates is

\begin{equation}q_z = - k_{zz} \frac{{\partial T}}{{\partial z}} - k_{xz} \frac{{\partial T}}{{\partial x}} - k_{yz} \frac{{\partial T}}{{\partial y}}.\end{equation}
qz=kzzTzkxzTxkyzTy.
(A13)

Taking the Fourier transform with respect to the in-plane coordinates and the Laplace transform with respect to time

\begin{equation}\hat q_z = - k_{zz} \frac{{\partial \hat T}}{{\partial z}} - i(k_{xz} \xi + k_{yz} \eta)\hat T,\end{equation}
q̂z=kzzT̂zi(kxzξ+kyzη)T̂,
(A14)

or in Feldman notation,

\begin{eqnarray}\hat q_z &=& - [ {k_{zz} u^ + + i(k_{xz} \xi + k_{yz} \eta)}]B^ +\nonumber\\&& - [ {k_{zz} u^ - + i(k_{xz} \xi + k_{yz} \eta)}]B^ -.\end{eqnarray}
q̂z=[kzzu++i(kxzξ+kyzη)]B+[kzzu+i(kxzξ+kyzη)]B.
(A15)

It is convenient to define the quantity

\begin{equation}\gamma ^ \pm \equiv k_{zz} u^ \pm + i(k_{xz} \xi + k_{yz} \eta).\end{equation}
γ±kzzu±+i(kxzξ+kyzη).
(A16)

Then the heat flux matching condition is

\begin{equation}\gamma _n^ + B_n^ + + \gamma _n^ - B_n^ - = \gamma _{n + 1}^ + B_{n + 1}^ + + \gamma _{n + 1}^ - B_{n + 1}^ -.\end{equation}
γn+Bn++γnBn=γn+1+Bn+1++γn+1Bn+1.
(A17)

Equations (A12) and (A17) define a systems of equations

\begin{eqnarray}\left[ {\begin{array}{c}{\hat T} \\[5pt]{\hat q} \\\end{array}} \right] &=& \left[ {\begin{array}{c@{\quad}c}1 & 1 \\[5pt]{ - \gamma _n^ + } & { - \gamma _n^ - } \\\end{array}} \right]\left[ {\begin{array}{c}{B_n^ + } \\[5pt]{B_n^ - } \\\end{array}} \right]\nonumber\\&=& \left[ {\begin{array}{c@{\quad}c}1 & 1 \\[5pt]{ - \gamma _{n + 1}^ + } & { - \gamma _{n + 1}^ - } \\\end{array}} \right]\left[ {\begin{array}{c}{B_{n + 1}^ + } \\[5pt]{B_{n + 1}^ - } \\\end{array}} \right],\end{eqnarray}
T̂q̂=11γn+γnBn+Bn=11γn+1+γn+1Bn+1+Bn+1,
(A18)

which can be inverted to define a transfer matrix across the boundary,

\begin{equation}\left[ {\begin{array}{*{20}c}{B_n^ + } \\{B_n^ - } \\\end{array}} \right] = \Gamma _{n + 1 \to n} \left[ {\begin{array}{*{20}c}{B_{n + 1}^ + } \\{B_{n + 1}^ - } \\\end{array}} \right],\end{equation}
Bn+Bn=Γn+1nBn+1+Bn+1,
(A19)

where

\begin{equation}\Gamma _{n + 1 \to n} = \frac{1}{{\gamma _n^ - - \gamma _n^ + }}\left[\begin{array}{c@{\quad}c}{(\gamma _n^ - - \gamma _{n + 1}^ +)} & {(\gamma _n^ - - \gamma _{n + 1}^ -)} \\[6pt]{( - \gamma _n^ + + \gamma _{n + 1}^ +)} & {( - \gamma _n^ + + \gamma _{n + 1}^ -)} \\\end{array} \right].\end{equation}
Γn+1n=1γnγn+(γnγn+1+)(γnγn+1)(γn++γn+1+)(γn++γn+1).
(A20)

We can calculate the temperature and heat flux state at any point in the multilayer by using a series of transfer matrix operations if the state of the system is known at one point. The surface condition is calculated by iteration starting from the top of layer N (the bottom-most layer) to the top surface (n = 1) through successive application of Eqs. (A10) and (A19):

\begin{equation}\left[ {\begin{array}{c}{B_n^ + } \\[3pt]{B_n^ - } \\\end{array}} \right] = \Pi (h_n)\Gamma _{n + 1 \to n} \left[ {\begin{array}{c}{B_{n + 1}^ + } \\[3pt]{B_{n + 1}^ - } \\\end{array}} \right].\end{equation}
Bn+Bn=Π(hn)Γn+1nBn+1+Bn+1.
(A21)

If the bottom-most layer is semi-infinite, then in order for the temperature to be bounded at z → ∞,

\begin{equation}\left[ {\begin{array}{c}{B_N^ + } \\[3pt]{B_N^ - } \\\end{array}} \right] = \hat T_{S,N} (\xi,\eta,\omega)\left[ {\begin{array}{c}0 \\[3pt]1 \\\end{array}} \right].\end{equation}
BN+BN=T̂S,N(ξ,η,ω)01.
(A22)

The objective of the analysis is to calculate the surface temperature response of the system, |$\hat G = ( {\hat T_1 /\hat q_1 })_{surface} $|Ĝ=(T̂1/q̂1)surface which is then

\begin{equation}\hat G(f,\xi,\eta) = - \left( {\frac{{B_1^ + + B_1^ - }}{{\gamma _1^ + B_1^ + + \gamma _1^ - B_1^ - }}} \right).\end{equation}
Ĝ(f,ξ,η)=B1++B1γ1+B1++γ1B1.
(A23)

Note that the factor, |$\hat T_{S,N} (\xi,\eta,\omega)$|T̂S,N(ξ,η,ω), in Eq. (A22) always cancels as a result of taking the ratio in Eq. (A23) and so without loss of generality, it can be set to unity during evaluation.

1.
A. J.
Schmidt
,
X. Y.
Chen
, and
G.
Chen
,
Rev. Sci. Instrum.
79
(
11
),
114902
(
2008
).
2.
D. G.
Cahill
,
Rev. Sci. Instrum.
75
(
12
),
5119
(
2004
).
3.
K.
Kang
,
Y. K.
Koh
,
C.
Chiritescu
,
X.
Zheng
, and
D. G.
Cahill
,
Rev. Sci. Instrum.
79
(
11
),
114901
(
2008
).
4.
Y. K.
Koh
,
S. L.
Singer
,
W.
Kim
,
J. M. O.
Zide
,
H.
Lu
,
D. G.
Cahill
,
A.
Majumdar
, and
A. C.
Gossard
,
J. Appl. Phys.
105
(
5
),
054303
(
2009
).
5.
Y. K.
Koh
,
Y.
Cao
,
D. G.
Cahill
, and
D.
Jena
,
Adv. Funct. Mater.
19
(
4
),
610
(
2009
).
6.
R. M.
Costescu
,
M. A.
Wall
, and
D. G.
Cahill
,
Phys. Rev. B
67
(
5
),
054302
(
2003
).
7.
J. P.
Feser
and
D. G.
Cahill
,
Rev. Sci. Instrum.
83
(
10
),
104901
(
2012
).
8.
A. J.
Schmidt
,
R.
Cheaito
, and
M.
Chiesa
,
Rev. Sci. Instrum.
80
(
9
),
094901
(
2009
).
9.
A.
Feldman
,
High Temp.-High Pressure
31
(
3
),
293
(
1999
).
10.
J. M.
Powers
,
J. Heat Transfer
126
(
5
),
670
(
2004
).
11.
G.-M.
Choi
,
B.-C.
Min
,
K.-J.
Lee
, and
D. G.
Cahill
,
Nat. Commun.
5
,
4334
(
2014
).
12.
M. W.
Chase
,
NIST-JANAF Thermochemical Tables
, 4th ed. (
National Institute of Standards and Technology
,
Washington, DC/New York
,
1998
).
13.
G. R.
Stewart
,
Rev. Sci. Instrum.
54
(
1
),
1
(
1983
).
14.
Y. S.
Touloukian
,
Thermal Conductivity: Nonmetallic Solids
(
IFI/Plenum
,
New York
,
1970
), pp.
xxix
.
15.
Y. K.
Koh
and
D. G.
Cahill
,
Phys. Rev. B
76
(
7
),
075207
(
2007
).
16.
A. J.
Minnich
,
J. A.
Johnson
,
A. J.
Schmidt
,
K.
Esfarjani
,
M. S.
Dresselhaus
,
K. A.
Nelson
, and
G.
Chen
,
Phys. Rev. Lett.
107
(
9
),
095901
(
2011
).
17.
K. T.
Regner
,
D. P.
Sellan
,
Z. H.
Su
,
C. H.
Amon
,
A. J. H.
McGaughey
, and
J. A.
Malen
,
Nat. Commun.
4
,
1640
(
2013
).
18.
A. A.
Maznev
,
J. A.
Johnson
, and
K. A.
Nelson
,
Phys. Rev. B
84
(
19
),
195206
(
2011
).