In this work, we develop a numerical fitting routine to extract multiple thermal parameters using frequency-domain thermoreflectance (FDTR) for materials having non-standard, non-semi-infinite geometries. The numerical fitting routine is predicated on either a 2D or 3D finite element analysis that permits the inclusion of non-semi-infinite boundary conditions, which cannot be considered in the analytical solution to the heat diffusion equation in the frequency domain. We validate the fitting routine by comparing it with the analytical solution to the heat diffusion equation used within the wider literature for FDTR and known values of thermal conductivity for semi-infinite substrates (SiO2, Al2O3, and Si). We then demonstrate its capacity to extract the thermal properties of Si when etched into micropillars that have radii on the order of the pump beam. Experimental measurements of Si micropillars with circular and square cross sections are provided and fit using the numerical fitting routine established as part of this work. Likewise, we show that the analytical solution is unsuitable for the extraction of thermal properties when the geometry deviates significantly from the standard semi-infinite case. This work is critical for measuring the thermal properties of materials having arbitrary geometries, including ultra-drawn glass fibers and laser gain media.

Measurements of thermal transport properties in nanoscale thin films are conventionally made using optical pump–probe thermoreflectance techniques, principally due to the non-contact nature with which they are able to interrogate nanoscale thermal transport characteristics for nearly any material type.1 In contrast to Raman spectroscopy2 and 3-ω3 techniques, thermoreflectance-based measurements can separate the impacts of thermal boundary conductance (G) across interfaces and thermal conductivity (κ) within individual material layers4 and have overwhelmingly served as the thermal characterization technique of choice for nanoscale material systems over the course of the last decade.1,5–8 However, current models used to extract the thermal properties of nanoscale materials limit the geometries that can be interrogated to those that are (1) semi-infinite in the radial direction and (2) have finite or semi-infinite thickness transverse to the direction of the applied heat source.9 In this work, we establish a finite element-based numerical fitting routine to extend the utility of thermoreflectance techniques for use with any planar geometry having finite dimensions.

The two thermoreflectance systems mostly used over the course of the last decade include time-domain thermoreflectance (TDTR)10 and frequency-domain thermoreflectance (FDTR).11 Both techniques rely on two separate laser beams to (1) heat the sample surface and (2) probe the reflectivity (i.e., temperature) on the sample surface. The beams that heat and probe the surface are referred to as the “pump” and “probe” beams, respectively. Typically, a metal transducer (50–150 nm of Au or Al) is deposited on the sample surface to convert the optical energy of the pump beam to thermal energy prior to the sample surface and due to a well-established wavelength-dependent relationship between the reflectivity of the transducer and its surface temperature (often referred to as the coefficient of thermoreflectance).12,13 The pump beam is modulated at a single frequency (TDTR) or across a range of frequencies (FDTR) such that we can use a lock-in amplifier to detect small changes in the reflectance as the heat penetrates into the sample. In TDTR, a pulsed laser source is used and split into two different pump and probe paths. A delay stage is used to physically delay the arrival of the probe beam relative to the arrival of the pump beam to monitor changes in reflectance (i.e., temperature) at the sample surface over time. On the other hand, FDTR utilizes two separate continuous wave (CW) lasers to establish pump and probe beams, where the pump beam is modulated over a range of frequencies. Modulating frequency allows for corresponding changes to the penetration depth of the heat deposited by the pump beam and thus establishes sensitivity to multiple thermal properties and/or the thermal properties of several underlying layers of material in a multilayer stack.9 

Recent advances in thermoreflectance-based techniques include extensions to a steady-state system to gain sensitivity to the so-called “buried interfaces,”14 the use of a magneto-optical Kerr effect (MOKE) to gain sensitivity to interfaces having large thermal conductance15 and phonon-magnon coupling effects16 and the development of a transient thermo-transmission technique to measure the thermal boundary conductance of nanoparticles suspended in transparent media.17 Recently, transient grating spectroscopy has also been used to probe non-diffusive thermal regimes in films that are suspended over micrometer-sized regions.18–20 However, conventional measurements are still limited to geometries that are semi-infinite, which still limits the utility of the technique to material systems that can be fabricated in such a way. For instance, these characterization techniques are blind to the thermal properties of material systems whose thermal properties are expected to change with geometry, such as ultra-drawn glass fibers (i.e., fiber-optic components), strained polymers, and laser gain media. In this work, we integrate a finite element-based numerical simulation (constructed in COMSOL Multiphysics v. 5.5) into the conventional fitting routine for frequency-domain thermoreflectance measurements such that the thermal properties of nanoscale and microscale material systems having non-standard geometries can be accurately extracted. We take measurements on standard substrates having well-known thermal properties (SiO2, Al2O3, and Si) to validate the numerical model against values obtained analytically and to those that exist within the wider literature. We then demonstrate the difference in the numerical solution to the phase lag (i.e., temperature response) at the sample surface due to changes in the radial geometry of Si in the form of Si micropillars with varying height. Finally, we obtain the thermal conductivity of Si when in micropillar form using the numerical fitting routine and an experimental measurement made using FDTR.

