A numerical fitting routine for frequency-domain thermoreflectance measurements of nanoscale material systems having arbitrary geometries

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.¹ In contrast to Raman spectroscopy² and 3-$ω$³ techniques, thermoreflectance-based measurements can separate the impacts of thermal boundary conductance (⁠$G$⁠) across interfaces and thermal conductivity (⁠$κ$⁠) within individual material layers⁴ 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.⁹ 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)¹⁰ and frequency-domain thermoreflectance (FDTR).¹¹ 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.⁹ 

Recent advances in thermoreflectance-based techniques include extensions to a steady-state system to gain sensitivity to the so-called “buried interfaces,”¹⁴ the use of a magneto-optical Kerr effect (MOKE) to gain sensitivity to interfaces having large thermal conductance¹⁵ and phonon-magnon coupling effects¹⁶ and the development of a transient thermo-transmission technique to measure the thermal boundary conductance of nanoparticles suspended in transparent media.¹⁷ 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.

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.¹¹ 

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.

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.¹¹ 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.

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

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

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

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

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

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

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

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

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$⁠,

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

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.

The thermal properties obtained above are consistent with those widely reported in the scientific literature.¹⁴ The thermal boundary conductance for all three samples is high relative to other values reported in the literature;¹⁴ 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

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

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,²³ 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.

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:

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 $−k∇T$ (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.

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.

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.

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, $R≈3$⁠. The pillar height is also measured by SEM as $h≈31μ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=35MWm−2K−1$ (consistent with Ref. 24), $κ=144Wm−1K−1$ and $Cv=1.65MJm−3K−1$ (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$⁠.

As shown in Fig. 10, the phase lag in the plot on the left is drastically altered as $rpillar→w0$ 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 reduced²⁵ 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 $rpillar→w0$⁠.

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$⁠).

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 $κ≈145Wm−1K−1$⁠, 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).

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

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 formulations¹¹ 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).

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.

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.

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 fibers^26,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.