Anisotropic thermoreflectance thermometry: A contactless frequency-domain thermoreflectance approach to study anisotropic thermal transport

The study of thermal anisotropy in solids has recently drawn considerable scientific and technological attention in different fields.^1–9 The ability to control heat transport through intrinsic or purposely engineered thermal anisotropy has the potential to open new routes to develop novel concepts toward smart heat manipulation.^10–13 In thermally anisotropic materials, thermal anisotropy is evidenced by the different elements of the thermal conductivity tensor (κ_ij), leading to the tensorial expression of Fourier’s law, q_i = κ_ij(∂T/∂x^j), where q_i represents each of the components of the vectorial heat flux, T is the temperature, and x^j is the spatial coordinate. The development of novel experimental methodologies to study anisotropic thermal transport has recently become a relevant research objective. A considerable number of experimental techniques and methodologies^9,14–27 based on the variations of the 3-omega method,^28,29 time-domain thermoreflectance,³⁰ and frequency-domain thermoreflectance³³ have been developed for this purpose, demonstrating their capability to obtain the components of κ_ij. The main differences between these approaches are the dimensionality of the heat source (line or spot) and their contact or contactless fashion (electrical resistor or focused optical spot). The techniques developed in Refs. 15 and 16 are based on the geometry used in the 3-omega method, i.e., a long and narrow line-shaped heat source. This geometry provides a main advantage, simultaneous thermal sensitivity to only two crystallographic directions, which arises from the one-dimensional (1D) geometry of the heat source. In other words, the temperature distribution in the direction along the heat source (long narrow line) is uniform; hence, no heat flow occurs in this direction, considerably simplifying the data analysis process. The main drawbacks of the approaches given in Refs. 15 and 16 are (i) that electrical contacts must be deposited and contacted, and (ii) the technical difficulties in studying electrical conductors, since the transducer must be electrically insulated from the sample by an intermediate insulating layer to avoid leakage of the electrical current used to heat the resistor. On the other hand, contactless techniques such as those reported in Refs. 9 and 17–27 are based on small Gaussian or ellipse-shaped focused spots, i.e., sensitive to all crystallographic directions simultaneously. Although this is not an intrinsic impediment to obtain κ_ij, it substantially complicates the analysis of the measured data with respect to the case of an elongated line-shaped heat source (sensitive to only two crystallographic directions simultaneously), as it is evident from Refs. 9 and 17–27. In addition, for Gaussian or ellipse-shaped small focused spots, precise determination of the spatial intensity distribution of the focused laser spot is required, which in some situations can be challenging due to the highly asymmetric shape of the heater spot.^22,23 We note that for elongated line-shaped geometries of the focused spot, the precise shape of the spot must not be necessarily known to obtain the thermal properties of the samples (“slope method”). The main objective of the present work is to provide a new technical approach that combines the key advantages of the experimental methods reported in Refs. 9 and 15–27 but without their drawbacks, i.e., based on a one-dimensional heat source configuration, which is focused on the surface of the samples in a contactless fashion.

Here, we developed a novel experimental approach to measure the thermal conductivity tensor in a contactless fashion, which is advantageous to study isotropic and anisotropic materials and thin films. This technique, which we refer to as “Anisotropic Thermoreflectance Thermometry” (ATT), is based on the geometry used in the 3-omega method, but instead of using an electrical resistor, we use a focused laser beam to define a long and narrow line-shaped heat source, hence taking profit from the many advantages offered by this geometry, such as, e.g., the determination of the thermal conductivity through the so-called “slope method.” We derive the fundamental equations that lead to the solution of the problem for bulk, thin films, and multilayer systems. The excellent experimental performance of this methodology is demonstrated in isotropic bulk materials such as glass, Si, and Ge, as well as for the case of thin films (SiO₂ thin films on Si substrates). Finally, we address its applicability to study the thermal conductivity tensor in anisotropic materials such as β-Ga₂O₃ and Kapton tape, showing that the present approach requires minimum sample processing and data modeling, but keeping at the same time a rather high experimental accuracy and temperature sensitivity.