In this work, FDTR is used to characterize the thermal properties of bulk substrates and Si micropillars. Previous studies provide a detailed description of the FDTR system used in this effort.21,22 Briefly, FDTR measures the phase lag at the sample surface relative to the imposed phase applied by the modulated heating event. One useful analogy (despite occurring in the time domain) is the lag in the temperature response of water relative to the temperature of a stove-top’s heater. Because our experiment applies a sinusoidal modulation to the heating event, we can observe the temperature response at the sample surface by measuring the phase lag at the sample surface. A representative schematic of the mechanism used to track temperature on the sample surface is provided in Fig. 1.

FIG. 1.

(a) Time-domain representation of modulated heat source (blue) and temperature response (green), where the peak-to-peak difference is the phase lag measured by the lock-in amplifier (ϕLI) and (b) phase lag vs measured beam intensity. We note that the temperature difference between the probe beam and the probe reference (i.e., the temperature rise/fall at the transducer surface in response to the modulated heating event) is very small, requiring the use of a lock-in amplifier).

FIG. 1.

(a) Time-domain representation of modulated heat source (blue) and temperature response (green), where the peak-to-peak difference is the phase lag measured by the lock-in amplifier (ϕLI) and (b) phase lag vs measured beam intensity. We note that the temperature difference between the probe beam and the probe reference (i.e., the temperature rise/fall at the transducer surface in response to the modulated heating event) is very small, requiring the use of a lock-in amplifier).

Close modal

In Fig. 1(a), the blue curve represents the power of the pump laser, which is modulated between an upper and lower power (1 and 1, respectively). The corresponding temperature response of the surface is shown in green, which lags behind the applied temperature and fluctuates between the upper and lower temperature (again 1 and 1, respectively). The difference in the phase between the pump and the probe is the phase lag measured by the lock-in amplifier, ϕLI, and this provides information on the thermal properties of the material when used in tandem with a well-established multilayer analytical model.11 

Figure 1(b) displays the intensity of the pump beam as well as the intensity of the reflected probe. The intensity and apparent modulation of the probe beam comes from the temperature response of the material induced by the pump. The intensity of the reflected probe depends on the surface reflectivity, which is proportional to the change in temperature. The probe reference that did not interact with the surface—and therefore remains unmodulated—is shown by the dotted black line. By subtracting the probe signal from the reference and taking into account the coefficient of thermoreflectance, the temperature response of the sample can be determined.

A schematic of the FDTR system used in this work is shown in Fig. 2. The system contains two separate continuous-wave lasers that act as the “pump” (405 nm Coherent OBIS CW laser) and the “probe” (532 nm Coherent OBIS CW laser) respectively. We note that at these wavelengths, the transducer absorbs a significant amount of the pump beam and we are highly sensitive to changes in the reflectivity of the sample surface by the probe beam.

FIG. 2.

Schematic of the frequency-domain thermoreflectance system built at USNA. Note the following acronyms: PBS (polarizing beam splitter), WP (waveplate), and EOM (electro-optic modulator).

FIG. 2.

Schematic of the frequency-domain thermoreflectance system built at USNA. Note the following acronyms: PBS (polarizing beam splitter), WP (waveplate), and EOM (electro-optic modulator).

Close modal

The pump beam is first split into two separate paths: (1) through the electro-optic modulator (EOM, Conoptics Model 350-160) and (2) into the balanced photodetector (Thorlabs, PDB450A-AC) that tracks the reference phase via a 1% beam splitter. The small amount of pump beam that is immediately directed into the photodetector is used as a reference to subtract any coherent noise within the laser. The portion of the pump beam that is routed through the EOM is modulated using the built-in waveform generator in our lock-in amplifier (Zurich UHFLI). After exiting the EOM, the pump beam is again split into separate “primary” (downward direction) and “reference” (leftward direction) paths using a polarizing beam splitter. The primary path is steered into an objective lens (Mitutoyo 50x) using a dichroic mirror and focused onto the sample surface. The reference path is used to track the applied phase of the pump beam at the sample surface. In order to match the phase at the sample surface, we modulate at the highest frequency used in our measurement (ωmax=20MHz) and match the numerical phase of the pump beam that leaks into the primary balanced photodetector (shown on the left of the image) while blocking the probe beam. This allows us to measure the phase response at the sample surface (measured with the probe beam) and the imposed phase at the sample surface (applied via the pump beam) simultaneously, which is unique to our implementation of FDTR. By subtracting the two signals, we obtain ϕLI as described previously.

Provided with a measurement of ϕLI, one can obtain the underlying thermal properties across an interface (i.e., the thermal boundary conductance, G) or within an individual layer of a multilayer stack (e.g., the thermal conductivity, κ). The entire formulation of the analytical expression used to extract thermal properties is described in detail elsewhere.11 However, it is useful to describe several features of the analytical solution to the frequency-domain version of the heat diffusion equation to demonstrate its geometric limitations. The analytical solution is constructed within the framework of a multilayer material stack whose substrate is semi-infinite in the radial direction, and most often semi-infinite in the through-thickness direction. A schematic of the general multilayer material stack used in the development of the analytical solution is provided in Fig. 3.

