The recent push for the “materials by design” paradigm requires synergistic integration of scalable computation, synthesis, and characterization. Among these, techniques for efficient measurement of thermal transport can be a bottleneck limiting the experimental database size, especially for diverse materials with a range of roughness, porosity, and anisotropy. Traditional contact thermal measurements have challenges with throughput and the lack of spatially resolvable property mapping, while non-contact pump-probe laser methods generally need mirror smooth sample surfaces and also require serial raster scanning to achieve property mapping. Here, we present structured illumination with thermal imaging (SI-TI), a new thermal characterization tool based on parallelized all-optical heating and thermometry. Experiments on representative dense and porous bulk materials as well as a 3D printed thermoelectric thick film (∼50 *μ*m) demonstrate that SI-TI (1) enables paralleled measurement of multiple regions and samples without raster scanning; (2) can dynamically adjust the heating pattern purely in software, to optimize the measurement sensitivity in different directions for anisotropic materials; and (3) can tolerate rough (∼3 *μ*m) and scratched sample surfaces. This work highlights a new avenue in adaptivity and throughput for thermal characterization of diverse materials.

## I. INTRODUCTION

The modern “materials by design” paradigm rests on the three pillars of computation, synthesis, and characterization, integrated in a high-throughput (H-T) manner.^{1,2} While the first two pillars have seen major advances with the development of high-performance computing and artificial intelligence,^{2–6} and autonomous robotic synthesis platforms,^{7–9} when thermal properties of materials are of interest, the H-T characterization step can be the key bottleneck. This is especially the case for the study of thermal conductivity (*k*),^{10,11} which is critical for many energy technologies from thermal management to thermoelectric energy conversion.^{12,13} To date, only a few hundreds of the ∼10^{5} laboratory-synthesized materials documented in the Inorganic Crystal Structure Database have *k* values experimentally measured, let alone manufactured materials in their real-world forms with defects, porosity, anisotropy, and/or complex morphology.^{7,12–14} This limited experimental dataset size poses an important challenge for *k* prediction with machine learning^{15,16} and further hinders the development of new thermal functional materials.^{11} Therefore, innovation in thermal metrology is strongly needed to complete the H–T materials research infrastructure.

An ideal H–T thermal metrology to facilitate materials screening should allow parallelized *k* measurements for property mapping and measuring multiple samples simultaneously, be able to resolve arbitrary thermal anisotropy (different *k* in different sample directions), require only simple sample preparation, and tolerate some level of sample roughness. As summarized in Table S1, these criteria are challenging to realize with current techniques. Contact-based thermal metrologies are often slow, need tedious sample preparation (e.g., steady state methods^{17,18}), and/or are not spatially resolvable (e.g., transient 3ω methods^{19,20} and transient plane source methods^{21,22}). On the other hand, non-contact optical methods are appealing because they require little or no microfabrication or handling of sensors. They also are attractive because of their potential to be integrated into autonomous robotic experimental platforms (e.g., 3D printing and combinatorial manufacturing),^{23,24} which fabricate batches of samples of small volume but similar geometry, for synergistic H–T synthesis and characterization. Traditional optical methods such as pump-probe laser-based techniques (e.g., time or frequency domain thermoreflectance: TDTR or FDTR)^{25,26} are capable of property mapping with micrometer-scale resolution and can measure thin film materials (<1 *μ*m)^{27–29} but suffer from the need for mirror smooth sample surfaces which can be difficult to prepare.^{30,31} Among other optical techniques, typical commercial laser flash method^{32} and general photothermal radiometry^{33} are convenient for bulk solid materials with uniform thickness and tolerate sample surface roughness to a certain extent^{34} though difficult for property mapping, often requiring an inconvenient angled configuration for the laser source and IR camera; and Raman scattering thermometry is well suited for 2D materials but is otherwise limited in application^{35,36} (see supplementary material Sec. 1 for more discussion).

In the pursuit of H-T measurement, although TDTR and 3ω have been modified to study a combinatorial materials library,^{37,38} these measurements are conducted in a series manner with samples or regions of interest (ROIs) scanned/studied one by one. In fact, current optical-based thermal metrologies, from laser flash, TDTR or FDTR,^{25,26,30} to probe beam defection^{39} and gas-microphone photoacoustic method,^{40} all rely on single-point measurement which is sometimes inefficient for large area raster scans. To the best of our knowledge, none of the current thermal techniques have demonstrated a parallelized measurement of more than one sample in a single shot which is ultimately important for the whole H–T experimental infrastructure to allow efficient batch screening of a large number of materials.

Another important challenge is measuring samples with anisotropic *k*, which is common in various materials of contemporary interest that have complex unit cells and/or microstructures. Extending existing techniques to anisotropic *k* measurements is nontrivial and often time-consuming, e.g., requiring microfabrication of extra heaters in electrothermal methods^{41} and reconfiguration of optics, beam scanning, or detector modification in optical methods.^{30,42}

To meet these key needs for parallelized, non-contact, measurements of multiple samples with anisotropic *k*, here we present the structured illumination with thermal imaging (SI-TI) method. This all-optical approach uses adaptive heating patterns with software-controllable shape and frequency, combined with large area temperature mapping by thermoreflectance imaging. We show how these software-reconfigurable heating patterns enable measurements of the *k* tensor of highly anisotropic materials, as verified on muscovite mica. We also demonstrate how SI-TI enables parallelized measurement of multiple samples for a batch of three bulk samples and an independent measurement of nine regions on a single sample, both in a single imaging field. This technique is also more tolerant to sample surface roughness than typical transient thermoreflectance methods, as shown in measurements of a rough frosted glass sample, a porous aerogel, a rough yttrium stabilized zirconia sample, and a 3D printed thermoelectric thick film (∼50 *μ*m). Therefore, SI-TI represents a new avenue to enhance the adaptivity and throughput of thermal characterization to facilitate experimental research on materials by design.

## II. EXPERIMENTAL

### A. Principles of SI-TI

The basic SI-TI thermal characterization platform is depicted in Fig. 1(a). It uses transient thermoreflectance (TR) imaging to map the evolution of sample surface temperature, $\Delta T(x,y,t)$, by monitoring changes in the sample's reflectivity $\Delta R(x,y,t)=dRdT\Delta T(x,y,t)$, where $dRdT$ is the TR coefficient.^{43} Note that $dRdT$ is nearly constant from 300 to 400 K for gold, which is used as the transducer layer in our experiments below.^{44} We note the distinction between optical techniques based fundamentally on single-point TR thermometry, such as thermal metrology using TDTR, FDTR, and photothermal/thermal wave microscopy, all of which rely on a monolithic photodetector,^{25,26,30,31,45,46} and techniques using large-field (kilopixel to megapixel) TR imaging. The latter has been used for temperature mapping of microelectronics,^{47,48} with some extension to the 1D heat transfer characterization of thin films and nanowires using Joule heating, or composites with a single-point-type laser pump-probe method.^{29,49–51} Here, we show how TR thermal imaging can be adapted for thermal property measurements with microscopic spatial resolution (see supplementary material Sec. S2.1) and overcome the fundamental limitation of single-point detection in traditional techniques, which require raster scanning for surface property mapping or measurement of multiple samples (more discussion in supplementary material Sec. 1).

Referring to the schematic of SI-TI in Fig. 1(a), the two main components of the system are as follows: (1) a digital micromirror device (DMD, Mightex Polygon400^{®}) coupled with an intensity modulated pump light to achieve spatial-temporal control of the optical heating at the sample plane, $PH(x,y,t)$; and (2) an infinity corrected optical microscope with a 2 Mpixel CMOS camera for TR imaging. DMDs are widely used in digital projectors, stereolithography,^{52} and structured illumination microscopy of biological systems^{53} but have not been used for thermal metrology to our knowledge. The key function here is that every pixel of the micromirror array can be independently controlled to either reject or reflect the pump light to the sample plane, thereby providing the “on” or “off” state of every pixel of $PH(x,y,t)$ (see Appendix and supplementary material Sec. S2.1 for more details).