In this section, we derive the solution for the temperature response of the system upon a line-shaped thermal excitation. In particular, we provide the mathematical expressions necessary to obtain the thermal conductivity for different geometries of the sample from the frequency-dependent thermal response. The problem we aim at studying is the temperature response of a heterostructure subjected to a long and narrow line-shaped heat source. This geometry is similar to that used in the 3-omega method,^28,31 where a narrow (≈10 μm) and long (≈1 mm) metallic resistor is deposited on the surface of the sample, which is used simultaneously as the heater and thermometer. However, in our case, the heater is optically defined by focusing a Gaussian laser on a line-shaped geometry with uniform power distribution, which is achieved by diffractive beam profile reconstruction. A schematic illustration of the studied geometry is shown in Fig. 1, and the technical details are provided in Sec. III. The intensity distribution of the heat source in the direction across the line is Gaussian, whereas it is almost uniform along the long axis (flat-top profile). We start by providing the solution for the semi-infinite medium case, which we generalize for the case of an anisotropic multilayer system using the formalism developed by Borca-Tasciuc et al. in Ref. 31. We also provide the limit for a large thermal penetration depth, which leads to the so-called “slope method” for a semi-infinite anisotropic sample. We recall that this method is particularly attractive because it allows an accurate and rather simple determination of the thermal conductivity since κ ∝ [∂T/∂  ln(ω)]⁻¹; hence, a linear fit of the experimental data is sufficient to obtain the thermal conductivity with no additional assumptions and/or numerical modeling. Finally, the solution for the case of a thin film on a semi-infinite substrate is also provided. We note that although the mathematical treatment of this problem is similar to that presented in Refs. 28 and 31, in our case, the heat source intensity distribution in the direction perpendicular to the long axis is Gaussian and not a step function as in the case of the 3-omega method,²⁸ which must be taken into account in order to obtain the correct thermal response of the system.

We summarize the equations that apply to each case depending on the geometry of the sample as follows, including a short description of the key assumptions:

1. Isotropic substrate (general solution) → Eq. (7):

   The solution can be obtained by numerical integration and is valid for the whole frequency range provided that the substrate can be considered as semi-infinite.

2. Isotropic substrate (Λ_p = |q|⁻¹ ≫ σ) → Eq. (10):

   The solution is valid when the thermal penetration depth (|q|⁻¹) is much larger than the width, σ, of the Gaussian heat source. The thermal conductivity is obtained from the slope method with κ = [P_p/(2πl_p)](∂T/∂ ln(ω))⁻¹.

3. Anisotropic substrate (general solution) → Eq. (11):

   The solution can be obtained by numerical integration and is valid for the whole frequency range provided that the substrate can be considered as semi-infinite.

4. Anisotropic substrate (Λ_p = |q|⁻¹ ≫ σ) → Eq. (15):

   The solution is valid when the thermal penetration depth (|q|⁻¹) is much bigger than the width, σ, of the Gaussian heat source. The effective thermal conductivity is obtained from the slope method as κ_eff = [P_p/(2πl_p)](∂T/∂ ln(f))⁻¹, and it is related to the components as $κeff=κ⊥κ∥$⁠.

5. Thin films on substrate → Eq. (16):

   This expression is valid when the thermal conductivity of the substrate is larger than the thermal conductivity of the thin film; rule of thumb → κ_substrate ⪆ 10κ_TF.