FIG. 3.

Arbitrary multilayer material stack with semi-infinite boundary conditions in r- and z-directions. Pump (blue) and probe (green) beams are depicted above the Au transducer.

FIG. 3.

Arbitrary multilayer material stack with semi-infinite boundary conditions in r- and z-directions. Pump (blue) and probe (green) beams are depicted above the Au transducer.

Close modal

In the time domain, the equation governing heat diffusion is expressed as

Cvθt=κrrr(rθr)+κz2θz2.
(1)

The above expression accounts for 2D heat flow in the radial (r) and through-plane (z) directions. A Fourier transform is applied to obtain the frequency-dependent heat diffusion equation, which is written as

κz2θ(ω,k,z)z2=(κrk2+Cviω)θ(ω,k,z),
(2)

where q is defined for a layer of material n and thickness d as

q2=κrk2+Cviωκz.
(3)

We can relate the temperature to the heat flux at the top surface (subscript t) of a slab made of a certain material in the frequency domain with the bottom surface (subscript b) using

[θn,bfn,b]=[cosh(qd)1κzqsinh(qd)κzqsinh(qd)cosh(qd)][θn,tfn,t].
(4)

The temperature and heat flux between the bottom surface of material n are connected to the top of material n+1 via

[θn+1,tfn+1,t]=[1G101][θn,bfn,b],
(5)

where G is the thermal boundary conductance between the two layers. The heat flux boundary condition of the top, ft, can be found with

ft=A02πexp(k2wo28),
(6)

which is the Hankel transform of a Gaussian spot with a power of A0 and a 1/e2 radius of w0. If there are multiple layers, the solution can be found with

[θbfb]=MnMn1M2M1=[ABCD][θtft],
(7)

where Mn is the matrix of the bottom layer. If the bottom layer is treated as adiabatic or semi-infinite, the surface temperature can be found using

θt=DCft.
(8)

The final frequency response, H(ω), is found by taking the inverse Hankel transform of Eq. (2) and weighting it with a Gaussian spot with a 1/e2 radius of w1,