The pump system [red geometries and dashed lines Fig. 1(a)] is based on a continuous wave LED at 470 nm wavelength (Mightex, ∼17 nm FWHM around the nominal 470 nm) with its intensity modulated temporally as square or triangular wave at a frequency *f*_{0} (see below). The pump light then passes through a liquid light guide, the DMD, a 475 nm short pass filter, a free space focusing lens system, a dichroic mirror with a cut-on wavelength of 490 nm, and a fluorite semi apochromat infinity corrected objective lens to form a sharply defined heating pattern on the sample surface (supplementary material Sec. S2.1). In most cases, a 20× magnification objective lens is used to ensure enough optical heating power density at the sample (∼40 mW mm^{−2}) for sufficient temperature rise, while at the same time a moderate numerical aperture (NA) of 0.4 is selected to ensure enough depth of field to tolerate the out-of-plane variations of sample surface morphology (see Sec. II D below). The corresponding field of view (FOV) is 572 × 357 *μ*m^{2}. In the demonstration of the spatially resolvability of SI-TI [Fig. 5(a)], we also used a 5× objective lens with NA = 0.1 and a FOV of 1427 × 2286 *μ*m^{2} which yields a smaller power density of 3.1 mW mm^{−2}.

Similar to typical LED or laser sources in microscopy, the pump light is not perfectly spatially uniform, but rather exhibits mild non-uniformities across the FOV. Such non-uniformity is calibrated and well approximated by a 2D quadratic function (Fig. S1, with typical RMS error <5%) and is taken into account in our model (see Subsection 3 of the Appendix and details in supplementary material Sec. S2 1 and S4). This calibration of the pump intensity distribution is a one-time step which does not have to be repeated from sample to sample [or pattern to pattern in certain cases (Fig. S4)].

We also have explored the effects of defocus in the case that the sample surface does not lie precisely in the focal plane. The fitting parameters for the pump distribution quadratic function vary by only a few percent when sample height is changed by up to ±10 *μ*m for the 20× lens, or ±40 *μ*m for the 5× lens, leading to no more than 1% change in the fit *k* result in the worst case. See also Figs. S2 and S3 and Table S2. For the applications, we envision for SI-TI such as thick film materials from batch synthesis [Figs. 4(c) and 4(g)], the height is typically uniform to within ±10 *μ*m.