Figure 2 displays the numerical solution of Eq. (7) for the case of glass and Si substrates as well as the numerical solution of Eq. (11) for a SiO₂/Si thin film using 10 mW of heater power. We plot the frequency-dependent maximum temperature rise obtained using the present method (ATT) as well as the corresponding solution for the 3-omega geometry, which would correspond to a hypothetical illumination beam with a square section in the sort axis [see the inset of Fig. 2(a)]. Figure 2(a) displays the case of a glass substrate coated with 60 nm of an Au transducer. As expected, the temperature response approaches a linear dependence (in a logarithmic frequency scale) for rather low frequencies (large thermal penetration depth). Within the lower frequency range, the thermal conductivity can be directly obtained through the “slope method,” i.e., κ ∝ [∂T/∂ ln(ω)]⁻¹. For higher frequencies, the thermal penetration depth decreases and becomes comparable with the width of the heat source, hence leading to deviations from the linear behavior (in the logarithmic frequency scale) observed at lower frequencies. Figure 2(b) displays the amplitude of the temperature oscillations for a Si substrate and for a 285 nm thick SiO₂ thin film deposited on a Si substrate. For Si, the maximum temperature depends linearly on ln(ω), which is a direct consequence of the large thermal penetration depth originated from its large thermal conductivity (or, equivalently, thermal diffusivity $→α=κCp−1ρ−1$⁠). It is interesting to note that in all studied cases, the temperature response is similar to that obtained from the 3-omega method (dashed lines in Fig. 2); however, for similar excitation spot sizes, the magnitude of the temperature response of the present approach is always larger as compared to the case of the 3-omega method. This originates from the fact that the power density is larger by a factor $8/π$ for a Gaussian intensity distribution with radius σ, as compared to a step-function distribution with half width, b, when σ ≈ b. In conclusion, the thermal response obtained through the herein studied geometry (long line with a Gaussian cross section) is similar to that obtained through the geometry used in the 3-omega method (long line with a square cross section). The only difference between these approaches in terms of the temperature response is that the effective area of the heater region is different in each case (Gaussian vs square cross section). However, note that the key advantage of our approach is that it allows us to access an arbitrary in-plane direction due to its contactless fashion.

The experimental setup used in the present experiments is shown in the schematic illustration of Fig. 3. The experiments are based on a frequency-domain pump-and-probe concept, where the wavelength of the pump laser was set to λ^pump = 405 nm (Omicron A350) and the probe laser wavelength was λ^probe = 532 nm (Cobolt SAMBA series), whereas the pump laser power was modulated with an external harmonic wave generator (Rigol DG5332) between 63 Hz and 100 kHz, and the incident probe laser intensity was constant in time. The power of the probe laser was kept low in order avoid any heating effects, typically of the order of hundreds of μW. On the other hand, the pump laser power was increased until a temperature rise of the order of several degrees was observed, depending on the thermal conductivity of the sample. The typical power used for the heater laser was in all cases ≈50 mW; however, this quantity was precisely measured for each sample and optical alignment conditions. Both the lasers were coupled into the main optical axis using several beam splitters (BSs) and dichroic beam splitters. The diameter of the pump and probe beams at the laser output was 1100 and 800 μm, respectively, and their intensity distribution was Gaussian. Two beam expanders were placed in the optical path of the lasers because the beam diameter must be adjusted precisely for successful diffractive beam reconstruction, which was achieved through a custom designed holographic diffractive optical element (DOE) purchased from HOLO/OR Ltd. This element was placed just before an achromatic focusing lens (f = 50 mm). The DOE is one of the key elements of the experimental setup because it is responsible for homogenizing the intensity of the pump and probe lasers, i.e., the Gaussian distribution is converted through diffraction into a line-shaped spot with uniform intensity distribution along the long axis and with Gaussian distribution in the perpendicular direction. Figure 3 also displays an optical image of the pump and probe lasers after intensity reconstruction using the DOE. We note that although there are several ways to convert a Gaussian beam into a flat-top distribution, such as, e.g., using a fiber with a square core to induce mode mixing or simply using an optical diffuser in combination with a cylindrical lens, we have observed that using a DOE leads to the best results in terms of spot shape and size, and output pump laser power. The back reflection of the lasers after focusing on the sample was coupled into the detection arm using a BS, and the laser wavelength was selected by two mechanically controlled notch filters. As the detector, we used a large area (5 × 5 mm²) AC-coupled balanced detector (Thorlabs PDB210A/M-AC-SP). A sample of the probe laser (532 nm) was independently coupled to the second detector input. The output of the AC-coupled balanced detector is sensitive to the difference between both the inputs. The purpose of such a detector is reducing the intrinsic noise within the probe laser, typically 0.2%. When the sample is excited by the harmonic pump laser, the output of the detector was a frequency modulated signal arising from the thermal oscillations on the surface of the sample. We note that using such a balanced detector scheme is convenient mostly in the lower frequency range, f < 5 kHz, because the total noise is dominated by f⁻¹ noise. In order to minimize the laser noise of the optical path of the signal, reference lasers were kept similar within 10 cm, and the input power on both the detectors was balanced. The pump laser was also coupled to the detection system in order to obtain a reference phase for the thermal signal, i.e., first, the pump laser signal (amplitude and phase) was measured by the lock-in amplifier after which the mechanical notch filter was switched, and the amplitude and phase of the thermal signal were obtained. This process was sequentially repeated for each excitation frequency.