H(ω)=A02π0k(DC)exp(k2(w02+w128)dk.
(9)

The thermal model for H(ω) is then fitted to the lock-in phase data. By changing the parameters of the thermal model to fit the lock-in data, the thermal properties can be determined. The lock-in phase data measured are given by

ϕLI=tan1(H(ω))(H(ω))+ϕext,
(10)

where (H(ω)) is the out-of-phase signal, (H(ω)) is the reference signal, and ϕext is the external phase shift caused by other aspects not caused by changes in reflectively, such as the optical path length, driving electronics, and photodetectors. Thermal properties are extracted by fitting Eq. (10) to measurements of the phase lag via FDTR for SiO2, Al2O3, and Si, shown in Fig. 4.

FIG. 4.

Measured data (open black circles) and corresponding analytical fits for SiO2 (solid blue line), Al2O3 (solid red line), and Si (solid purple line) substrates. Pump and probe diameters are 5.7 μm and 3.7 μm, respectively.

FIG. 4.

Measured data (open black circles) and corresponding analytical fits for SiO2 (solid blue line), Al2O3 (solid red line), and Si (solid purple line) substrates. Pump and probe diameters are 5.7 μm and 3.7 μm, respectively.

Close modal

The thermal properties obtained above are consistent with those widely reported in the scientific literature.14 The thermal boundary conductance for all three samples is high relative to other values reported in the literature;14 however, we utilize an 5nm Ti layer between the Au transducer and the bulk substrates, which improves adhesion and therefore enhances thermal transport across the interface.

To determine whether a particular thermal property (κ, Cv, and G) can be extracted across the range of modulation frequencies applied on the sample surface, we can determine the sensitivity of each parameter to small perturbations in the measured phase lag. The phase sensitivity to a particular thermal property, x, at a given frequency can be found with

S(ω)=ϕ(ω)lnx.
(11)

The sensitivity to the thermal properties shown in Fig. 4 is provided in Fig. 5.

FIG. 5.

Sensitivity (S) to the extracted thermal properties (where each property is represented by the subscript “x”) for each substrate shown in Fig. 4.

FIG. 5.

Sensitivity (S) to the extracted thermal properties (where each property is represented by the subscript “x”) for each substrate shown in Fig. 4.

Close modal

A two-dimensional (2D) finite element model is developed to directly fit for the thermal properties of geometries with non-semi-infinite boundary conditions. (Note that the reader can find the full Matlab fitting routine and details of its link to COMSOL in the supplementary material.) The most palatable system to demonstrate the utility of the finite element model is one in which the geometry is confined in the radial direction and one that has been widely fabricated in laboratory environments. Consequently, we choose to develop the model based on Si micropillar arrays. The arrays we interrogate have geometries that vary in both the radial and through-plane directions on the order of single-digit μm to 10s of μm. As these lengths remain larger than the mean free path of phonons in Si,23 the geometric confinement should not result in any change in κSi. However, it is likely that the phase lag (ϕLI) does change due to the confinement of heat (recall that the phase lag represents the response in the temperature on the top of the transducer relative to the modulated signal of the applied heating event).

The finite element model is created in COMSOL Multiphysics v. 5.5. We note here that COMSOL is used with its LiveLink module in order to communicate with Matlab. Without the LiveLink module, COMSOL cannot communicate with Matlab, which would not allow us to build the simulation or use COMSOL’s iterative solver to converge on a solution numerically. The Si micropillar is modeled using the computational domains shown in Fig. 6. We note that we incorporate the Si substrate and the Si micropillar within our model. An 80 nm Au transducer is constructed above the pillar and a finite value for thermal boundary conductance is applied at the Au/Si micropillar interface (in our model, this remains a free parameter). The Si micropillar and the Si substrate are “continuous” in the sense that there is no applied thermal boundary conductance at the interface between domains. This is physically appropriate given the nature of the fabrication process; the Si micropillars themselves are etched and are never physically separated/reattached during the process.

FIG. 6.

Schematic of the numerical model built in COMSOL Multiphysics. hpillar represents the pillar height, tSi is the thickness of the Si substrate (approximated as semi-infinite, but modeled as 500 μm in the computational domain), rpillar is the radius of the micropillar, rSi is the radius of the Si substrate (approximated as semi-infinite, but modeled as 300 μm), and w0 and w1 are the pump and probe radii, in μm. Purple boundaries represent those boundaries that are insulated, blue boundaries represent boundaries which are held at a constant initial temperature, T0= 300 K, and red boundaries are symmetric.

FIG. 6.

Schematic of the numerical model built in COMSOL Multiphysics. hpillar represents the pillar height, tSi is the thickness of the Si substrate (approximated as semi-infinite, but modeled as 500 μm in the computational domain), rpillar is the radius of the micropillar, rSi is the radius of the Si substrate (approximated as semi-infinite, but modeled as 300 μm), and w0 and w1 are the pump and probe radii, in μm. Purple boundaries represent those boundaries that are insulated, blue boundaries represent boundaries which are held at a constant initial temperature, T0= 300 K, and red boundaries are symmetric.

Close modal

COMSOL features a built-in frequency-domain perturbation solver. In this case, COMSOL is solving the frequency-domain representation of the heat equation in the following form:

iωρCpT+q=Q+Qp,
(12)

where ω is the applied angular frequency, ρ is the material density, Cp is the specific heat capacity, T is the temperature (where ρ, Cp, and T are assumed to be at constant pressure), q is the vector quantity described by kT (assuming here that k is isotropic, though it does not need to be in COMSOL), Q is the average value of the heat load applied to the surface, and Qp is the harmonic perturbation around Q. This allows us to compute the harmonic temperature variations around an equilibrium temperature in response to the selection of our boundary heat source as a perturbation load. The use of the above expression results in a linear quasi-steady problem in the frequency domain, which is less computationally expensive than the equivalent problem in the time domain.

The numerical model itself is meshed using a graded grid in each independent sub-domain shown in Fig. 6, including the Au transducer layer, the Si micropillar, the region immediately below the Si micropillar (i.e., from r = 0 to r=rpillar in the Si substrate), and the remainder of the Si substrate. An image of the meshed sub-domains is provided in Fig. 7. Note that a mesh independence study was completed to ensure that the solution was independent of the number of nodal points in each sub-domain.

FIG. 7.

Mesh for Si micropillar on the Si substrate with the Au transducer. The entire domain contains 64 000 elements split equally among four distinct sub-domains (numbered 1–4).

FIG. 7.

Mesh for Si micropillar on the Si substrate with the Au transducer. The entire domain contains 64 000 elements split equally among four distinct sub-domains (numbered 1–4).

Close modal

As shown, the mesh size (i.e., the size of each mesh element) is increased downward (z) and to the right (r) in sub-domains 2, 3, and 4. However, the mesh size decreases in the downward direction (z). This is done in an effort to capture the relevant thermal transport physics at the transducer/pillar interface (i.e., the thermal boundary conductance, G). This is particularly important when (1) we are sensitive to G and (2) we need to extract κ independent of G (as we do here). In this case, a graded mesh with 40 mesh elements in the radial direction and 20 mesh elements in the through-plane direction was used for each of the four solid domains that were modeled. We note that there is a less than 1% difference between the values of kz, G, and Cp obtained via the analytical and numerical solutions and that a mesh independence study was done to ensure that the solution was not dependent on the mesh itself. Prior to using the numerical simulation described here to fit any data, we use it to fit the data in Fig. 4.

The numerical model described in Sec. II is validated by fitting the measured data in Fig. 4 and comparing it with both the resulting analytically determined thermal properties and those available within the wider literature. The only changes to the models shown in Figs. 2 and 7 include (1) the removal of the pillar domain and (2) the transducer is split into two separate domains above the Si substrate [one from r = 0 to r=2×(w0+w1) and the other from r=2×(w0+w1) to r]. This effectively models the semi-infinite nature of the substrates shown in Fig. 4, which are blanket coated with 80 nm Au and 5 nm Ti. Numerical fits to the data are provided in the supplementary material and match well within 1% of the analytical fits provided in Fig. 4. Consequently, we consider the model used in this work to be valid for fits to the thermal properties of Si when the geometry is confined in the radial direction (e.g., Si micropillars).

The micropillars were fabricated using the following procedure (depicted schematically in Fig. 8). AZ5214E photoresist was spun on to the surface of a double-side polished silicon wafer [MTI Corporation, Si (100), undoped] with a spin rate of 4000 rpm, leading to a photoresist thickness of 1.7 μm, which was soft baked at 120 °C for 45 s. The patterned features were then written using a Heidelberg VPG200++ laser writer. The laser power was kept low enough for image reversal to successfully occur. The wafer was then baked at 110 °C for 60 s to cross-link the photoresist on the written features, and the entire wafer was flood exposed with a UV lamp with ample dosage to successfully complete the image reversal process. The wafer was then developed in a 4:1 diluted AZ400K photoresist developer for 30 s. The photoresist that was not exposed during the laser write, but that was exposed during the flood exposure washed away in the developer, whereas on the written features, the photoresist remained. The wafer was then cleaned and placed in a deep silicon etch tool (PlasmaTherm DSE), where the silicon was etched away to a depth of 29.5±1μm everywhere but where the photoresist remained. Then, the wafer was soaked in acetone to remove the photoresist from the written features, and the surface was cleaned in solvent rinse. Finally, a blanket-coating of 5 nm Ti and 80 nm Au was deposited by electron beam evaporation on the surface to serve as the FDTR transducer layer.

FIG. 8.

Schematic representation of the process flow to prepare the pillar arrayed investigated in this work.

FIG. 8.

Schematic representation of the process flow to prepare the pillar arrayed investigated in this work.

Close modal

SEM images of the micropillars used in this work are provided in Fig. 9. We note that while we only use a single pillar for characterization in this proceedings, an array of pillars was fabricated for use in future work. As shown, square pillars were also created but are not reported in this work.

FIG. 9.

Scanning electron microscope image of a pillar array used for this work. Note: only the three smallest circular pillars are shown, with diameters of (a) 20 μm, (b) 40 μm, and (c) 100 μm.

FIG. 9.

Scanning electron microscope image of a pillar array used for this work. Note: only the three smallest circular pillars are shown, with diameters of (a) 20 μm, (b) 40 μm, and (c) 100 μm.

Close modal

For this study, pillar C (rpillar=50μm) is used and measured with pump and probe radii of 16.8 μm and 12.8 μm, respectively. The ratio of pillar radius to pump radius (R=rpillar/w0) is, therefore, R3. The pillar height is also measured by SEM as h31μm (see supplementary material). Finally, electron beam (e-beam) evaporation is used to deposit an 80 nm Au transducer above the surface of the pillars for thermal characterization with FDTR.

In this section, we provide quantitative results that suggest the phase lag (i.e., temperature) response at the sample surface is substantially different for the confined geometry case (e.g., when r) across a wide range of frequencies. After we demonstrate this quantitative difference, we use our numerical fitting routine to determine the thermal conductivity of Si when in micropillar form in order to demonstrate the utility of the technique and distinguish it from the analytical solution. Consequently, we show that (1) the analytical solution is insufficient to capture the phase lag (and therefore the temperature response) at the sample surface for confined geometries and (2) a unique numerical fitting routine can be used to capture the physics that govern thermal transport in novel geometries, paving the way for the thermal characterization of micro- and nanoscale materials with thermal properties that are expected to depend on material geometry, such as ultra-drawn glass fibers and laser gain media.

We computationally model a Si pillar with a 20 μm height to highlight the expected effects that pillar diameter has on the phase lag at the sample surface. The pump and probe diameters we use are 7.55 μm and 1.55 μm, respectively, which is consistent with those used in Ref. 24 (and therefore provides secondary verification of the numerical results). The thermal boundary conductance is fixed at G=35MWm2K1 (consistent with Ref. 24), κ=144Wm1K1 and Cv=1.65MJm3K1 (both consistent with values reported in Ref. 14 for Si). To demonstrate the power of the numerical model, we plot the phase lag as a function of modulation frequency in Fig. 10 as rpillar is varied from rpillar=w0 to 6×w0.

FIG. 10.

(Left) Expected phase lag at the transducer surface as a function of modulation frequency for different values of the pillar radius (relative to the pump radius, w0). (Right) Cross-sectional view of the temperature rise within the sample at a modulation frequency of ω=4kHz (relative to the initial temperature with no applied power, T0) with an applied power of A0=50mW.

FIG. 10.

(Left) Expected phase lag at the transducer surface as a function of modulation frequency for different values of the pillar radius (relative to the pump radius, w0). (Right) Cross-sectional view of the temperature rise within the sample at a modulation frequency of ω=4kHz (relative to the initial temperature with no applied power, T0) with an applied power of A0=50mW.

Close modal

As shown in Fig. 10, the phase lag in the plot on the left is drastically altered as rpillarw0 and begins to converge with the analytical solution for a semi-infinite Si substrate (black dashed line) at high frequencies. This makes physical sense in that, at higher frequencies, the thermal penetration depth is reduced25 and therefore becomes less influenced by the confinement of heat at the pillar boundary. At low modulation frequencies (e.g., ω=4kHz, as used for the average temperature images in the right of Fig. 10), one can see that the influence of the boundary on the temperature distribution becomes more pronounced as rpillarw0.

The height of the pillar is also expected to alter the phase lag at the sample surface relative to the case of a semi-infinite substrate. Here, the expected phase lag is simulated with heights that range from hpillar=1μm to 50 μm. Figure 11 provides the relative impact of pillar height on phase lag across a range of pillar radii (rpillar=w0 to 4×w0).

FIG. 11.

Phase lag vs modulation frequency for (top left) hpillar=1μm, (top right) hpillar=2.5μm, (bottom left) hpillar=6.3μm, and (bottom right) hpillar=50μm.

FIG. 11.

Phase lag vs modulation frequency for (top left) hpillar=1μm, (top right) hpillar=2.5μm, (bottom left) hpillar=6.3μm, and (bottom right) hpillar=50μm.

Close modal

Figure 11 reveals that for comparatively low pillar heights and large pillar radii, the phase lag at the sample surface approaches what we would expect to see for measurements on bulk substrates. This is evident in the top left plot of Fig. 11, where pillar radii greater than 2×w0 sit neatly on the phase lag curve for semi-infinite Si (green dashed line). On the other hand, as the Si micropillar height becomes more pronounced, the phase lag at the sample surface begins to deviate drastically from the phase lag in the semi-infinite case. This is consistent with what we would expect physically; that is, as the pillar height becomes smaller, the semi-infinite substrate below it has an increasing impact on the temperature response at the sample surface. As before, a decreasing pillar radius also results in increasing deviation from the phase lag. Because these deviations are so apparent, we expect to be sensitive to the geometry of the pillar itself, permitting accurate experimental measurements of thermal transport properties via FDTR.

Here, we provide several measurements of the phase lag at the sample surface for Si micropillars that are between 28 and 33 μm tall and have varying radii. For this work, we utilize one of the circular micropillars shown in Fig. 9 (pillar C with rpillar=49.98μm and hpillar=29μm as measured via SEM). We note that we can vary the pump and probe radii used in our FDTR measurements by changing the objective lens magnification. This study utilizes a 20× objective to establish pump and probe radii of 16.8 μm and 12.8 μm, respectively.

As previously described, our numerical technique is used to fit the thermal properties of the pillar material below the transducer. The phase lag is fit to measured FDTR data with an applied pump power of A0=7mW, which produces a ΔT<1K at the transducer surface for our the lowest modulation frequency (1 kHz). As shown, the extracted thermal conductivity of the Si micropillar measured in this work is κ145Wm1K1, which remains consistent with the wider literature and within the measurement uncertainty for the conditions described in this work (Uκ=9.5%) relative to the measurements made on semi-infinite substrates.

We note here that an attempted fit of the data using the analytical model was not possible (shown in Fig. 12 is the expected phase lag when the sample is semi-infinite in order to clearly identify that the micropillar geometry has an impact on the phase lag at the transducer surface).

FIG. 12.

Numerical fit to the thermal properties of the Si micropillar and comparison to the expected phase lag for a completely semi-infinite geometry.

FIG. 12.

Numerical fit to the thermal properties of the Si micropillar and comparison to the expected phase lag for a completely semi-infinite geometry.

Close modal

A comparison between the obtained thermal properties for the semi-infinite geometry and the micropillar geometry is provided in Table I.

TABLE I.

Comparison of thermal properties for semi-infinite (SI) samples using analytical and numerical fitting routines to micropillar (MP) numerical fitting routine. G is not shown due to the lack of sufficient sensitivity.

PropertySI analyticalSI numericalMP numerical
κ (W m−1 K−1139 ± 7 139.2 ± 7.1 145 ± 10 
Cv (MJ m−3 K−11.5 ± 0.3 1.6 ± 0.2 1.59 ± 0.2 
PropertySI analyticalSI numericalMP numerical
κ (W m−1 K−1139 ± 7 139.2 ± 7.1 145 ± 10 
Cv (MJ m−3 K−11.5 ± 0.3 1.6 ± 0.2 1.59 ± 0.2 

Table I suggests that the formulation described here appropriately captures the thermal properties of materials whose geometries are confined and are quite different from those solutions which use the standard analytical fitting routine. Consequently, current analytical formulations11 for resolving the thermal conductivity of materials having non-semi-infinite geometries are insufficient to resolve their thermal properties.

We also note that there are significant alterations to the sensitivity of each of the parameters that we fit in our fitting routine, as shown in Fig. 13. Here, we show that large deviations in the curvature of our sensitivity distributions occur as a direct consequence of the geometry of the pillar. This could perhaps be used as an incentive for creating unique geometries that improve our sensitivity to individual thermal parameters in multilayer material systems whose thermal properties are otherwise difficult to extract (e.g., the thermal boundary conductance at SiO2/SiO2 interfaces).

FIG. 13.

Sensitivity to various thermal parameters as a function of modulation frequency for semi-infinite, planar Si (dashed lines) and 29 μm tall Si micropillars for the case in Fig. 12 (solid lines).

FIG. 13.

Sensitivity to various thermal parameters as a function of modulation frequency for semi-infinite, planar Si (dashed lines) and 29 μm tall Si micropillars for the case in Fig. 12 (solid lines).

Close modal

To provide further evidence that the fitting routine is able to capture the effects of the adiabatic boundary on the phase response at the upper surface of the pillar, we use our fitting routine to determine the thermal conductivity of Si micropillars having different diameters. We also utilize different pump sizes than those used in Fig. 12 (in this case, w0=21.2μm and w0=35.7μm). We first compare the impact of the relative size of the pump beams on the phase lag for micropillar diameters of 100 μm and 200 μm to highlight the sensitivity of the phase lag to the presence of confined boundaries, as shown in Fig. 14.

FIG. 14.

Phase lag vs frequency for micropillars with dpillar=100μm (blue) and 200 μm (red) and pump beam diameters of w0=21.2μm (circles) and 35.7 μm (diamonds). Solid and dashed lines represent the numerical fits to our data for the 21.2 μm and 35.7 μm diameter pump beams, respectively.

FIG. 14.

Phase lag vs frequency for micropillars with dpillar=100μm (blue) and 200 μm (red) and pump beam diameters of w0=21.2μm (circles) and 35.7 μm (diamonds). Solid and dashed lines represent the numerical fits to our data for the 21.2 μm and 35.7 μm diameter pump beams, respectively.

Close modal

In Fig. 14, the discrepancy between the phase lag for each of the two pump diameters is more pronounced for the pump beam with the larger diameter. This suggests that as the size of the pump beam is increased relative to the size of the pillar, the boundaries start to play a major role in the temperature at the surface of the transducer. We observe the same phenomenon when we hold the pump diameter constant and reduce the pillar cross section, as shown in Fig. 15.

FIG. 15.

Phase lag vs frequency for micropillars of varying size with w0=35.7μm.

FIG. 15.

Phase lag vs frequency for micropillars of varying size with w0=35.7μm.

Close modal

In Fig. 15, the phase lag is shown to be dependent on the ratio of the pump diameter to the diameter of the pillar. As the pump beam becomes more confined by the geometry of the pillar, its impact on the phase lag in the low frequency regime is more pronounced. We note that all fits resulted in a thermal conductivity of κ=145W/mK, which is consistent with that reported in the wider literature. This fitting routine is, therefore, critical when characterizing materials that are expected to have thermal properties that are dependent on finite geometries, such as ultra-drawn glass and polymer fibers26,27 and laser gain media, provided that care is taken to ensure that the pump beam is sized in such a way that we obtain sufficient sensitivity to the parameters of interest.

In this work, we develop a numerical model that is capable of fitting the phase lag of FDTR measurements on samples having non-semi-infinite geometries. We show that, for a known geometry, we can extract accurate values of thermal conductivity, volumetric heat capacity, and the thermal boundary conductance across the transducer/material interface for a non-typical geometry (in this case, a Si micropillar). The impact of the substrate is included to incorporate the impact of thermal spreading at low pump beam modulation frequencies, which we clearly demonstrate in a numerical parametric analysis. Given this new fitting paradigm, we expect thermal scientists to be able to characterize materials having geometry-dependent thermal properties, to include ultra-drawn glass fibers and laser gain media.

See the supplementary material for additional information on the numerical simulations, Matlab code for our numerical fitting routine for semi-infinite substrates, and our measurements and calculations of the pillar diameter and height.

R.J.W. and A.A.W. contributed equally to this work.

R.J.W. and B.F.D. would like to acknowledge primary financial support from Mr. Peter Morrison and the Office of Naval Research (ONR) (No. N001420WX00381). We also thank Dr. Mark Spector and ONR for equipment support (No. N0001420WX01170). N.M. thanks ONR for support (No. N000141612625).

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

1.
P.
Jiang
,
X.
Qian
, and
R.
Yang
,
J. Appl. Phys.
124
,
161103
(
2018
).
2.
R.
Yan
,
J. R.
Simpson
,
S.
Bertolazzi
,
J.
Brivio
,
M.
Watson
,
X.
Wu
,
A.
Kis
,
T.
Luo
,
A. R.
Hight Walker
, and
H. G.
Xing
,
ACS Nano
8
,
986
(
2014
).
3.
T.
Tong
and
A.
Majumdar
,
Rev. Sci. Instrum.
77
,
104902
(
2006
).
4.
K.
Kang
,
Y. K.
Koh
,
C.
Chiritescu
,
X.
Zheng
, and
D. G.
Cahill
,
Rev. Sci. Instrum.
79
,
114901
(
2008
).
5.
J. P.
Feser
,
J.
Liu
, and
D. G.
Cahill
,
Rev. Sci. Instrum.
85
,
104903
(
2014
).
6.
Z.
Ge
,
D. G.
Cahill
, and
P. V.
Braun
,
Phys. Rev. Lett.
96
,
186101
(
2006
).
7.
H.-K.
Lyeo
and
D. G.
Cahill
,
Phys. Rev. B
73
,
144301
(
2006
).
8.
M. D.
Losego
,
M. E.
Grady
,
N. R.
Sottos
,
D. G.
Cahill
, and
P. V.
Braun
,
Nat. Mater.
11
,
502
(
2012
).
9.
D. G.
Cahill
,
Rev. Sci. Instrum.
75
,
5119
(
2004
).
10.
J.
Zhu
,
D.
Tang
,
W.
Wang
,
K. W.
Holub
, and
R.
Yang
,
J. Appl. Phys.
108
,
094315
(
2010
).
11.
A. J.
Schmidt
,
R.
Cheaito
, and
M.
Chiesa
,
Rev. Sci. Instrum.
80
,
094901
(
2009
).
12.
T.
Favaloro
,
J.-H.
Bahk
, and
A.
Shakouri
,
Rev. Sci. Instrum.
86
,
024903
(
2015
).
13.
R.
Wilson
,
B. A.
Apgar
,
L. W.
Martin
, and
D. G.
Cahill
,
Opt. Express
20
,
28829
(
2012
).
14.
J. L.
Braun
,
D. H.
Olson
,
J. T.
Gaskins
, and
P. E.
Hopkins
,
Rev. Sci. Instrum.
90
,
024905
(
2019
).
15.
J.
Kimling
,
A.
Phillipi-Kobs
,
J.
Jacobsohn
,
H. P.
Oepen
, and
D. G.
Cahill
,
Phys. Rev. B
95
,
184305
(
2017
).
16.
X.
Chen
,
K.
Jarvis
,
S.
Sullivan
,
Y.
Li
,
J.
Zhou
, and
L.
Shi
,
Phys. Rev. B
95
,
144310
(
2017
).
17.
C. J.
Szwejkowski
,
A.
Giri
,
R.
Warzoha
,
B. F.
Donovan
,
B.
Kaehr
, and
P. E.
Hopkins
,
ACS Nano
11
,
1389
(
2017
).
18.
J. A.
Johnson
,
A.
Maznev
,
J.
Cuffe
,
J. K.
Eliason
,
A. J.
Minnich
,
T.
Kehoe
,
C. M. S.
Torres
,
G.
Chen
, and
K. A.
Nelson
,
Phys. Rev. Lett.
110
,
025901
(
2013
).
19.
T.
Kim
,
D.
Ding
,
J.-H.
Yim
,
Y.-D.
Jho
, and
A. J.
Minnich
,
APL Mater.
5
,
086105
(
2017
).
20.
A. J.
Minnich
,
Phys. Rev. Lett.
109
,
205901
(
2017
).
21.
R. J.
Warzoha
,
B. F.
Donovan
,
N. T.
Vu
,
J. G.
Champlaign
,
S.
Mack
, and
L. B.
Ruppalt
,
Appl. Phys. Lett.
115
,
023104
(
2019
).
22.
B. F.
Donovan
,
R. J.
Warzoha
,
R. B.
Venkatesh
,
N. T.
Vu
,
J.
Wallen
, and
D.
Lee
,
ACS Appl. Nano Mater.
2
,
6662
(
2019
).
23.
Z.
Wang
,
J. E.
Alaniz
,
W.
Jang
,
J. E.
Garay
, and
C.
Dames
,
Nano Lett.
11
,
2206
(
2011
).
24.
J.
Yang
,
E.
Ziade
,
C.
Maragliano
,
R.
Crowder
,
X.
Wang
,
M.
Stefancich
,
M.
Chiesa
,
A. K.
Swan
, and
A. J.
Schmidt
,
J. Appl. Phys.
116
,
023515
(
2014
).
25.
J. L.
Braun
and
P. E.
Hopkins
,
J. Appl. Phys.
121
,
175107
(
2017
).
26.
S.
Shen
,
A.
Henry
,
J.
Tong
,
R.
Zheng
, and
G.
Chen
,
Nat. Nanotechnol.
5
,
251
(
2010
).
27.
M.
Matsuo
,
C.
Sawatari
, and
T.
Ohhata
,
Macromolecules
21
,
1317
(
1988
).

Supplementary Material