Time-domain thermoreflectance and frequency-domain thermoreflectance (FDTR) have been widely used for non-contact measurement of anisotropic thermal conductivity of materials with high spatial resolution. However, the requirement of a high thermoreflectance coefficient restricts the choice of metal coating and laser wavelength. The accuracy of the measurement is often limited by the high sensitivity to the radii of the laser beams. We describe an alternative frequency-domain pump-probe technique based on probe beam deflection. The beam deflection is primarily caused by thermoelastic deformation of the sample surface, with a magnitude determined by the thermal expansion coefficient of the bulk material to measure. We derive an analytical solution to the coupled elasticity and heat diffusion equations for periodic heating of a multilayer sample with anisotropic elastic constants, thermal conductivity, and thermal expansion coefficients. In most cases, a simplified model can reliably describe the frequency dependence of the beam deflection signal without knowledge of the elastic constants and thermal expansion coefficients of the material. The magnitude of the probe beam deflection signal is larger than the maximum magnitude achievable by thermoreflectance detection of surface temperatures if the thermal expansion coefficient is greater than 5 × 10^{−6} K^{−1}. The uncertainty propagated from laser beam radii is smaller than that in FDTR when using a large beam offset. We find a nearly perfect matching of the measured signal and model prediction, and measure thermal conductivities within 6% of accepted values for materials spanning the range of polymers to gold, 0.1–300 W/(m K).

## I. INTRODUCTION

Pump-probe optical techniques based on time-domain thermoreflectance (TDTR)^{1–4} and frequency-domain thermoreflectance (FDTR)^{5} have been developed for non-contact measurement of thermal conductivity of thin film and bulk materials. TDTR and FDTR are widely used for the measurement of anisotropic thermal conductivity^{6–12} and thermal conductivity mapping with micrometer-scale spatial resolution.^{13–19} In these techniques, a metal coating serves as the thermometer by generating a thermoreflectance signal; i.e., the signal derives from the change of reflectivity *R* with temperature *T*. Combinations of the metal and laser wavelength that give large thermoreflectance coefficient $1RdRdT$ and thus the large magnitude of the signal are typically used to maximize the signal-to-noise ratio; common examples are Au at a laser wavelength of 532 nm and Al at a laser wavelength of 785 nm.^{20} However, for accurate measurements of in-plane thermal conductivity or high spatial resolution measurements of materials with low thermal conductivity, a metal with low thermal conductivity is needed to suppress lateral heat spreading in the metal coating.^{21,22}

The accuracy of TDTR and FDTR is often limited by the systematic error propagated from the uncertainties in the pump and probe beam radii.^{23} This problem is compounded when working at high spatial resolution because of the small depth of focus and the difficulty in accurately determining the intensity profile on small length scales. In a thermal conductivity mapping experiment, a drift in the position of the sample surface relative to the focal point of the laser beams propagates into systematic errors in thermal conductivity as a function of position.

In a TDTR or FDTR measurement scheme, the deflection (i.e., a small change in direction) of the probe beam that is offset with respect to the pump beam is an alternative signal for detecting the change of temperature. The beam deflection is created by thermoelastic deformation of the sample surface. Utilizing the probe beam deflection signal, Opsal *et al.* developed a technique for measuring the thickness of both transparent and opaque thin films with high spatial resolution.^{24}

We have previously used a frequency-domain probe beam deflection (FD-PBD) method to measure the thermal diffusivity of transparent polymers attached to an Al-coated silica substrate. With the laser beams passing through the polymer,^{25,26} the probe beam deflection signal is generated by the temperature dependence of the index of refraction *n*. This approach is limited to transparent materials and materials with small elastic constants so that the thermoelastic stresses in the material do not significantly deform the silica substrate.

In this paper, we describe a new FD-PBD method for measurement of thermal conductivity of bulk materials, which removes these limitations. Basically, the probe beam deflection signal is generated by the thermoelastic deformation of the bulk material to be measured. Metallic materials can be measured directly; other materials are first coated with a metal film. The signal is measured as a function of the modulation frequency of the pump beam, and the data are fitted to a model to extract the thermal conductivity of the material.

The rest of the paper is structured as follows: In Sec. II, the model is derived in detail. In Sec. III, the experimental measurement of probe beam deflection is explained. In Sec. IV, analysis of the magnitude of the signal, discussion of simplification of the model, sensitivity and uncertainty analysis, and fitting results of the experimental data are presented. We highlight several features of this method as compared with FDTR, including a larger magnitude of signal, flexible choices of the metal coating and laser wavelength, and smaller uncertainty propagated from the laser beam radii. We validate this work by demonstrating a nearly perfect matching of the model with the measured data and an overall agreement of measured thermal conductivities within 6% of accepted values in the range of 0.1–300 W/(m K).

## II. MODEL OF PROBE BEAM DEFLECTION

### A. Isotropic free thermal expansion model

This simplified model neglects all mechanisms of probe beam deflection except the deformation of the surface of the bulk material due to free thermal expansion, i.e., without mechanical constraint by the metal coating. The material is assumed to be a linearly elastic semi-infinite solid with isotropic elastic and thermal expansion properties. The heating is due to the absorption of a Gaussian pump beam by the metal coating. The solution of probe beam deflection for this situation has been reported in our previous work on probe beam deflection,^{27,28} which is based on a solution of surface displacement derived with the non-Fourier law of heat conduction equation that involves a relaxation time.^{29} Here, we obtain the temperature field with Fourier law of heat conduction equation instead, while the derivation of probe beam deflection given the temperature field essentially follows previous work^{27–29} and is elaborated here to fill in some intermediate steps in the derivation that are not included in the original paper.^{29}