In order to convert the voltage measurement (V units) to temperature (K units), several factors must be taken into account. First, the amplitude of the reflectivity oscillations measured with the lock-in amplifier (V_AC) must be normalized by the DC reflectivity at a probe wavelength of 532 nm (V_DC). This process must be carefully done because the balanced detector and the low-noise amplifier (LNA) involve several amplification stages of the modulated signal (thermal signal). Note that whereas V_DC is measured by the DC-coupled monitor outputs of the detector (not amplified), V_AC is obtained after signal subtraction with the reference optical input (RF output), which involves an internal amplification stage. In addition, V_AC was also externally amplified by a low-noise amplifier (LNA) to 40 dB (two orders of magnitude). The AC (thermal signal) to DC (probe signal) amplification factor was determined for each setting of the external LNA and ranges between 80× and 8000×. Hence, after correcting V_AC by the amplification factor, the resulting voltage ratio can be related to the optical reflectivity as V_AC/V_DC = ΔR/R₀. Finally, the reflectivity change, ΔR/R₀, is converted to temperature by using the temperature coefficient of reflectivity of the Au transducer. The photothermal response of the Au transducer was calibrated studying a glass control sample whose thermal conductivity was independently measured using the 3-omega method, obtaining a thermal conductivity of 0.9 W m⁻¹ K⁻¹. This value was used as the input calibration to obtain the temperature coefficient of the reflectivity of the Au transducer of the same glass sample through the application of the slope method using Eq. (10). We have obtained (1/R₀)[∂R/∂T] = 3 × 10⁻⁴ K⁻¹ in good agreement with the most reported values for Au.³⁴ So far, we have not observed an appreciable variation of this coefficient between samples fabricated in different Au evaporation cycles because we always use the same evaporation conditions. We have estimated the temperature uncertainty of the method using a glass substrate coated with a 60 nm thick Au transducer by measuring the minimum detectable signal as compared to the calibrated temperature rise in the samples. We have found that it is possible to measure temperatures as low as ≈3 mK. This is possible mainly due to the balanced detection scheme used for signal detection, which provides a noise floor suppression of more than two orders of magnitude (50 dB). Regarding the accuracy in the determination of the thermal conductivity, we have estimated an experimental accuracy of $≈5%$⁠, which was statistically obtained by 30 successive measurements. We recall that the glass sample was independently measured using the 3-omega method. A detailed study of the various error sources and of the sensitivity for different frequencies will be considered for future study.

Regarding the phase lag of the optical reflectivity with respect to the thermal harmonic excitation, it can be shown by solving Eq. (7) that for substrates with rather low thermal conductivity such as glass, the phase lag is almost independent of the thermal boundary resistance (TBR) between the Au transducer and the substrate. However, for substrates with higher thermal conductivity such as Si, the phase lag is affected by the TBR, and hence, this method can be also used to determine the TBR. Note that in both the previous cases, the geometry of the heat source also affects the phase lag, and thus, its dimensions must be precisely measured. With respect to the absolute value of the temperature rise, it can be also shown through the solution of Eq. (7) that its logarithmic frequency derivative remains approximately unaffected, as expected for the “slope method.” Finally, the thermal boundary resistance between layers was modeled following the approach developed by Cahill,³⁰ where the heat diffusion equation is solved in the presence of a multilayered system, and each layer has a finite thermal conductivity, heat capacity, and thickness. The case of an interface is modeled assuming a rather low heat capacity and a thickness of 1 nm, as described in Ref. 30. A comprehensive review of the different models that can be used to account for the thermal boundary resistance can be found in Ref. 35.