A continuous wave 530 nm green LED [Thorlabs, ∼15 nm FWHM, green overlays in Fig. 1(a)] is used as the probe light to provide diffuse illumination for TR temperature mapping, which avoids the undesirable parasitic interference under high magnification if a more coherent laser source had been used.^{54} A total probe optical power of ∼1.9 mW mm^{−2} using the 20× lens (and 0.15 mW mm^{−2} using the 5× lens is obtained on the sample plane which is attenuated by a neutral density (ND) filter to avoid saturating the camera. The probe wavelength is chosen to maximize the TR coefficient of gold to enhance the signal-to-noise ratio (SNR, see supplementary material Sec. S2.5). A 530 nm longpass edge filter in front of the CMOS camera admits the green probe light while ensuring that the pump light leakage into the camera is negligible, <0.5% of the typical measured TR signal intensity. The TR signal carried by the probe light is collected by a 12-bit CMOS camera with a 1920 × 1200 resolution and a frame rate of 60 fps. With 20× magnification, the microscope imaging resolution is determined by the diffraction limit, ∼660 nm (supplementary material Sec. S2.1). The temperature resolution of the CMOS camera (supplementary material Sec. S2.3) can be characterized by the noise equivalent temperature difference (NETD), ∼18 mK (for 12 min averaging of a typical ROI of ∼10^{4} pixels), which is consistent with prior reports.^{55} We estimate a temperature sensitivity of 0.44 K Hz^{−1/2} based on the calibration shown in Fig. S5. Such a temperature resolution is enabled by dithering-based stochastic resonance in TR which provides count precision below the minimum quantization level of bit depth imaging (see supplementary material Sec. S2.2).^{56}

Samples for SI-TI should be flat but not necessarily mirror smooth, and as detailed below, we have successfully measured rough and/or porous samples including an aerogel, a scratched ceramic, and a 3D printed thin film (see Sec. IV for Fig. 4). For a consistently good optical absorption (*α *= 0.62 at 470 nm^{57}) and a high TR coefficient ($|dRdT|$ = 2.1 × 10^{−4}–2.36 × 10^{−4} K^{−1} at ∼530 nm^{44,58}), all samples are coated with a Au thin film (∼60–200 nm) transducer layer by e-beam evaporation (with 10 nm Ti adhesion layer, Subsection 4 of the Appendix). The thickness and *k* of the Au layer are obtained by a profilometer and a combination of four-probe electrical resistivity measurement and the Wiedemann–Franz law (see supplementary material Sec. S3). Heat conduction in the air domains is included in the thermal model (supplementary material Sec. S4). Convective and radiative heat losses are neglected because they are estimated to make <3% contribution to the results according to our numerical simulations (supplementary material Sec. S5.2). The current system is not designed for high temperature measurements, which would require a heating stage with an optically transparent window. Gold thin films tend to agglomerate at ∼600–650 °C,^{59} and at higher temperatures, metallic nitrides such as HfN^{60} and VN^{61} might be more appropriate.

### B. Data acquisition and analysis

Figure 1(b) depicts the concept of a typical SI-TI measurement for two samples, A and B. Two different rectangular heating patterns are used to enhance the measurement sensitivity to different components of the *k* tensor: for sample A, the heating pattern is elongated along the *x* direction, so that the dominant direction of transient heat flow is along the *y* and *z* directions; conversely, sample B uses a heating pattern elongated in the *y* direction, so that the heat diffusion is more sensitive to the *x* (and *z*) components of *k*. Periodic optical heating by the pump light at frequency *f*_{0} creates a thermal wave in the sample, which is monitored by the CMOS camera through lock-in thermography, using a series of 210 images captured at 60 fps over a time window of *t*_{p} = 3.5 s. This data acquisition process is then repeated and averaged over tens to hundreds of cycles of such windows to further improve the signal-to-noise ratio (SNR) (see Subsection 2 of the Appendix for details).^{62} The processed TR image, $\Delta Rx,y,tR0,$ where $R0$ is the average reflectivity of the Au transducer at room temperature, is related to the sample surface temperature response $\Delta Tx,y,t$ by the normalized TR coefficient $Cth=dRR0dT$ (see Subsection 2 of the Appendix). The obtained $\Delta Tx,y,t$ is averaged spatially over selected regions of interest (ROI) $Ri$ to obtain $\Delta Tit$ and further Fourier transformed to the frequency domain for analysis (see below and supplementary material S2 4). For the 530 nm probe LED, our calibration yields $Cth=\u22122.57\u2009\xb1\u20090.24\xd710\u22124\u2009\u2009K\u22121$, which is consistent with a literature value of $2.5\xd710\u22124\u2009K\u22121$ (see supplementary material Sec. S2.2).^{62} With the small (<5 K) temperature changes in our experiment and the temperature-stable $dRdT$ for gold, $Cth$ can be reasonably treated as a constant here. Notably, by analyzing the absolute phase lags of the regions with respect to the heating wave, or the ratios of the frequency domain temperature responses among different regions, several experimental parameters cancel out (namely, *α*$,\u2009Cth,$ and the heating power). Thus, typical SI-TI measurements do not need to calibrate these parameters in most cases, as will be shown later.

To obtain information about multiple frequencies from a single trial for higher measurement throughput, rather than sine-wave modulation of the pump, we use square or triangle wave profiles because they automatically contain higher order harmonics at odd multiples of the heating frequency *f*_{0} (e.g., *f*_{0}, 3*f*_{0}, and 5*f*_{0}). 0.286 < *f*_{0} < 4.000 Hz is used (hence, the maximum of 5*f*_{0} = 20.000 Hz, as exemplified in supplementary material S8.4) to ensure a sufficiently high temperature rise of the sample to TR imaging, which is also well below the Nyquist frequency limit of 30 Hz set by the camera frame rate. Although a higher *f*_{0} can be achieved by methods such as heterodyne lock-in thermography^{47} or compressed sensing^{63,64} to surpass the Nyquist limit, it is not necessary for our experiment since the current range is sufficient to provide measurement sensitivity (see Fig. 3 and supplementary material Sec. S6), and undesirable drops of $\Delta T$ at high *f*_{0} undermine the SNR. Because the current SI-TI apparatus has relatively small optical power density, the demonstrations below focus on materials with relatively low *k* (<5 W m^{−1} K^{−1}) to ensure a reasonable temperature rise in the sample (see Sec. IV and supplementary material S1 and S2 5).

To analyze the data, the spatially averaged $\Delta Ti(t)$ of selected ROIs are converted to $\Delta Ti(f)$ in the frequency domain using the fast Fourier transform (FFT). To avoid spectral leakage and scalloping losses of the signal, *f*_{0} is chosen to be *n*/*t*_{p}, where *n* is an integer such that the harmonic signals all fall into FFT bins exactly (see supplementary material Sec. S2.4). With temporal averaging for each frame, the general SNRs of the SI-TI measurement, defined as the peak to average noise floor ratio in the frequency domain, with the 20× lens are >50 and >15 for the first and third harmonic signals, respectively (the minimum SNR is for the ROI farthest from the heater for the sample with the highest *k*; see supplementary material Sec. S2.4). The SNRs of the 5× lens measurement of glass in Fig. 5 (direct heated region) are >30 and >10 for the first and third harmonic signals, respectively.

The resulting frequency domain data, $\Delta Ti(fj)$, where $fj$ are discrete frequencies being odd multiples of *f*_{0}, are analyzed by an analytical model developed for layered geometries^{25,65,66} generalized for the rectangular-shaped heater and accounting for spatial non-uniformity of the pump heating (see Subsection 3 of the Appendix and details in supplementary material Secs. S2.1 and S4). The analytical model was substituted by a COMSOL Multiphysics^{®} finite element modeling (FEM) analysis for the paralleled measurement of three samples [Figs. 5(b) and 5(d)] to account for the sample lateral boundaries. In nonlinear least squares regression, the sample(s) *k* is derived by fitting the experimental data with *k* as the only adjustable parameter. The sample heat capacity *C* is fixed at an accepted literature value or determined from separate differential scanning calorimetry (DSC) measurements (supplementary material Sec. S3.5 and Table S3).

Two approaches are used for the uncertainty estimation. For most analysis, we used the traditional method considering error propagation from all fixed parameters in the analytical model and experimental noise approximated from the variance of repeated measurements following the method in Ref. 67, generalized here to complex data (see supplementary material Sec. S7.1). To understand the error propagation through model parameters and the effects of reducing the total quantity of data on the uncertainty, a Markov Chain Monte Carlo (MCMC) method was used under the framework of Bayesian parametric inference, where the uncertainties and mean values of all parameters are assigned in the corresponding priors (see supplementary material Sec. S7.3). This MCMC method was applied to the glass and the yttrium stabilized zirconia sample. While our primary technique with least squares fitting uses measurements at many different $fj$ [e.g., data in Figs. 2(e) and 2(f) and 3(d), 3(h), and 3(l) below], we also have demonstrated the use of only a single heating frequency with the MCMC method which can enable a higher throughput at the cost of larger uncertainty in the fit *k* (supplementary material Sec. S7.3).

## III. RESULTS

### A. Validation of SI-TI with standard silica glass

We now present several experiments to validate SI-TI for thermal property characterization. We begin with a material with isotropic *k*, namely, a silica glass microscope slide. The sample surface is cleaned with ethanol and coated with a 184-nm-thick gold film. The sample is heated using a simple rectangular-shaped pattern [Figs. 2(a) and 2(b)] which is driven by a triangular waveform in time. A square wave time-forcing (used for most other samples, Fig. S31) was also tested for three randomly selected frequencies and gave consistent results, which are omitted here for brevity. Driven by this heating, Fig. 2(c) shows a representative thermal image for a particular heating frequency and time delay and using the literature TR coefficient *C*_{th0} = −2.5 × 10^{−4} K^{−1}.^{62} The TR temperature mapping is also validated by a hybrid electrical heating TR thermometry measurement of the same type of glass sample in supplementary material Sec. S8.1 (Fig. S28). The heat diffusion beyond the original heating pattern (R1) is clearly visible. Notably, the partial vignetting (reduced image intensity) seen near the left edge of Fig. 2(b) is almost completely gone in the thermal image of Fig. 2(c). This robust feature of SI-TI arises from the TR imaging approach with pixel-by-pixel normalization (see Subsection 2 of the Appendix and supplementary material Sec. S2) which makes the probe intensity irrelevant to leading order, at least until the TR signal becomes so low that SNR issues become a concern (see SNR discussion in supplementary material Sec. S2 5).

To extract thermal properties from the raw TI data, we next calculate the spatially averaged temperatures, $\Delta Tit|f0$, at time *t* and for various triangle-wave fundamental heating frequencies $f0$, of the five ROIs [*i* = R1–R5, marked in Fig. 2(c)]. For one representative frequency, the resulting $\Delta TR1$, $\Delta TR3$, and $\Delta TR5$ are shown by the black circles in Fig. 2(d), where $\Delta TR2$ and $\Delta TR4$ are omitted for visual simplicity. Note that the absolute phase of the thermal wave is not important as phase differences between ROIs are used in the subsequent analysis. Then, for every heating frequency $f0$, the $\Delta Tit|f0$ time series is Fourier transformed into the frequency domain, yielding a complex $\Delta Tif$ with amplitude $Ai(f)$ and phase $\Phi i(f)$ (supplementary material Sec. S2 4). Only signals at *f *=* f*_{0} and 3*f*_{0} corresponding to the first and third harmonic of the Fourier series of the heating waveform are retained for analysis. Next, to facilitate a calibration-free measurement (eliminating the need to measure $Cth$ and the optical absorptivity of the transducer layer) and to minimize any phase uncertainties from the electrical input and signal delay, we use a relative approach which refers all results to R1. Specifically, in Fig. 2(e), we plot the amplitude ratios of R2,…,R5 each with respect to R1, denoted $Ai/1\u2261ARi/AR1$, and similarly in Fig. 2(f), the phase differences $\Delta \Phi i\u22121\u2261\Phi Ri\u2212\Phi R1$. Since the forcing function is a triangle rather than sine wave, we are able to extract and plot meaningful FFT information for both the first and third harmonic components, depicted with black and blue points, respectively.

Finally, the frequency domain data corresponding to Figs. 2(e) and 2(f) are fit with a one-parameter thermal model (see Subsection 3 of the Appendix and supplementary material Sec. S4), which yields *k*_{glass} = 1.40 ± 0.12 W m^{−1} K^{−1} (*C*_{glass} is fixed at 2.12 MJ m^{−3} K^{−1} as determined by a separate DSC measurement; see supplementary material Sec. S3 5). As shown in Figs. 2(d)–2(f), agreement between the model (lines) and the experiment (points) is excellent. This fit value of *k*_{glass} agrees very well (within 3%) with our independent *k* measurements of this type of glass sample using both a 3ω method and the above-mentioned hybrid electrical-heating-optical-thermometry technique (supplementary material Sec. S8 1); and these values are also all within the range of literature reports, 1.3–1.4 W m^{−1} K^{−1}.^{68,69} The estimated uncertainty of the fit *k*_{glass} here is obtained from the traditional error propagation mentioned above (supplementary material S7 1). To further understand the error propagation from model parameters and experimental data size, a Bayesian inference based on MCMC was also utilized to analyze the measurement of another, nominally identical glass sample which yields a posterior median with 95% (±2σ) credible interval of $kglass\xb12\sigma =1.36\u22120.180.19$ W m^{−1} K^{−1} (see supplementary material Sec. S7 3 1 and Table S4). This MCMC method is also used to investigate the effect of fitting using data from a single *f*_{0} (see Sec. IV). The excellent agreement between SI-TI measured and independently measured values of *k*_{glass} is also depicted below in Fig. 6, which summarizes all *k* measurements from this study.

### B. Measurement of anisotropic materials using adaptive heating patterns

We further demonstrate the unique advantage of SI-TI's adaptive heating patterns by independently determining all three principal thermal conductivities for the *k* tensor of a thermally anisotropic material, muscovite mica (grade-1 natural mica from SPI^{®}). Mica's layered structure^{70,71} leads to a relatively large $kmica,xx=kmica,yy$ in the in-plane directions^{72} but a small cross plane $kmica,zz$ with a thermal anisotropy ratio of ∼9. To effectively measure the *k* tensor, we need to create three different heating patterns [Figs. 3(a), 3(e), and 3(i)], each optimized to maximize the measurement sensitivities [Figs. 3(c), 3(g), and 3(k)] to a different component of the thermal conductivity tensor: $kmica,xx$, $kmica,yy$, and $kmica,zz$, respectively (see also Figs. S10 and S14 for other parameters exemplified with the case of glass). The sensitivity of the SI-TI measurement to a parameter $\theta $ in the thermal model is defined as $S\theta Y=\u2202ln(Y)\u2202ln(\theta )$, where *Y* is the signal. Since the phase is dimensionless and invariant with the addition of an integer multiple of 360°, no logarithm in the numerator is taken when *Y* = phase difference. For the unknown parameters being fit for, e.g., $\theta =kmica,yy$ of the sample, a large $S\theta Y$ is preferred to minimize the uncertainty in the fit value of *θ*. We emphasize that the novelty of the SI-TI method here is that these three heating patterns can be easily and rapidly changed purely using *software*, without any manual reconfiguration of the optics or moving the laser as in other photothermal methods^{30,73} or microfabricating multiple electrical heaters and thermometers as in electrothermal methods.^{41}

Square waves at seven different *f*_{0} are used to optically heat the thin mica sheet (274 *μ*m thick) sample, using the three heater patterns depicted in Figs. 3(a), 3(e), and 3(i). For each heating pattern, the mica *k* measurements in Fig. 3 follow a similar procedure as the glass measurements of Fig. 2. Representative thermoreflectance snapshots for a specific frequency are given in Figs. 3(b), 3(f), and 3(j), which clearly demonstrate the effect of the three different heating patterns. The utility of the different heating shapes is further quantified in the sensitivity analyses of Figs. 3(c), 3(g), and 3(k), which are consistent with qualitative intuition. For example, for pattern 1, Fig. 3(c) shows that the temperature rise of the top region R5 is mainly determined by the in-plane heat flow upward (*y* direction) and hence is most sensitive to $kmica,yy$, since the calculations show $|Skyy\Delta \Phi 5\u22121|>|Skxx\Delta \Phi 5\u22121|$, $|Skzz\Delta \Phi 5\u22121|$ for all measured frequencies. Similarly, Fig. 3(g) quantifies how measurements using pattern 2 are most sensitive to $kmica,xx$. For pattern 3, we use the absolute temperature rise rather than phase in order to maximize the sensitivity to $kmica,zz$ [Fig. 3(k)], since the phase difference between these two ROIs here is small and also not sensitive to $kmica,zz$. Figures 3(c), 3(g), and 3(k) also show that a single type of heating configuration, which is typical in previous thermal imaging-based methods,^{29,49,50,74} cannot provide effective measurements of all 3 *k* components (e.g., pattern 1 has almost no sensitivity to $kmica,xx$).

Finally, the frequency-domain results of Figs. 3(d), 3(h), and 3(l) are used to fit the three components of *k* using an iterative approach (see supplementary material Sec. S8.2 with the corresponding flow chart). The inputs to the fit are the complex temperature ratio from Patterns 1 and 2 [with phase data in Figs. 3(d) and 3(h)], as well as the absolute temperature rises amplitude of all three patterns (those for patterns 1 and 2 are omitted from Fig. 3 for brevity). This is a four-parameter fit yielding $kmica,\u2009xx$, $kmica,yy,\u2009kmica,zz,$ and a lumped optical constant (the product of *α*, *C*_{th}, and the incident heating power) with the ratio data in Patterns 1 and 2 providing $kmica,xx$, and $kmica,yy$, and the absolute value data yielding the rest. We emphasize that $kmica,xx$ and $kmica,yy$ are fit independently to demonstrate this capability of SI-TI; namely, we do not enforce $kmica,xx=kmica,yy$ even though this is known for mica based on its crystal symmetry.

The resulting converged fit values for the anisotropic thermal conductivity tensor of this muscovite mica are *k _{mica,}*

_{xx}= 4.04 ± 0.39 W m

^{−1}K

^{−1},

*k*

_{mica,}_{yy}= 4.02 ± 0.48 W m

^{−1}K

^{−1}, and

*k*

_{mica,}_{zz}= 0.46 ± 0.066 W m

^{−1}K

^{−1}. These values are in excellent agreement (better than 5%) with a literature report of 4.05 and 0.44 W m

^{−1}K

^{−1}for in- and cross-plane thermal conductivity in the same material,

^{72}as summarized in Fig. 6.

A detailed sensitivity analysis of different combinations of heater and ROI geometries and heating frequencies is given in supplementary material Sec. 6.2 (Fig. S17), which shows that the configurations used in our experiment are reasonably close to optimal for the FOV size and heating frequency range in the current setup. The dependence of the sensitivity on the aspect ratio or size of the heaters and ROIs is generally not strong except for some extreme cases. A relatively high-frequency measurement combined with large ROI shift using either amplitude ratio or phase is most effective when fitting for a single in-plane *k* (i.e., *k _{xx}* or

*k*), while the amplitude data are most helpful for extracting the cross plane

_{yy}*k*

_{zz}. Within the range of our experimental dimensions and frequencies, it is in fact not necessary to use a large heater as in Fig. 3(i) for

*k*

_{zz}since the amplitude sensitivity to

*k*

_{zz}depends only weakly on the heater/ROI geometry and frequency. The

*xy*spatial resolution of

*k*

_{zz}mapping is largely determined by the in-plane thermal diffusion (penetration) length of the measurement, defined as $dpi=\alpha si\pi f0$, where $\alpha si$ is the sample in-plane thermal diffusivity (see Sec. IV below).

### C. Demonstration of SI-TI for challenging regimes: Low-*k*, porous, and rough materials

Having established and validated the basic capabilities of the SI-TI technique above, we now demonstrate its extension into several areas, which can be challenging for traditional methods. First, we consider aerogels and similar porous low-*k* materials, which are challenging for traditional steady state or transient methods in which the parasitic heat losses can lead to large errors. As shown in the SEM image in Fig. 4(d), we used SI-TI to measure *k* of an 11-mm-thick sample of commercial aerogel (Airloy^{®} x56), which we have measured previously using hot disk and guarded hot plate methods [see Fig. S30(a)–S30(d)].^{22} With SI-TI, thanks to the small thermal diffusion length in the aerogel (tens of micrometers for *f*_{0} > 1.4 Hz used here, see also Sec. IV B), the thermally probed zone is relatively small which minimizes parasitic heat losses due to radiation and convection, estimated as <3% in the measurement here (see supplementary material Sec. S5 2). Undesirable heat loss to the sensor encapsulation in contact methods (e.g., hot disk)^{22} is also eliminated.

This aerogel sample has very large porosity and significant surface roughness, ∼93 nm, both of which would be challenging for measurements by traditional TDTR and FDTR, which require a more continuous, specularly polished surface finish. A representative time domain dataset for this sample is shown in Fig. S31(a), which demonstrates high SNR. For SI-TI, we found the analysis is more robust if we first remove dark spots and defect pixels from the raw image by applying an intensity high pass filter, set to eliminate pixels with intensity <10% of the maximum intensity of the image (this typically eliminates ∼6% of the pixels caused by adsorbed particles or defects in an ROI). Combined with a low measurement frequency to probe large area, this rough aerogel can still be analyzed reasonably with SI-TI (see next paragraph and supplementary material S8.3 for the analysis of YSZ with even higher roughness). The thermal model used to analyze the raw data are bidirectional, accounting for heat conduction into both the sample and the static air above it (see supplementary material Sec. S4 for the model and S8 3 for results). Least square fitting of the frequency domain data gives *k*_{aerogel} = 0.020 ± 0.0034 W m^{−1} K^{−1}, which agrees with our previous measurements of the same material using a guarded hot plate method to within 7%.^{22} This agreement is also depicted graphically in Fig. 6, and more details are given in supplementary material Sec. 8 3.

To further demonstrate the ability of SI-TI to tolerate sample roughness, we next measured an imperfectly polished piece of nanocrystalline yttrium stabilized zirconia: 8 mol. % Y_{2}O_{3} (YSZ) and a highly pitted frosted glass sample [Figs. 4(b), 4(e), and 4(f)] both with much higher roughness than the aerogel sample. As seen in Fig. 4(e), the YSZ sample has an average surface roughness of 0.85 *μ*m which can causes a great deal of light scattering and non-uniform optical contrast, making it very challenging for measurements using traditional laser-based TR methods (TDTR, FDTR) and probe beam deflection methods.^{30,39,75} The sample was prepared by the well-established Current Activated Pressure Assisted Densification (CAPAD) method^{76} and is ∼1.1 mm thick (see Subsection 6 of the Appendix). Microstructure analysis shows that the average grain size is 184 ± 48 nm (mean ± std. dev; see Fig. S7). We used a relatively thick transducer layer (210 nm Au) for better coverage of some of the roughness length scales. Similar to the aerogel image processing, we apply an intensity high pass filter to the raw TR images of YSZ with a high pass threshold of ∼23% of the maximum intensity in the ROI, which eliminates ∼10% of the pixels generally associated with the particle-like defects and black border regions (supplementary material Sec. S8.3 and Fig. S32). In other words, regions with poor transducer conformity that appear black are excluded from the data processing. Despite the visually apparent contrast in Figs. 4(e) and 4(f), the spatial variations of the absorption of the Au transducer may be less since both diffuse and specular reflections may correspond to similar absorptivity yet look very different in these near-normal-incidence reflection maps.

Further investigation of the influence of image processing and image filtering on the TR map is performed by removing a random fraction of pixels from the ROI and systematically varying the filter threshold as shown in supplementary material Sec. S8 3 and Figs. S33 and S34 for the Bi_{2}Te_{2.7}Se_{0.3} sample. These studies indicate that for rough samples with relatively uniform morphology, the SNR is mostly degraded by pixels with low intensity. In addition, randomly removing even 90% of the pixels was found to cause only a small impact on the extracted amplitude and phase. This indicates that removing only several tens of percent of the pixels, as done in the actual image processing for the main experiments, does not cause significant artifacts in the measurement. We emphasize that such pixel removal in postprocessing impacts the thermoreflectance thermometry, i.e., the probe, but has no impact on the heat flux which arises from the pump.

Another key point that allows SI-TI to tolerate rough samples is that its cross-plane thermal penetration length $dpc=\alpha sc\pi f0$, where $\alpha sc$ is the sample cross-plane thermal diffusivity, is 2–3 orders of magnitude longer than $dpc$ in TDTR and FDTR, since SI-TI's characteristic *f*_{0} is ∼4–6 orders of magnitude smaller (on the other hand, SI-TI's large in-plane $dpi$ negatively impacts its spatial resolution, as discussed later). The large $dpc$ in SI-TI means that it can much better retain the assumption of a homogeneous flat sample domain in the thermal model, which is expected to hold as long as the roughness scales are small compared to $dpc$. Thus, with its large area averaging, image post processing, and relatively long thermal diffusion length, SI-TI is more tolerant to sample surface roughness as compared to other laser-based TR methods. Representative time domain datasets of YSZ and frosted glass are shown in Figs. S31(b) and S31(c).

Finally, our measurements of the YSZ sample of Fig. 4(e) yield *k*_{YSZ} = 2.06 ± 0.18 W m^{−1} K^{−1}, in good agreement with *k*_{YSZ} = 2.0–2.5 W m^{−1} K^{−1} reported previously^{77} for dense nanocrystalline YSZ fabricated using a similar method and with a similar grain size. To further confirm the reliability of the measurement, we polished the top surface of the same YSZ sample and repeated the SI-TI study using an intentionally different heating pattern (supplementary material Sec. S7 3 2). The measurement on the smooth surface with roughness of 35 nm yields *k*_{YSZ} = 1.96 ± 0.14 W m^{−1} K^{−1} from least squares fitting and $kYSZ\xb12\sigma =1.95\u22120.220.22$ W m^{−1} K^{−1} from the MCMC Bayesian inference (95% credible interval) which confirms the direct SI-TI measurement on the rough surface (see Sec. IV and supplementary material Sec. S7 3.2 for Bayesian inference analysis).

To explore the maximum sample roughness measurable by SI-TI, a frosted glass sample is studied as shown in Fig. 4(f). This sample has a pitted surface with high roughness of 3.4 ± 0.9 *μ*m [see supplementary material Fig. S30(m)]. Note that this average roughness also happens to be close to the microscope depth of field under the 20× lens for the 530 nm probe light, which is ∼3.3 *μ*m estimated based on Ref. 78. Thus, some portions of the image can appear slightly out of focus. The acquired SI-TI temperature data from the heated region still shows a reasonable signal-to-noise ratio and we were able to obtain meaningful thermal conductivity data, while in contrast the thermoreflectance response for shifted ROIs was too noisy to be readily used. Meanwhile, the absolute amplitude is inconvenient with the surface absorption expected to be different from the smooth Au/glass samples. Thus, we only used the absolute phase lag with respect to the heating waveform of the heated area for analysis [Fig. S30(o)]. By fitting the phase data extracted from two areas (∼5 mm away from each other) which are found to be consistent with each other, we obtained $kfrosted\u2009glass$ = 1.48 ± 0.09 W m^{−1} K^{−1} which is consistent with the measurement of the smooth region of this glass slide [1.42 W m^{−1} K^{−1} for Fig. 5(a)] and typical values of silica glass. These results demonstrate that the current SI-TI can in some circumstances tolerate roughness of at least ∼3.3 *μ*m. We speculate this roughness limit is mainly determined by the microscope depth of field, in which case it could be increased by using a lower NA lens with more forgiving depth of field, but we have not explored this experimentally. We note also that photothermal radiometry techniques can tolerate even higher sample roughness though typically require mm to cm scale sample size.^{79}

Moving beyond bulk samples, another attractive application of SI-TI is for screening studies of emerging materials often with small quantities and dimensions prepared by various H–T techniques like additive manufacturing and combinatorial synthesis. As one example, historically the thermal conductivity measurement of printed thermoelectric thick films (tens of *μ*m) has been challenging since bulk techniques are not applicable and due to the large surface roughness. As a proof of concept to overcome such challenges, we have used SI-TI to measure *k* of an extrusion printed thermoelectric Bi_{2}Te_{2.7}Se_{0.3} (BTS) thick film (see Subsection 5 of the Appendix for synthesis details), shown in Fig. 4(g). This 49-*μ*m-thick film is composed of Bi_{2}Te_{2.7}Se_{0.3} nanoplates with lateral dimensions of approximately one *μ*m and random orientations (Figs. S8 and Fig. S30) and is being pursued for flexible thermoelectric energy harvesting devices.^{24,80,81} Note that SI-TI would not be suitable for films much thinner than the minimum cross-plane thermal penetration length $dpc$, >135 *μ*m in this case. Figure 4(g) shows a micrograph of the film after printing on a muscovite mica substrate (the same type as studied in Fig. 3) and coating with a 99-nm-thick gold transducer layer, and also shows a representative heating pattern. The average surface roughness is ∼0.31 *μ*m. The inset shows a macroscopic photo of the sample.

As detailed in supplementary material Fig. S30, fitting the SI-TI data gives *k*_{BTS} = 0.65 ± 0.07 W m^{−1} K^{−1}, which is consistent with a literature report of 0.7 ± 0.1 W m^{−1} K^{−1} for BTS synthesized with a similar method.^{82} The sensitivity analysis of this thick film sample is shown in supplementary material Sec. 6.1 (Figs. S15 and S16). This demonstration highlights how SI-TI can be an effective tool for characterizing rough film materials on substrates, which are hard to measure with bulk steady state methods, TDTR/FDTR, or the laser flash method (which needs uniform thickness and centimeter scale lateral dimensions). We note that the current $dpc$ in SI-TI is on the order of 100 *μ*m which is much larger than TDTR/FDTR methods which can provide $dpc$<1 *μ*m. This currently limits SI-TI to films of thickness ≫1 *μ*m. This restriction originates from practical factors including the relatively small optical power and camera frame rate (60 fps) which are not fundamental limitations (see Sec. IV B). As a non-contact optical method that requires relatively simple sample preparation and small sample volume, it is expected that SI-TI can function synergistically with 3D printing to enhance the efficiency in discovery and deployment of new thermoelectric materials.

### D. Multiplex and parallel measurement in a single field of view

To demonstrate the multiplexing and spatial resolvability of SI-TI, we performed a measurement on a silica glass sample [smooth part of the same sample as Fig. 4(f)] using nine separate but simultaneous pump heating patterns to individually study the local *k* of each region. See Figs. 5(a) and 5(c). In this measurement, the ROIs are selected to overlap with their respective heater patterns. The SNR exceeds 10 for the third harmonic signal [Fig. S35(a)]. To ensure the locality of each thermal measurement, a 5× lens is used which gives a 4 times larger FOV (16× larger by area) to allow enough distance [“blank streets” in Fig. 5(c)] between adjacent heating patterns. Frequencies above 3 Hz are used to ensure the in-plane thermal penetration lengths ($dpi<265\u2009\mu m$) are similar to or smaller than this inter-heater gap (∼230–280 *μ*m). By comparing the model result for the full model including all nine heaters and a hypothetical case with only a single heater, it is found that once *f*_{0} > 3 Hz, the amplitude of the temperature responses for each ROI from the two cases overlap. This indicates that our analysis of each ROI using amplitude data above 3 Hz is local enough to be considered independent of each other. Thus, all 9 ROIs are measured simultaneously and independently, using the amplitude data. All nine results are mutually consistent and yield a mean and standard deviation of 1.42 ± 0.08 W m^{−1} K^{−1} as shown in Fig. S35 and summarized in Table S5. See more details in supplementary material Sec. S8 4. Similar to other transient optical thermal metrologies,^{29} the spatial resolution of such *k* mapping measurement is determined by the thermal diffusion length. As such, the current thermal resolution of SI-TI is largely limited by the low heating frequency (see Sec. IV B).

Finally, we highlight the capability of SI-TI to simultaneously measure multiple samples in a single FOV. To the best of our knowledge, such parallelized thermal transport measurements have not been demonstrated in any prior techniques, which have instead relied on sequential measurements and serial raster scanning.^{29,37,38} We demonstrate this capability using three distinct samples placed within a single FOV, with small gaps as shown in Fig. 5(d), which in principle could be part of a large batch of similar samples depicted in Fig. 5(b). The samples are low-*k* single crystals with 0.5 mm thickness and 5 × 5 mm^{2} lateral dimension procured from MSE Supplies LLC, specifically β-PbF_{2}, Bi_{4}Ge_{3}O_{12} (abbreviated BGO), and Pb(Mg_{1/3}Nb_{2/3})O_{3}-PbTiO_{3} (with 28–30 mol. % PbTiO_{3}, abbreviated PMN-PT). As shown in Fig. 5(d), each sample is irradiated with its own rectangular-shaped optical pump, simultaneously. Square wave modulation is used in the time domain for all three samples. To analyze the raw data, we used COMSOL FEM software to more accurately account for the finite sample extents in-plane, i.e., the lateral sample boundaries seen as dark regions in Fig. 5(d). See the detailed data fitting results in supplementary material Sec. S8 4 (Fig. S36) and sensitivity analysis in Fig. S18.

BGO and β-PbF_{2} have cubic lattice structures and hence isotropic *k*. PMN-PT has weak crystal anisotropy and a single value of *k* is typically reported in the literature,^{83,84} and thus, we also treat it as isotropic in the FEM simulations. As summarized in Fig. 6, the FEM-based fitting results agree well with prior measurements in literature: 1.44 ± 0.12 W m^{−1} K^{−1} measured here for PMN-PT (compared to 1.38 W m^{−1} K^{−1} in Refs. 83 and 84), 2.80 ± 0.31 for BGO W m^{−1} K^{−1} (compared to 2.6 W m^{−1} K^{−1} in Ref. 86), and 1.35 ± 0.15 W m^{−1} K^{−1} for PbF_{2} (compared to 1.40 W m^{−1} K^{−1} in Ref. 85). As further verification, we also conducted separate SI-TI measurements of these three samples individually, placing the heater pattern in the center region of each sample with ∼5 mm lateral dimension as in Fig. 2 so that the semi-infinite assumption of the analytical model is justified. These additional sample-by-sample measurements yield nearly the same *k* results, within 5%, as the paralleled measurement (see the representative result for BGO in supplementary material Sec. S8 4).

FEM is not always necessary for a paralleled SI-TI measurement even when sample edges lie within the FOV; sometimes, the simpler and faster analytical model, which assumes a semi-infinite sample, may also give acceptable results. To explore this, we applied the analytical model to the data from the parallelized measurement of each sample separately. As shown in Fig. S37 and discussed in supplementary material Sec. 8.5, the model fits yield relatively small RMS errors for the lower-*k* samples (namely, PMN-PT and PbF_{2}), yielding *k* values that are lower but close to those from FEM: *k*_{PMN-PT} = 1.32 ± 0.12 and *k*_{PbF2} = 1.27 ± 0.14 W m^{−1} K^{−1}. On the other hand, the fitting error becomes much larger for BGO and the result is unreliable. We believe this is because BGO has larger *k*, resulting in larger $dpi$ so that the edge effects which are neglected in the analytical model are no longer negligible. It is expected that if the heater is placed further away from the edge, and/or a higher heating frequency is used to suppress the in-plane thermal penetration length, the analytical model should still be applicable to the BGO measurement of Fig. 5(d). Furthermore, although beyond the scope of this work, we note that in principle the analytical transfer matrix model itself could be further modified to include finite sample size effects.^{87,88}

## IV. DISCUSSION

### A. Demonstrated advantages of SI-TI

The results above demonstrated SI-TI's capability for paralleled measurements of a wide variety of samples, including rough and anisotropic samples. Looking forward, we believe the SI-TI technique also holds promise for H-T thermal metrology, which will be important to complement H-T theory and synthesis, such as samples from combinatorial manufacturing as shown in Fig. 4(c), in the materials by design development cycle.^{1,2,89} From this perspective, one strength of SI-TI is the relatively simple sample preparation, which requires only general flattening to surface roughness of a few micrometers [recall the rough samples of Figs. 4(e) and 4(f)], rather than the more stringent requirement for optically flat surfaces in methods like TDTR and FDTR. Furthermore, SI-TI samples require a single metal evaporation step without patterning needed, which is significantly easier than the microfabrication requirements of the electrothermal 3ω method.^{20,30} In principle, the metal transducer layer might even be omitted entirely for certain samples as has been recently demonstrated in another TR method.^{90,91} SI-TI also works for small, mm-scale samples, which is an advantage compared to some traditional methods which require centimeter scale samples.^{21,22,92}

Another strength of SI-TI is the adaptivity of the heating patterns to probe thermally anisotropic samples, as demonstrated in Fig. 3 for mica. Thus, the method can resolve the *k* tensor purely by software control of the heating patterns, without any need to physically change the optics to modulate the beam shape as in TDTR^{65} and Raman thermometry methods, or rotating the sample carefully in TDTR/FDTR-based beam-offset methods.^{30,93} In principle, SI-TI could be extended to measure far more than three samples in parallel, by using lower-magnification optics to increase the FOV and installing a correspondingly higher power optical heating source. Extending the ideas of Fig. 5 would also enable spatial mapping of $kx,y$ in inhomogeneous materials like composites or diffusion couples.

In the current setup, if maximum throughput is the priority, supplementary material Sec. S7 3 demonstrates *k* extraction using data from just one heating frequency using Bayesian parametric inference, specifically the highest frequency because it has the greatest sensitivity to *k*. Table S4 compares the tradeoffs between the number of cycles averaged and the number of unique heating frequencies, for a glass sample. This comparison shows that analyzing single frequency data from 250 cycles of averaging gives results for *k* very similar to that from measurements at five different frequencies each with 250 cycles (1250 cycles total), namely, $kglass\xb12\sigma =1.43\u22120.310.34$ and $kglass\xb12\sigma =1.42\u22120.250.29$ W m^{−1} K^{−1}, respectively. In other words, this major reduction in the total measurement time by a factor of 5 only changed the median value of *k* by 1%, and only slightly increased the total breadth of the ±2σ uncertainty band, by 20% (0.11 W m^{−1} K^{−1}). In addition, with rich information from thermal imaging, we show in supplementary material Sec. S7 3 2 using the polished YSZ as an example that analyzing data of 10 ROIs can further improve the accuracy of single frequency measurements, yielding $kYSZ\xb12\sigma =2.07\u22120.230.26$ W m^{−1} K^{−1}, as compared with the result using 5 ROIs, $kYSZ\xb12\sigma =1.98\u22120.270.30$ W m^{−1} K^{−1}.

### B. Limitations of the current system and potential extensions

It is important to note several limitations of the SI-TI method as demonstrated here. The greatest limitations arise fundamentally from the low heating power reaching the sample, which is due to the relatively modest power of the LED source and several challenging transmission losses in the microscope system. This weak power then necessitates a high magnification objective lens for higher power density, low measurement frequencies, and low-*k* materials, all of which tend to increase the amplitude of the temperature response; we have found that ∼1 K peak-to-peak is generally adequate. Note the low heating frequencies also results in a relatively large in-plane thermal penetration length $dpi$. This then leads a secondary limitation on the property mapping spatial resolution of $kx,y$, which is far coarser than the intrinsic resolution of the optical imaging system (supplementary material Sec. S2 1).

This work demonstrated paralleled measurement of nine thermally independent ROIs (which can be regarded in a sense as nine different samples) in Fig. 5(c) and supplementary material Sec. S8 4. In principle, the total number of thermally independent measurements zones could be extended much higher. For a best-case scenario of the current setup, a highly favorable material would have ultralow *k*, e.g., *k _{s}* = 0.02 W m

^{−1}K

^{−1}and

*C*

_{s}= 0.2 MJ m

^{−3}K

^{−1}. With a SI-TI measurement configuration similar to Fig. 5 and

*f*

_{0}= 10 Hz, we can achieve a short $dpi$∼ 56

*μ*m with a reasonable ∼1 K amplitude of the temperature oscillation. We assume square heaters of size $dpi$ × $dpi$, separated by gaps of width $dpi$, organized in a square array. For the 20× lens, this yields 15 thermally independent heating patterns within its 357 × 572

*μ*m FOV; and similarly for the 5× lens, up to 260 independent measurements could be obtained simultaneously within a single FOV, giving an approximate upper limit for paralleled measurements in the current SI-TI system.

Total measurement time is another important factor for measurement. Considering their diverse approaches to thermometry, operational protocol, and sample preparation, it is difficult to quantitatively compare the measurement throughput of SI-TI to that of alternative methods,^{17–22,27,28,39,95,96} though we provide a semi-quantitative comparison in supplementary material Sec. S1 2. Additional analysis of the experimental noise in SI-TI (supplementary material Sec. S2 5) reveals that the SNR and thus throughput is currently limited mainly by the relatively low optical pump and probe power density at the sample plane. For the measurements reported above on samples with *k*∼ 1 W m^{−1} K^{−1}, we used relatively long averaging times to minimize this experimental noise. As shown in supplementary material S7 2, averaging 250–300 cycles (each of duration *t*_{p} = 3.5 s, thus totaling ∼14–17 min) can increase the third harmonic SNR to >10 (see supplementary material Sec. S3 4) for the glass sample. For materials with *k* $\u226a$ 1 W m^{−1} K^{−1}, such as the aerogel of Fig. S30, ∼20 cycles (∼1 min) are enough to achieve a SNR > 30.

Looking ahead, in a next generation of SI-TI, the pump power could in principle be increased by one or two orders of magnitude by combining a higher-power light source with free space optics. This would enable much higher heating frequencies while still maintaining a reasonable temperature response of ∼1 K amplitude and thus improve the measurement throughput and *k* mapping resolution. A high-speed camera or heterodyne detection^{94} might also be needed for the optical thermometry at higher frequency. Assuming a heating frequency of 1000 Hz which is typical for a DMD refresh rate, the $dpi$ and $dpi$ will be ∼15 *μ*m for glass and ∼160 *μ*m for silicon, two materials which approximately bound many other common materials (e.g., ceramics and composites). Assuming a square matrix geometry for the pump spots with the heaters size $dpi$ × $dpi$ and a gap width of $dpi$ [similar to Fig. 5(c)], with the requirement of at least 1 K temperature oscillation ($\Delta T$) for measurement, the optical power density at the sample will need to be approximately $\Delta T2\pi fksCs\u223c$ 0.14 W mm^{−2} for glass and 1.2 W mm^{−2} for Si. With a low magnification lens such that the FOV is ∼(5 mm),^{2} the corresponding total optical power required at the sample ranges from ∼1.4 W for glass to 12 W for Si, considering ¼ of the area is illuminated and an Au absorption *α *= 0.62. A total power of 12 W is viable for some existing commercial DMDs. With thermal spatial resolution of $2dpi$, the corresponding total number of thermally independent pixels that can be interrogated in parallel ranges from ∼240 for Si to over 27 700 for glass. Higher-frequency DMDs coupled with high power pump power could further increase these values. With reduced $dpc$, it is also reasonable to envision the application of SI-TI in studying film materials with thickness down to <5 *μ*m. Another extension of SI-TI would be to generate different heating frequencies (waveforms) for different pump pixels which can enhance the throughput and spatial resolution (see supplementary material Sec. S1 2).

## V. CONCLUSION

We report a new all-optical thermal property characterization technique based on structured illumination with thermal imaging, which combines the advantages of a digital micromirror array for spatial control of the heating pattern with thermoreflectance for full-field thermometry. We demonstrate the accuracy and unique power of SI-TI through measurements of diverse forms of materials, from dense isotropic and anisotropic crystals to porous 3D printed films. SI-TI overcomes challenges with many traditional thermal metrologies if applied for next-generation flexible and high-throughput thermal measurements, by enabling sample-adaptive reconfigurable heating and unprecedented parallelized measurement capability. Together with its relatively high tolerance of sample roughness and simple sample preparation, this method can open a new avenue in adaptivity and throughput for thermal characterization of diverse materials.

## SUPPLEMENTARY MATERIAL

See the supplementary material for extensive additional information about the experimental method and data analysis.

## ACKNOWLEDGMENTS

This work was supported principally by the Building Technologies Program in the Office of Energy Efficiency and Renewable Energy of the U.S. Department of Energy under Contract No. DEAC02–05CH11231. Portions of this work were performed at the San Diego Nanotechnology Infrastructure (SDNI) of UCSD, a member of the National Nanotechnology Coordinated Infrastructure (NNCI), which is supported by the National Science Foundation (NSF) Grant No. ECCS-1542148. Partial funding by NSF Grant No. 1545852 to UCSD is also acknowledged. The thermoelectric film printing performed at the University of Notre Dame was funded by the U.S. Department of Energy under Award No. DE-EE0009103.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

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

### APPENDIX: METHODS

##### 1. Optical components in the SI-TI system

The SI-TI system here is composed of a MTIR120 Lock-In Thermal Imaging System thermoreflectance imaging system from Microsanj which is built upon a top imaging (i.e., epi-illumination) industrial Olympus BXFM Modular Microscope with infinity corrected optics. The whole system is placed on an optical table to damp vibration. An Olympus LMPlanFL N 20x/0.40 objective lens is used to heat the sample and collect the thermoreflectance image. A Mightex high-power LED collimator source with 470 nm wavelength, LCS-0470-50-22, is used for the pump light. A liquid light guide with NA = 0.59 is used to carry the pump light from the source to the DMD. The microscope field of view is calibrated using an NBS 1963A resolution target. No noticeable distortion of the image is observed near the edge of the image.

A Mightex Polygon400 DMD illuminator with 684 × 608 pixels with vertical and horizontal pitch size of 7.64 *μ*m is used to modulate the LED heating light with software controlled patterns. The on or off state of each DMD pixel is determined by controlling an array of opto-mechanical aluminum micromirrors which can individually switch between two discrete angular positions of ±12° with respect to the surface normal to change the direction of the light deflection.

Three singlet spherical focal lenses with focal lengths of 10, 10, and 15 cm and a diameter of 10 cm are used after the DMD output to modulate the shape of the pump light rays before they enter the microscope system. This improves the power density of the pump light heating while maintaining sharp edges of pump patterns. The pump power of each heating pattern (or a component of it) is measured by a Thorlabs S120C photodiode power sensor.

A Thorlabs M530L3 LED at 530 nm driven in the CW mode is used as probe light for TR signal detection. The thermoreflectance signal is collected by a 1920 × 1200 pixel Si CMOS camera with 5.86 *μ*m pitch size which is mounted onto the microscope using a Olympus U-TV1X-2 C-mount adaptor with a telecentric tube lens. A Rigol DG4062 function generator is used to provide the phase-locked external trigger pulse for the CMOS camera and the modulation of the pump light (or electrical excitation in supplementary material Sec. S8 1). A Semrock 519 nm cut-on blocking edge BrightLine^{®} long-pass filter and a 532 nm cut-on EdgeBasic^{®} long-pass filter are used to further filter the reflected pump light at long wavelength.

##### 2. Temperature calculation in thermoreflectance thermal imaging

The temperature map $\Delta T(x,y,t|f0)$ is derived from the raw image $I(x,y,t|f0)$ at time *t* from the 2M pixel CMOS camera by the TR effect using the following equations ($f0$ is omitted for brevity below). $I(x,y,t)$ corresponds to a frame among the series of 210 frames in the 3.5 s window which is averaged with repeated measurement in a time window of same length (see supplementary material Sec. S2 5). First, the processed image intensity $Ip(x,y,t)$ of the TR light is obtained through the normalization of $I(x,y,t)$ as

where $I\u2009x,y\xaf$ is the time-averaged image of all 210 frames over the time window, $Rx,y,t$ is the sample surface reflectance, and $PI(x,y)$ is the probe illumination power which is constant in time and cancels out in the calculation. Thus, $Ip(x,y,t)$ is eventually proportional to $\Delta Rx,y,t$, the difference between $Rx,y,t$ and the average reflectance of all frames in this time window, $Rx,y\xaf$, all normalized by $Rx,y\xaf$. Next, through the TR effect, $Ip(x,y,t)$ is related to the sample surface temperature response $\Delta Tx,y,t$,

where $R0$ is the room temperature reflectance of the gold transducer layer (that is, we take $R(x,y\xaf)\u2248R0\u2248const.$ over the entire measurement range, which is nearly exact since the steady state temperature rise in our measurement is <5 K and the gold TR coefficient is approximately constant),^{44,58} and $Cth$ is the normalized TR coefficient (see calibration in supplementary material Sec. S2 2).

##### 3. Data analysis and thermal model

We briefly introduce the analytical model here, with more details given in supplementary material Sec. S4. The details of the numerical finite element model are given in supplementary material Sec. S5.

The analytical thermal model used in the analysis is based on the frequency domain solution of the boundary condition problem of heat diffusion in multilayered structure with an arbitrary number of surface heaters (1, 2, 3, …., *N*), any *j*th among which has a finite power only within a rectangular area with half width $wj$ and half height $hj$.^{25,65,66} The analytical form of the pump power distribution $pjx,y,t$ is derived from a quadratic fitting of the experimentally measured results as detailed in supplementary material Sec. S2 1. In brief, the spatially averaged temperature response for a rectangular-shaped ROI_{i} (*i *=* *1, 2, 3, …, M) with half width *w*_{i} and half height *h*_{i} and a center shift from the *j*th heater center of $\Delta xij,\Delta yij,$ denoted $\Delta Ti\Delta xij,\Delta yij,f$, is given by the inverse spatial Fourier transform of the product of the spatially non-uniform *j*th pump heating $Pju,v,f=F(pjx,y,t)$ and the sample Green's function $G\u0302u,v,f$, where *u* and *v* are Fourier vectors corresponding to *x* and *y*, as follows:

Here, $G\u0302u,v,f$ is derived from a transfer matrix calculation of the multilayer structure which is a function of the *k, C*, and thickness *L* of each layer including the interface(s). The frequency domain expression of the *periodic* pump heater $Pju,v,f$ is a modulated Dirac comb in the temporal frequency space with non-zero values only at discrete frequencies of $fn=nfj0$, *n *=* *1, 2, 3, …,.

Define

and

Then, the temperature response of the *i*th ROI at a certain FFT pickup frequency

Thus, the temperature of a certain ROI is the summation of the individual contributions from all heaters. The summation in $\Delta Ti$ for $\Delta xij,\Delta yij$ over *j* is omitted for conciseness although it is a function of distances from all the heaters. A Clenshaw–Curtis quadrature based on Boyd's method^{97} is used to numerically calculate the integral over an infinite domain. Since we only use data at discrete frequencies after the FFT, $\Delta Ti,fn$ is also represented as $\Delta Ti(f)$ for simplicity in the main text.

The thermal properties are found by fitting the frequency domain signals $\Delta Ti(f)$, using experimentally determined heater and ROI geometries (*w*_{j}, *h*_{j}), (*w*_{i}, *h*_{i}), and the center shift $(\Delta xij,\Delta yij)$. Although not used for the property fitting, when the model curve is depicted in the time domain [as in Fig. 2(d)], to avoid the Gibbs phenomenon^{98} we sum the Fourier terms using Fejér's theorem.^{99}

##### 4. Transducer layer deposition for SI-TI measurement

Gold thin film deposition is carried out by electron-beam physical vapor deposition method using the CHA Solution E-beam Evaporator in the UC Berkeley Marvell Nanofabrication Laboratory. For each deposition, a 10 nm Ti layer is deposited first to promote the adhesion of gold on the sample. The transducer layer thickness is measured using a control silicon wafer placed next to the samples in the deposition process by stylus profilometry and verified by a laser confocal microscopy. The total *k* of the layer is obtained by a combination of four-probe electrical resistivity measurement and Wiedemann–Franz law and literature value of lattice thermal conductivity (see supplementary material Sec. S3).

##### 5. Synthesis and extrusion printing of thermoelectric Bi_{2}Te_{2.7}Se_{0.3} thin film

First, the Bi_{2}Te_{2.7}Se_{0.3} nanoplates are prepared using a modified wet-chemistry approach.^{100–102} In a typical synthesis, 18 mmol bismuth (III) nitrate pentahydrate (Bi(NO_{3})_{3}·5H_{2}O), 27 mmol sodium tellurite (Na_{2}TeO_{3}), 5 g sodium hydroxide (NaOH), and 2.4 g polyvinylpyrrolidone (PVP) are added to 150 ml ethylene glycol. The mixture is heated under reflux overnight before the particles are collected by centrifugation at 5000 rpm. The final products are further cleaned by washing with ethanol three times, followed by drying under vacuum for 48 h.

To form the ink for extrusion printing, 1.5 g of Bi_{2}Te_{2.7}Se_{0.3} nanoplates is crushed using a mortar and then dispersed in ethylene glycol, forming 40 wt. % Bi_{2}Te_{2.7}Se_{0.3} inks. A small amount of 2-mercaptoethanol (0.2 wt. %) is added to improve the performance of Bi_{2}Te_{2.7}Se_{0.3} inks. The resulting dispersions are sonicated for 15 min, followed by constant homogenization using a vortex machine (VWR^{®} Analog Vortex Mixer, USA). Once the ink is homogenized, a custom-built pneumatic 3D extrusion printer is used to fabricate Bi_{2}Te_{2.7}Se_{0.3} films. To ensure the film uniformity, a systemic optimization of the extrusion pressure is performed using a 22G nozzle, where the film size is set to be 15 × 2 mm^{2} under a 3.5 mm/s speed. The printed TE film is then dried in a tube furnace at 250 °C for 1 h, followed by the cold pressing twice at 20 MPa for 5 min. Finally, the pressed TE film is then sintered at 440 °C for 1 h in an inert environment.

##### 6. Fabrication of yttria-stabilized zirconia (YSZ) sample

Commercial 8 mol. % yttria-stabilized zirconia (YSZ) powder procured from Tosoh Corporation (TZ-8Y, Tokyo, Japan) was used to fabricate bulk ceramic samples via the Current-Activated Pressure-Assisted Densification (CAPAD) processing technique.^{76} Approximately 1.50 g of as-received powder was packed into a cylindrical graphite die (with inner and outer diameters of 19 and 45 mm) which was then placed inside the CAPAD chamber, secured between two graphite and copper electrodes, and evacuated to a vacuum of 10^{−3 }Torr prior to densification. Densification was performed using a two-step pressure ramp, consisting of an initial uniaxial external applied stress of 106 MPa, followed by a second ramp to a terminal applied stress of 141 MPa. The green body was subject to an average heating rate of 260 °C min^{−1} to a target temperature of 1200 °C, and a dwell-time of 10 min.

## References

*Conduction of Heat in Solids*, 2nd ed. (Oxford University Press, London,