Since heating by the pump beam has cylindrical symmetry, the temperature field *T* and displacement field ** u** also have cylindrical symmetry. The calculation can thus be greatly accelerated by selecting a cylindrical coordinate system and making use of Hankel transforms. As shown in Fig. 1(a), we choose

*r*= 0 at the center of the pump beam,

*z*= 0 at the surface of the bulk material, and $z\u0302$ pointing toward the interior of the sample. Via Helmholtz’s theorem, the displacement field $u=urr\u0302+uzz\u0302$ can be written as

where *φ* is a scalar dilatational potential, $\psi =\psi rr\u0302+\psi \theta \theta \u0302+\psi zz\u0302$ a vector potential, giving

where *ψ* is defined by $\psi \theta =\u22121r\u2202\psi \u2202r$. The governing equations of elasticity are given by^{29}

where $\gamma =\alpha T3\lambda +2\mu /\lambda +2\mu $, *λ* and *μ* the Lamé parameters, and *α*_{T} the linear coefficient of thermal expansion. With the mechanical constraint by the metal coating neglected, Eq. (2) is subjected to the stress-free condition at the interface between the metal coating and the bulk material,

Via Hooke’s law, the strain-displacement relation^{30} and Eq. (1), *σ*_{zr} and *σ*_{zz} in Eq. (3) can be expressed by *φ* and *ψ*, giving

After Hankel transform [$gr$ to $g\u0304k$] and Fourier transform [$ht$ to $h\u0303\omega $], the governing equations Eq. (2) become

while the boundary conditions Eq. (4) become

The temperature field $T\u0304\u0303$ is governed by the Fourier law of heat conduction after Hankel and Fourier transform,^{31}

where Λ_{z} and Λ_{r} are the cross-plane and in-plane thermal conductivity, respectively, and *C* = *ρc*_{p} the volumetric heat capacity, giving the temperature field in the bulk material,

where the surface temperature $Tbs\u2261T\u0304\u0303|z=0$ is solved by transfer matrix approach for heat conduction in semi-infinite air/metal coating/semi-infinite material system as detailed in the literature.^{2,31} Basically, Eq. (7) applies to each layer given the thermal conductivity and volumetric heat capacity. Adiabatic boundary condition applies at *z* → ±∞, while heat flux boundary condition applies at the air/metal interface with the heat flux equals

where *P*_{0} is the amplitude of the fundamental harmonic of the absorbed power of modulated pump laser, *w*_{0} the 1/e^{2} radii of the pump beam.

With the temperature field given by Eq. (8), the general solution of the governing equations Eq. (5) is

where the constants *a*_{φ} and *a*_{ψ} are uniquely determined by the boundary conditions Eq. (6) to be

The surface displacement of this “isotropic free thermal expansion” model, $Z\u0304\u0303iso\u2212freek,\omega \u2261\u2212u\u0304\u0303z|z=0$ (the minus sign makes displacement outward positive) is given by

The “low frequency approximation” is applicable is most cases,