We have studied five bulk samples of glass, such as Si(100), Ge(100), Ga₂$O3(2̄01)$⁠, Kapton tape, and a 285 nm thick SiO₂ thin film deposited on a Si(100) substrate (10–20 Ω cm). All samples were coated with a 60 nm thick Au thin film (transducer), which was thermally evaporated at a base pressure of 10⁻⁶ mbar. All measurements were done under ambient conditions. Of particular interest, we investigated a $(2̄01)$ oriented Ga₂O₃ substrate, which exhibits strong in-plane thermal anisotropy. We show in this section that the present method provides a straightforward means to determine the in-plane components of the thermal conductivity tensor.

Figure 4 displays the experimental measurements of the amplitude of the temperature oscillations as a function of the excitation frequency for the case of glass, Si, Ge, and for a 285 nm thick SiO₂/Si thin film. For the glass substrate, shown in Fig. 4(a), the temperature decays approximately linearly with ln(ω). However, as the frequency increases, the frequency-dependent temperature rise gradually deviates from the linear behavior observed at lower frequencies, in good agreement with the calculations shown in Fig. 2. We recall that the linear dependence observed at lower frequencies allows the determination of the thermal conductivity directly through the “slope method,” i.e., avoiding complex data modeling and numerical fitting of the experimental data. This is a key aspect because in order to quantitatively model the experimental data, it is necessary to know the precise width of the line-shaped optical heater. However, the slope of the linear dependence observed at lower frequencies is independent to some extent of the width of the heater. A linear fit to the data points in the frequency range up to ≈1 kHz results in a thermal conductivity of 0.9 ± 0.05 W m⁻¹ K⁻¹.

In Fig. 4(b), we show the frequency-dependent temperature response for the case of Si, Ge, and a 285 nm thick SiO₂ thin film deposited on a Si(100) substrate. Note that the used frequency range is different as compared to the case of glass [see Fig. 4(a)], because for the lower frequency range (f ⪅ 1 kHz), the thermal penetration depth for Si and Ge becomes comparable with the thickness of the substrates; hence, the thermal response is influenced by the sample holder. The amplitude of the temperature rise is observed to be dependent linearly on ln(ω) for Si and Ge, which resembles the observations on the glass sample at lower frequencies. However, for Si and Ge, the thermal penetration depth, $Λp=1/|q|=α/(ω)$⁠, is substantially larger as compared to the case of glass since $Λp∝α$⁠, where we recall that α = κ/(C_pρ) is the thermal diffusivity of the sample. The thermal conductivity of the samples was obtained using the slope method as 137 ± 7 and 52 ± 3 W m⁻¹ K⁻¹ for Si and Ge, respectively, in good agreement with the generally accepted values for Si and Ge.³⁶ 

The frequency dependence of the temperature rise for the 285 nm thick SiO₂/Si thin film is similar to the case of pure Si as shown in Fig. 4(b), although a temperature offset arising from the presence of the SiO₂ thin film is observed. We recall that, similarly to the case of the 3-omega method, no frequency dependence is introduced by the presence of the thin film in the temperature rise. In other words, if the thermal penetration depth is much larger than the thickness of the thin films, no frequency dependence of the temperature rise is expected. The thermal conductivity of this sample was 1.2 ± 3 W m⁻¹ K⁻¹, which is a typical value for SiO₂/Si thin films.

Finally, we demonstrate the strengths of the present approach for the study of the thermal properties of anisotropic materials. As a model system, we chose the sesquioxide β-Ga₂O₃, which is the most stable of several Ga₂O₃ polymorphs and crystallizes in the monoclinic structure (space group C2/m) with pronounced optical, electrical, and thermal anisotropies. Here, we present measurements of a single crystal β-Ga₂O₃ substrate with the $(2̄01)$-plane. Figure 5(a) displays the temperature rise as a function of the excitation frequency for two perpendicular in-plane directions corresponding to the [010] and ⊥[010] crystallographic directions. A different temperature response is observed for these two directions, with a lower temperature rise along the [010] as compared to ⊥[010]. The temperature response is linear in the logarithmic scale up to ≈5 kHz, whereas for higher frequencies, a similar deviation as previously discussed for the case of glass is observed. Thus, the anisotropic thermal conductivity is once more obtained directly by the “slope method” in this frequency range, i.e., without the need for numerical fitting.

Finally, in Fig. 6, we show the optical and thermal response of a piece of Kapton tape. We aim at demonstrating that even if the quality of the surface is far from ideal (arising in this case from the fabrication process), it is still possible to determine in-plane thermal anisotropy. In particular, the fabrication process of Kapton tapes involves stretching of the material in one direction. Figure 6(a) displays a polarized optical image of the sample, where surface inhomogeneities, roughness, and bubbles can be appreciated. Moreover, the principal optical axis can be easily identified by optical polarization. All images were taken using a crossed polarization configuration for different rotation angles of the sample, as shown in the inset of Fig. 6. Accordingly, the thermal anisotropy scan exhibits anisotropy in a direction coincident with the optical axis of the material. We note that the dispersion of the points within Fig. 6(b) does not only arise from experimental errors (indicated in Fig. 6) but also from defects inherently present at the surface of the sample, which can affect the temperature distribution.

We have developed a novel contactless experimental approach suitable to study anisotropic thermal transport based on thermoreflectance thermometry. We demonstrated its applicability to the case of isotropic (Si and Ge), anisotropic (β-Ga₂O₃), and Kapton bulk samples, as well as to thin films (SiO₂/Si). The experimental accuracy of the method in the determination of the thermal conductivity is ≈5%, whereas the temperature uncertainty is ΔT ≈ 3 mK. The method combines all the advantages offered by non-contact approaches, but keeping mathematical simplicity in the data evaluation process. In fact, a simple linear fit to the frequency-dependent temperature response (“slope method”) is sufficient to obtain the thermal conductivity of the samples, similar to that in the case of the 3-omega method. The key advantage of this approach from the thermal perspective lies in the geometry of the heat source, which provides simultaneous sensitivity to two crystallographic directions (perpendicular to the heater line), whereas it is insensitive to the direction parallel to the heater. This ability, together with its contactless implementation, allows us to probe all in-plane directions by continuous rotation of the sample. Hence, the method is particularly convenient for measuring the projection of the thermal conductivity tensor onto the plane corresponding to the surface of the samples. We envisage that the present experimental approach provides a new alternative to accurately determine the thermal conductivity tensor in anisotropic materials.

The authors acknowledge financial support from the Spanish Ministerio de Economía, Industria y Competitividad, through Grant No. CEX2019-000917-S (FUNFUTURE) in the framework of the Spanish Severo Ochoa Centre of Excellence program and Grant Nos. PID2020-119777GB-I00 (THERM2MAIN), PDC2021-121814-I00 (COVEQ), PGC2018-095411-B-100 (RAINBOW), MAT2017-90024-P (TANGENTS)-EI/FEDER, and 2020AEP141; the Generalitat de Catalunya through Grant Nos. 2017SGR488 and AGAUR 2018 PROD 00191; and from the European Research Council (ERC) under Grant Agreement No. 648901. K. X. acknowledges a fellowship (No. CSC20180695006) from the China Scholarship Council and the Ph.D. program in Materials Science from Universitat Autònoma de Barcelona, in which he was enrolled.

The authors declare no conflict of interest.

L.A.P. and K.X. contributed equally to this work.

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