To see this, given a radii of pump beam *w*_{0} ∼ 10 *µ*m, consider a relatively extreme case of frequency as high as *f* = *ω*/2*π* ∼ 1 MHz and speed of sound as low as *v*_{T}, *v*_{L} ∼ 1000 m/s, then $\omega 2/vT2,\omega 2/vL2$ ∼ (2*π*/1000 *µ*m)^{2} ≪ (2*π*/*w*_{0})^{2} ∼ *k*^{2}. Note that the frequency of interest is typically smaller than 1 MHz because the frequency of interest shifts to smaller values with decreasing thermal diffusivity (as shown later) and 1 MHz is for a high thermal diffusivity material such as gold.

where the second equality is via the conversion formulas between elastic properties and *ν* is the Poisson’s ratio. In this low frequency approximation, the frequency dependence of the surface deformation on elastic properties is eliminated, which is desirable for frequency-domain measurement of thermal conductivity, heat capacity, and coefficient of thermal expansion.

An inverse Hankel transform is then applied to obtain the surface displacement in real space,

The surface slope is thus

If the surface slope is nearly constant in the range of *r* covered by the probe beam, the probe beam deflection angle would simply be twice the surface slope. This condition is not well satisfied for a typical experiment because the radii of the probe beam are comparable to the radii of the pump beam. Nevertheless, we previously showed that the following convolution of the probe intensity with the surface slope gives a good description of the probe beam deflection angle,^{32}

where *C*_{probe} is a material-independent constant on the order of unity, *w*_{0} the 1/e^{2} radii of the pump and probe beam, *r*_{0} the offset distance between the pump and probe beam.

Note that it is assumed that the probe beam deflection at modulation frequency *f* is not influenced by the steady-state surface deformation due to heating by the probe beam and dc component of the power of the pump beam, i.e., the linear superposition principle applies to surface deformations. This is valid if the product of the coefficient of thermal expansion and temperature excursion is small enough so that the radius of curvature of the sample surface is large compared with the focal length of the objective lens, which is typically true in our method.

### B. The full model

The deformation of the surface of the sample (metal coating/bulk material) can be different from the prediction of the isotropic free thermal expansion model due to two issues. The first issue is the effects of the metal coating on the deformation of the sample surface, including the thermal expansion of the metal coating, the elastic deformation of the bulk material created by the thermal expansion stress of the metal coating, and the constraint of the thermal expansion strain of the bulk material by the stiffness of the metal coating. The second issue is the anisotropy of elastic constants and thermal expansion coefficients of the metal coating and the bulk material. Furthermore, the isotropic free thermal expansion of the model neglects probe beam deflection created by mechanisms other than surface deformation.

To address these drawbacks of the isotropic free thermal expansion model, we present here “the full model” of probe beam deflection for metal-coated bulk material measured in air. The sample surface deformation is calculated for a finite thickness metal coating on a semi-infinite bulk material, where both the metal coating and the bulk material have the anisotropic thermal expansion, thermal conductivity, and elastic constants; i.e., thermal expansion and thermal conductivity are described by second rank tensors and the elastic constants are described by fourth rank tensors. We also include the probe beam deflection contributed by the changes of optical path length in air and changes in the phase of the reflection coefficient due to the temperature dependence of the refractive index of air and metal coating.

The Cartesian coordinate system is set up with *x* = *y* = 0 at the center of the pump beam, *z* = 0 at the surface of the bulk material, and $z\u0302$ pointing toward the interior of the sample. The governing equations of elasticity as expressed by the displacement $u=uxx\u0302+uyy\u0302+uzz\u0302$ are given by^{30}

where the stresses are written in terms of displacements by

via the inverted Duhamel-Neumann law and strain-displacement relation,^{30} where *α*_{x}, *α*_{y}, and *α*_{z} are the linear coefficient of thermal expansion along $x\u0302$, $y\u0302$, and $z\u0302$, respectively (assuming all off-diagonal components in the thermal expansion tensor being zero), the elastic constants $Cmn\u2032$ are a component of the forth-rank tensor of stiffness under the current coordinate system (the single number *m*, *n* = 1–6 is the abbreviated subscripts of the pairs 11, 22, 33, 23, 31, and 12). Since elastic constants are typically known under the coordinate system of crystal axes, the transformation of the tensor of stiffness needs to be performed considering the orientation of the crystal axes with respect to the current coordinate system^{33}.

The governing equations Eq. (18) are subjected to stress-free conditions at the surface of the metal coating and continuity conditions at the interface between the metal coating and the bulk material as follows:^{30}

where *L* is the thickness of the metal coating.

With the approximation of the heating by the laser entering from the sample surface, and assuming all off-diagonal components in the thermal conductivity tensor being zero, the temperature field *T* is governed by

After Fourier transforms [$hx,y,t$ to $h\u0302\eta ,\xi ,\omega $], the governing equations Eq. (18) and third, fourth, and fifth rows of Eq. (19) are rearranged to the following form,^{34}

where the matrices are given by

The boundary conditions Eq. (20) after Fourier transforms are given by

The governing equation of temperature field Eq. (21) after Fourier transforms is given by

which yields the following solution

where the subscript of Λ_{x}, Λ_{z}, *C*, and *ζ* being 1, 2, and 3 indicating air, the metal coating, and the bulk material, respectively, *G* the thermal boundary conductance of the interface of the metal coating and the bulk material, $T\u0302bs\u2261T\u0302|z=0$ and $T\u0302s\u2261T\u0302|z=\u2212L$ are solved by transfer matrix approach for heat conduction in semi-infinite air/metal coating/semi-infinite material system as detailed in the literature.^{2,30}

With the temperature field given by Eq. (28), the general solution of the governing equation Eq. (22) is given by

where *m* = 1–6, the vector ** J** is composed of constants to be determined by boundary condition; each column of the six-by-six matrix

**is an eigenvector of**

*Q*

*A*^{−1}

**, with the corresponding component of the vector**

*B***being the eigenvalue to which the eigenvector belongs; the vector**

*λ***=**

*U*

*Q*^{−1}

**. The subscript of**

*AD***,**

*Q***,**

*λ***, and**

*U***being 2 and 3 indicating the metal coating and the bulk material, respectively. Note that**

*J*

*Q*_{3}and

*λ*_{3}are arranged such that the real part of the first three components of

*λ*_{3}are negative, representing wave propagation along the positive

*z*direction. To make the determination of the sign of components of

*λ*_{3}free from influence of numerical error, we adopt the strategy of including a small imaginary part in the elastic constants, e.g., by multiplying (1 +

*i*× 10

^{−4}).

To find the sample surface displacement $Z\u0302ani\eta ,\xi ,\omega \u2261\u2212u\u0302z|z=\u2212L=$ −** S**(3)|

_{z=−L}, the constants

*J*_{2}are determined by solving the linear equations resulting from substitution of the general solution Eq. (29) for

**in the boundary conditions Eq. (26). The sample surface displacement is given by**

*S*The probe beam deflection due to changes in optical path length with temperature can be regarded as resulting from an effective surface displacement.^{27} The effective surface displacement due to changes of optical path length in the air with temperature is given by^{27}

where $dnairdT$ is the temperature derivative of the refractive index of air *n*_{air}. In addition, at the air/metal coating interface, changes in the phase of the Fresnel reflection coefficient *ϕ*_{r} with temperature give the effective surface displacement,^{27}

where *λ*_{1} is the wavelength of the probe beam, $d\varphi rdT$ the temperature derivative of the phase of the Fresnel reflection coefficient, a constant determined by

where *n*_{metal} and *κ*_{metal} are the real part and imaginary part of the refractive index of the metal coating, respectively. Equation (33) assumes that the metal coating is optically semi-infinite for the probe beam.

The probe beam deflection contributed due to each surface displacement $Z\u0302i\eta ,\xi ,\omega $ (where *i* refers to *ani*, *air*, or *r*) is then calculated in a similar way as described for the isotropic free expansion model. Inverse Fourier transforms are applied to obtain the surface displacement in real space,

By making the transforms of variables

the surface slope is found to be

Following the approximation of convolution of the probe intensity with the surface slope, the probe beam deflection angle is given by

where *C*_{probe} is a material-independent constant on the order of unity, *w*_{0} the 1/e^{2} radii of the probe beam, *r*_{0} the offset distance between the pump and probe beam, and *φ*_{0} the angle of the vector from the center of the pump beam to that of the probe beam in the current coordinate system. The MATLAB codes for the implementation of the above models are available for download at: https://zenodo.org/badge/latestdoi/529379931.

## III. MEASUREMENT OF PROBE BEAM DEFLECTION

The optics of the frequency-domain probe-beam-deflection (FD-PBD) system are shown in Fig. 1(b). The only significant difference between our FD-PBD system and a typical FDTR measurement system is our use of a quadrant-cell photoreceiver for detecting modulations in beam displacement in place of the single photodiode or balanced photoreceiver for detecting modulation of beam intensity. Two free-space diode lasers at a wavelength of 660 and 640 nm are used as the pump beam and probe beam, respectively. The two laser beams are focused by the objective lens and incident on the Al coated bulk material at normal incidence with a beam offset along the vertical direction. With the pump beam modulated at frequency *f*, the heating due to absorption by the Al coating creates elastic deformation of the sample surface and consequently an angular deflection of the probe beam along the vertical direction synchronous with the modulation frequency. After the reflected probe beam passes back through the objective lens, this angular deflection is converted to a transverse displacement along the vertical direction, which is measured by the difference signal (between the upper cells and lower cells) of the quadrant-cell photoreceiver, which has been aligned with the probe beam. With this signal as the input and the modulation signal of the pump beam as the reference, a lock-in amplifier yields the in-phase signal *X* and out-of-phase signal *Y* at frequency *f*.

Due to the finite bandwidth of the quadrant-cell photoreceiver, the frequency dependence of the measured signal contains significant contributions from the quadrant-cell photoreceiver at high frequencies, >1 kHz for the Newport^{™} 2901 quadrant-cell photoreceiver. To extract the true PBD signal, we divide the complex signal *X* + *iY* at each *f* by the response of the quadrant-cell photoreceiver measured at the same *f*. We obtain the response of the photoreceiver as follows: First, we block the probe laser and remove short-pass optical filter to detect the modulated pump laser using the quadrant-cell photoreceiver. Second, we collect the complex signal from the lock-in amplifier with the signal from quadrant-cell photoreceiver as the input and the modulation signal of the pump beam as the reference. Finally, the frequency response of the photoreceiver is obtained by normalizing this complex signal by its value in the low frequency limit.

To extend the measurements to frequencies *f* above the upper limit (100 kHz) of the quadrant-cell photoreceiver, we implement heterodyne detection and modulate the pump at frequency *f*, modulate the probe beam at frequency *f* + *f*_{dec}, and detect the signal at the difference frequency *f*_{dec}. We typically set *f*_{dec} = 5 kHz. The reference signal for the lock-in amplifier is generated by inputting the two modulation signals into a frequency mixer (Mini-Circuit® ZAD-6+) and extracting the difference signal at frequency *f*_{dec} with a lowpass filter (Thorlabs EF114).

The procedure of the experiment is as follows: First, we focus the pump and probe beams on the surface of the sample by moving the sample on the translation stage and focusing the optical system on defects or dusts on the sample surface using the video camera. Second, we align the quadrant-cell photoreceiver with the probe beam by adjusting the transverse position of the photoreceiver until the sum signal of all cells is a maximum and the difference signal between the upper cells and lower cells is zero. Third, we align the pump beam with the probe beam by rotating the beam-splitter near the objective lens until the two beam spots overlap as seen from the video camera and the lock-in amplifier yields zero signal. Fourth, we set the beam offset by rotating the beam-splitter using a differential micrometer. The proportionality factor between the beam offset distance and the position of the micrometer has been measured before the experiment. For a 10× microscope objective (focal length of 20 mm) and the Newport SL8A gimbal mount, beam offset is 4.9 per 10 *µ*m movement of the micrometer in the vertical direction. Finally, we collect the in-phase and out-of-phase signals at a series of modulation frequencies using the lock-in amplifier.

## IV. RESULTS AND DISCUSSION

### A. Analysis of calculated PBD signal

As shown in Fig. 2, the probe beam deflection angles (PBD) are calculated using the full model [Eq. (37) with *C*_{probe} = 1] for four representative samples (polystyrene, SrTiO_{3}, MoS_{2}, and Au) measured in air. The real part and imaginary part of the calculated complex number $\theta \u0303i$ are the in-phase and out-of-phase PBD (in unit of *μ*rad), respectively. Three sources of PBD are present, i.e., surface deformation (in red), changes of optical path length in air (in blue), and changes of the phase of the reflection coefficient (in green). Figure 2(a) shows PBD for polystyrene, an amorphous polymer with a low thermal conductivity of 0.156 W/(m K). Due to the high coefficient of thermal expansion *α*_{T} = 73 × 10^{−6} K^{−1}, the PBD due to surface deformation completely dominates the total PBD. Figure 2(b) shows PBD for cubic SrTiO_{3} (100), an oxide crystal with a thermal conductivity of 11.0 W/(m K). Due to the in-plane anisotropy of elastic constants, the PBD due to surface deformation is dependent on the in-plane orientation of beam offset. It is assumed here that the beam offset is along [010]. With *α*_{T} = 10.4 × 10^{−6} K^{−1}, the PBD due to surface deformation still dominates the total PBD. Figure 2(c) is for bulk hexagonal MoS_{2}, a material with highly anisotropic thermal expansion, thermal conductivity, and elastic constants. Due to the low coefficient of thermal expansion along the radial direction (5.8 × 10^{−6} K^{−1}), PBD due to surface deformation is no longer large enough to dominate the total PBD, with the PBD due to air also playing an important role. Figure 2(d) shows PBD for bare Au, a cubic crystal metal with a high thermal conductivity of 315 W/(m K). As exemplified by Au, metals can be measured without metal coating, since they are reflective and absorptive with the optical absorption depth much smaller than the laser beam radii. Calculation of PBD for bare metal can be done using the same formalism as that described above with the metal coating set to the same material as the bulk sample.

To gain insight on the dependence of the PBD signal on materials parameters, the frequency and value of the peak of the out-of-phase PBD are listed in Table I, together with the thermal diffusivity and coefficient of thermal expansion. The peak frequency of the PBD due to surface deformation scales, to a good approximation, with the in-plane thermal diffusivity of the bulk material. This makes it possible to measure the thermal conductivity without knowing the experimental parameters that only influence the magnitude of the signal, for example, the incident laser power and coefficient of optical absorption. In addition, as expected, the peak value of PBD due to air and reflection coefficient does not vary significantly with the type of material under investigation, while PBD due to surface deformation scales with the coefficient of thermal expansion. This observation motivates the following estimate of the order of magnitude of the PBD angle due to surface deformation,

where Δ*T* is ac temperature excursion (i.e., the amplitude of surface temperature rise at modulation frequency).

. | Polystyrene . | SrTiO_{3}
. | MoS_{2}
. | Au . | ||||
---|---|---|---|---|---|---|---|---|

f (kHz)
. | PBD (μrad)
. | f (kHz)
. | PBD (μrad)
. | f (kHz)
. | PBD (μrad)
. | f (kHz)
. | PBD (μrad)
. | |

Surface deformation | 0.23 | 56.79 | 6.7 | 7.57 | 60 | 2.69 | 212 | 13.02 |

Reflection coefficient | 2.36 | −0.09 | 10.8 | −0.12 | 108 | −0.12 | 302 | −0.01 |

Air | 1.63 | 0.41 | 9.6 | 0.53 | 37 | 0.56 | 51 | 0.54 |

Total | 0.23 | 57.06 | 6.7 | 7.97 | 53 | 3.12 | 212 | 13.34 |

D = Λ/C (×10^{−6} m^{2}/s) | 0.12^{a} | 4.0^{b} | In/through-plane 42/2.5^{c} | 127^{d} | ||||

α_{T} (10^{−6}/K) | 73^{e} | 10.4^{f} | In/through-plane 5.8/11.3^{g} | 14.2 |

. | Polystyrene . | SrTiO_{3}
. | MoS_{2}
. | Au . | ||||
---|---|---|---|---|---|---|---|---|

f (kHz)
. | PBD (μrad)
. | f (kHz)
. | PBD (μrad)
. | f (kHz)
. | PBD (μrad)
. | f (kHz)
. | PBD (μrad)
. | |

Surface deformation | 0.23 | 56.79 | 6.7 | 7.57 | 60 | 2.69 | 212 | 13.02 |

Reflection coefficient | 2.36 | −0.09 | 10.8 | −0.12 | 108 | −0.12 | 302 | −0.01 |

Air | 1.63 | 0.41 | 9.6 | 0.53 | 37 | 0.56 | 51 | 0.54 |

Total | 0.23 | 57.06 | 6.7 | 7.97 | 53 | 3.12 | 212 | 13.34 |

D = Λ/C (×10^{−6} m^{2}/s) | 0.12^{a} | 4.0^{b} | In/through-plane 42/2.5^{c} | 127^{d} | ||||

α_{T} (10^{−6}/K) | 73^{e} | 10.4^{f} | In/through-plane 5.8/11.3^{g} | 14.2 |

To check Eq. (38), we fix Δ*T* = 3 K in the low frequency limit (for each sample, the heating power has been chosen such that the amplitude of surface temperature oscillation in the low frequency limit is equal to 3 K), and therefore, Eq. (38) gives an estimate of 220, 31, 23, and 43 *µ*rad for polystyrene, SrTiO_{3}, MoS_{2} (with $\alpha T=23\alpha T,in+13\alpha T,through$), and Au, respectively. These estimates compare well in terms of order of magnitude with the calculations of the full model of 174, 18, 6, and 31 *µ*rad shown in Fig. 2 except for the highly anisotropic MoS_{2}. Note that the isotropic free thermal expansion model gives 177, 20, 5, and 27 *µ*rad, respectively, nearly identical to the results of the full model.

To compare the magnitude of the signals measured by the FD-PBD system with that of FDTR, we convert the probe beam deflection *θ* to the PBD lock-in-amplifier signal normalized by the sum signal of the photoreceiver. After the objective lens, the transverse displacement of the unfocused probe beam is

where *f*_{0} is the focal length of the objective. The intensity profile of the unfocused Gaussian probe beam is

where *P*_{1} is the total power, *I*_{1} the peak intensity, *w* and *w*_{0} the 1/e^{2} radii of the unfocused and focused probe beam, respectively. With *d* ≪ *w*, the ratio of interest, i.e., the difference of power (between the upper cells and lower cells) over the total power is given by

with *λ*_{1} = 640 nm and *w*_{0} = 8.3 *µ*m.

Using Eq. (38), this ratio, i.e., the normalized PBD signal, is estimated to be 65*α*_{T}Δ*T*. For most materials, *α*_{T} is larger than 5 × 10^{−6} K^{−1}, which gives a normalized PBD signal >3 × 10^{−4} K^{−1}Δ*T*. For comparison, the corresponding value in an FDTR measurement is given by $1RdRdT\Delta T$. Among the common combination of metal coating and wavelength of the laser, Au at the wavelength of 532 nm gives the maximum value of $1RdRdT\Delta T$ = 3 × 10^{−4} K^{−1}Δ*T*.^{20} Thus, the magnitude of the PBD signal is generally larger than that of FDTR, and more importantly, independent of the choice of metal coating and laser wavelength.

### B. Applicability of the isotropic free thermal expansion model

For extracting thermal conductivity by fitting to the measured PBD data, the isotropic free thermal expansion model is preferred compared with the full model, as the elastic constants and thermal expansion coefficient of the bulk material are not needed. This can be seen from Eqs. (14) and (17), where the Poisson ratio and coefficient of thermal expansion coefficient only enter the results as prefactors. To have a sense of the applicability of the isotropic free thermal expansion model, its error is evaluated by fitting it to the PBD data generated by the full model and then monitoring the error of the fitted thermal conductivity (Fig. 3). Due to the symmetry of the crystal structure, the PBD and consequently the fitted thermal conductivity of c-plane sapphire, MoS_{2} (001), and Bi_{2}Se_{3} (001) are independent of the in-plane orientation of beam offset. This is not the case for cubic crystals, for which a bar is used to indicate the range of fitted thermal conductivity with a circle indicating the average of these fitted values.

In this analysis, the materials are assumed to be measured in a vacuum to avoid the contribution coming from the PBD signal generated by the temperature field in the air above the sample. Otherwise, the fact that the simplified model neglects the PBD contribution from the air would have caused a noticeable error for materials with relatively small *α*_{T} such as borosilicate glass and MoS_{2}. Furthermore, we set the Al coating thickness to 32 nm to reduce the effects of the metal coating on the deformation of the sample surface while ensuring that less than 1% of the pump laser power enters the bulk material.

The error inherent in the simplified model is less than ∼6% for this selection of materials and measurement conditions. A wide range of applicability of the simplified model is also suggested by the calculation of its the errors as a function of material properties. Figure S1 shows that the error due to neglecting the PBD contribution from the air is within 2% for most polymeric materials (*α*_{T} > 50 × 10^{−6} K^{−1}),^{35} and within 6% for most non-polymeric materials (*D* > 0.5 mm^{2}/s) if *α*_{T} > 6 × 10^{−6} K^{−1}. Note that this error can be readily removed by performing experiments with the sample in a vacuum chamber. Figure S2 shows that the error due to neglecting effects of the Al coating on surface deformation is within 7% for most polymeric materials (*α*_{T} > 50 × 10^{−6} K^{−1}) if *E* > 100 MPa, and within 6% for most non-polymeric materials (*E* > 10 GPa) if *α*_{T} > 5 × 10^{−6} K^{−1}. A further consequence of this error analysis is that when using the full model, uncertainties in the elastic constants and thermal expansion coefficient have only a small effect on the accuracy of the fitted thermal conductivity.

### C. Sensitivity and uncertainty

To quantify the sensitivity of the PBD angle to the parameters in the model, we calculate the sensitivity of in-phase PBD angle, $Re(\theta \u0303iso\u2212free)$, and out-of-phase PBD angle, $Im(\theta \u0303iso\u2212free)$, to a parameter *x* in the isotropic free thermal expansion model [Eq. (17)] as defined by

where *In*_{max} and *Out*_{max} are the in-phase PBD angle and out-of-phase PBD angle with maximum magnitude in the frequency range studied. For an Al coated SrTiO_{3} sample, as shown in Fig. 4, *S*_{in} and *S*_{out} are calculated for parameter *x* as follows: in-plane thermal conductivity Λ_{r}, through-plane thermal conductivity Λ_{z}, volumetric heat capacity *C* of SrTiO_{3}; thickness *L*_{Al}, thermal conductivity Λ_{Al}, and volumetric heat capacity *C*_{Al} of Al; thermal boundary conductance of Al/SrTiO_{3} interface *G*; pump and probe beam radii *w*_{0}; and the beam offset *r*_{0}. By comparison of the first row [Figs. 4(a) and 4(b)] with *w*_{0} = 8.3 *µ*m, *r*_{0} = 19.6 *µ*m and the second row [Figs. 4(c) and 4(d)] with *w*_{0} = 8.3 *µ*m, *r*_{0} = 9.8 *µ*m, it is seen that with higher ratio of *r*_{0}/*w*_{0}, sensitivity to *w*_{0} is lower while sensitivity to *r*_{0} is higher. Since beam offset *r*_{0} can be measured at a higher accuracy than the beam radii *w*_{0}, *w*_{0} = 8.3 *µ*m, *r*_{0} = 19.6 *µ*m should reduce the total uncertainty due to the uncertainty of parameters. For comparison, the sensitivity of the FDTR signal of the same material under typical measurement conditions is shown in Fig. S3. Comparison with Figs. 4(a) and 4(b) shows that with a relatively large ratio of beam offset and beam radii, FD-PBD significantly improves the sensitivity to thermal conductivity as compared with sensitivity to beam radii.

Note that for SrTiO_{3}, Λ_{r} = Λ_{z}, but the sensitivity to Λ_{r} is significantly larger than the sensitivity to Λ_{z} at most frequencies. The reason is that PBD is primarily dependent on the radial derivative of temperature and is thus primarily sensitivity to the in-plane thermal diffusion. Partly due to the same reason, sensitivity to the thermal boundary conductance *G* is nearly zero. The reason for the low sensitivity to *G* is the large thermal penetration depth ∼ *d*_{p} = $D/\pi f$ (14 *µ*m at out-of-phase peak frequency *f* = 6.7 kHz).

To quantify the uncertainty of thermal conductivity due to the parameters, uncertainty due to each parameter *x* is calculated as follows: The PBD data are first calculated using the isotropic free thermal expansion model with the value of *x* increased by its uncertainty. The thermal conductivity is then extracted from this PBD data by fitting using the same model with the normal value of *x*. The uncertainty of thermal conductivity due to the uncertainty of *x* is evaluated as the error of the fitted thermal conductivity. With the values and uncertainties of parameters in Table II, the calculated uncertainties of thermal conductivity are shown in Fig. 5. A large uncertainty is assumed for thermal boundary conductance *G*, considering its significant variation with surface preparation.

C
. | Λ_{Al}
. | C_{Al}
. | L_{Al}
. | w_{0}
. | r_{0}
. | G
. |
---|---|---|---|---|---|---|

2 | 165 | 2.42 | 32 | 8.3 | 19.6 | 40 |

MJ/(m^{3} K) | W/(m K) | MJ/(m^{3} K) | nm | μm | μm | MW/(m^{2} K) |

3% | 5% | 3% | 5% | 8% | 2% | 50% |

C
. | Λ_{Al}
. | C_{Al}
. | L_{Al}
. | w_{0}
. | r_{0}
. | G
. |
---|---|---|---|---|---|---|

2 | 165 | 2.42 | 32 | 8.3 | 19.6 | 40 |

MJ/(m^{3} K) | W/(m K) | MJ/(m^{3} K) | nm | μm | μm | MW/(m^{2} K) |

3% | 5% | 3% | 5% | 8% | 2% | 50% |

Three different ways of fitting are considered, i.e., fitting to both in-phase and out-of-phase PBD (lines), only in-phase PBD (open circles), and only out-of-phase PBD (filled circles). The total uncertainty (in black) is the square root of the sum of the square of uncertainties of all parameters. Fitting to both in-phase and out-of-phase PBD gives a total uncertainty roughly constant at 5% when thermal conductivity varies within 0.1–300 W/(m K). The uncertainty of the beam offset *r*_{0} is identified as the dominant source of the total uncertainty, due to the highest sensitivity to *r*_{0} (when using *w*_{0} = 8.3 *µ*m, *r*_{0} = 19.6 *µ*m). If using *w*_{0} = 8.3 *µ*m, *r*_{0} = 9.8 *µ*m instead, the uncertainty of *w*_{0} will become the dominant source of total uncertainty and result in a larger total uncertainty due to the higher uncertainty of *w*_{0} than *r*_{0} (see Table II). We highlight the feature of FD-PBD in terms of reducing the uncertainty due to beam radii by comparing with FDTR, as shown in Fig. S4. Comparison with the solid blue curve in Fig. 5 shows that with a relatively large ratio of beam offset and beam radii (*w*_{0} = 8.3 *µ*m, *r*_{0} = 19.6 *µ*m), FD-PBD reduces the uncertainty due to beam radii by an order of magnitude.

### D. Fitting to measured PBD signal

To provide a real-world example, we plot as Fig. 6(a) the best fit of the measured PBD signal of Al coated SrTiO_{3} to the isotropic free thermal expansion model with the thermal conductivity as the only fitting parameter (during the fitting, in-phase and out-of-phase PBD signal are first divided by their respective maximum magnitude in the frequency range). In total, 50 data points are collected over a frequency range of a factor of 100, from a frequency ten times smaller than the peak in the out-of-phase signal to the frequency ten times larger than the peak. The parameters used in the fitting include Λ_{air} = 0.028 W/(m K), *C*_{air} = 1192 J/(m^{3} K) of air; *L*_{Al} = 85 nm, Λ_{Al} = 165 W/(m K), *C*_{Al} = 2.42 MJ/(m^{3} K); pump and probe beam 1/e^{2} radii *w*_{0} = 8.3 *µ*m, beam offset *r*_{0} = 9.8 *µ*m; *G* = 40 MW/(m^{2} K), *C* = 2.74 MJ/(m^{3} K) for SrTiO_{3}. Nearly perfect agreement between the fitting curve and the measured signal is found, yielding the thermal conductivity of SrTiO_{3} Λ = 10.4 W/(m K), only 5% smaller than the accepted value of 11.0 W/(m K).

The normalized PBD signal given by Eq. (41) is equal to the ratio of the amplitude of the PBD signal, *V*, over the sum signal of all cells of the quadrant-cell photoreceiver, SUM, thus giving the relation,

The experimental data in Fig. 6(a) are measured with SUM = 0.842 V. Therefore, a typical magnitude of PBD signal ∼100 *µ*V in Fig. 6(a) corresponds to a PBD angle *θ* ∼ 2 *µ*rad. This signal is accompanied by a steady state temperature rise due to pump laser ∼1.0 K and ac temperature excursion of Δ*T* ∼ 0.8 K at the lower end of the frequency range (as calculated using the measured absorbed power of pump laser *P*_{0} = 0.36 mW). With the time constant of the lock-in amplifier set to 1 s, data collection takes 5 s at each frequency and 250 s in total, while at *f* > 4 kHz, the noise of in-phase and out-of-phase signal (as measured by the standard deviation) is ∼0.1 *µ*V, i.e., ∼2 nrad according to Eq. (43). This noise is equal to that measured at the limit of the power of probe beam approaching zero, and is supposed to originate from the electronic noise of the quadrant-cell photoreceiver and the lock-in amplifier. At *f* < 1 kHz, the noise scales as 1/*f* and increases with the power of the probe beam. We conclude that the PBD signal has an excellent signal-to-noise ratio and is suitable for high throughput measurement.

In Fig. 6(b), we compare thermal conductivities measured by FD-PBD with accepted values for selected materials that span a factor of 3000 in thermal conductivity Measured thermal conductivities agree well with the accepted values except for borosilicate glass. Therefore, the full model is used to fit the same set of measured PBD data for borosilicate glass, and the accuracy of the measurement is greatly improved as shown by the triangle for borosilicate. We attribute the relatively large error of the isotropic free thermal expansion model for the borosilicate glass sample to the combination of the small thermal expansion coefficient, the contribution of air to the PBD signal, and the significant effects of the relatively thick Al film (87 nm) on the surface deformation. We also checked the full model for fitting the data for polystyrene. The full model gives essentially the same value of thermal conductivity as the isotropic free thermal expansion model. As an indicator of the overall accuracy of the measurement, the root means square of the relative deviation of the measured values from the accepted values is 5.7%.

## V. CONCLUSION

In summary, we developed a frequency-domain probe beam deflection (FD-PBD) method for the measurement of the thermal conductivity of bulk materials. We derived a simplified model for probe beam deflection, the isotropic free thermal expansion model, which allows measurement of thermal conductivity without knowing the elastic and thermal expansion properties of the material. This simplified model is shown to be widely applicable with small errors for measurement with relatively thin metal coating in a vacuum. We also derived the full model that incorporates the anisotropic elastic constants and thermal expansion coefficient of the material that is under investigation. We established two unique features of this method as compared with FDTR: (1) a larger signal independent of the type of metal coating and laser wavelength; (2) reduced uncertainty propagated from pump and probe beam radii if using a relatively large beam offset. Finally, we validated this method by matching the measured signal and model prediction and evaluating the agreement of measured thermal conductivities with accepted values in the range of 0.1–300 W/(m K).

## SUPPLEMENTARY MATERIAL

See the supplementary material for calculations and discussion are included concerning: error of isotropic free thermal expansion model as a function of materials properties; sensitivity and uncertainty analysis of FDTR; comparison with previous FD-PBD method in terms of sensitivity; and simultaneous extraction of two thermal properties by combining FDTR.

## ACKNOWLEDGMENTS

This research was supported by the NSF through the University of Illinois at Urbana-Champaign Materials Research Science and Engineering Center Grant No. DMR-1720633.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Jinchi Sun**: Conceptualization (supporting); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Guangxin Lv**: Data curation (supporting); Investigation (supporting); Methodology (supporting); Resources (supporting); Software (supporting); Supervision (supporting); Writing – review & editing (supporting). **David G. Cahill**: Conceptualization (lead); Data curation (supporting); Formal analysis (supporting); Funding acquisition (lead); Investigation (supporting); Methodology (equal); Project administration (lead); Resources (lead); Software (supporting); Supervision (lead); Validation (supporting); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.