1. Introduction

With the increased miniaturization and performance of electronics^[1-3], thermal management has become a critical issue. The local heat flux near hot spots can exceed the average heat flux by an order of magnitude^[4] and induce significant thermal stress, resulting in accelerated aging and performance degradation^[5]. To mitigate these challenges, thermal interface materials (TIMs) are employed to enhance the interfacial heat transport between the electronics and the heat sink^[6-8]. In recent years, the thermal conductivity of TIM has been significantly enhanced through microstructure engineering strategies, such as introducing interconnected filler networks^[9-12] and improving filler-matrix and filler-filler interfacial thermal transport via surface functionalization, covalent bonding^[13-15], and noncovalent modifications^[16-18]. Nevertheless, the lack of thermal property visualization tools has severely limited the understanding of the structure-property correspondence of TIMs with micro/nano-architectures, which is crucial for developing next-generation thermal management materials. Thermal property mapping can also serve as datasets for training machine learning models, for AI-assisted optimization and design of TIMs with superb heat transfer performance.

However, commonly used thermal metrology methods for characterizing TIMs, either the steady-state measurement based on ASTM D5470^[19-21], or transient methods such as laser flash analysis^[22] and photoacoustic measurements^[23,24], are primarily developed to measure the overall macroscale thermal conductivity or interfacial resistances without spatial resolution. While infrared microscopy methods can image temperature maps of the cross-sections of TIMs under heating load, extracting local thermal conductivity and interfacial resistances remains a very challenging inverse problem^[25,26]. Spatial thermal property mapping requires local measurements on micrometer scales. Pump-probe thermoreflectance techniques, including time-domain thermoreflectance (TDTR)^[27,28] and frequency-domain thermoreflectance (FDTR)^[29], rely on depositing a metal transducer layer on the sample surface, and using the pump and probe beams for laser heating and sensing of the thermoreflectance signal. TDTR and FDTR have been widely used to characterize thermal properties of bulk semiconductors^[30,31], thin films^[32,33], and van der Waals materials^[34], as well as TIMs^[35-38]. This noncontact, local heating and sensing scheme enables thermal property imaging through point-by-point scanning across the sample surface^[39-41]. However, when performing thermal property imaging, the sample surface must be smooth and flat with a roughness below 10 nm and free of surface wrinkles. This is extremely hard to achieve when depositing a transducer layer on TIMs with mechanical flexibility and fluidity.

In this work, we employ FDTR microscopy to visualize the microscale thermal-property heterogeneity of particle-filled TIMs. Unlike conventional TIM measurements that yield spatially averaged thermal resistance or thermal conductivity, our approach resolves the local thermal conductivity and interfacial thermal conductance within a sandwiched TIM layer. We observed that the local thermal properties in an Al₂O₃/ZnO-loaded thermal grease are highly non-uniform, exhibiting clustered high-conductivity domains and localized low-conductivity defects. Correlative micro-CT imaging further demonstrates that these thermal features are closely associated with the filler architecture, specifically, filler aggregation, filler-deficient regions, and pores. We also perform pressure-dependent FDTR imaging during loading and unloading to examine how averaged thermal properties, local conductive features, and interfacial contact evolve under compression. We further performed statistical analysis of the property image and find that the thermal conductivity values deviate from the normal distribution, and instead exhibit an asymmetric long-tailed probabilistic distribution. We develop a statistical model clarifying that the asymmetry in thermal-property distribution originates from the lognormal distribution of particle size and spatial randomness of filler particles. This work extends the capability of FDTR microscopy for profiling the thermal property of thermal management materials at microscales.

2. Methods

2.1 Measurement setup

Thermal property mapping of TIMs is performed using FDTR, as illustrated in Figure 1a. In FDTR measurements, a modulated pump laser generates a diffusive thermal wave, while a probe laser captures the harmonic response from the reflected beam. The pump beam with a wavelength of 488 nm for heating the sample is modulated with a sinusoidal signal using a Zurich Instruments HF2LI lock-in amplifier. The probe beam (with a 532 nm wavelength) is coaxially aligned with the pump beam and focused on the same spot on the sample surface for detecting the thermoreflectance signal. To enhance the signal-to-noise ratio in the low-frequency regime, a balanced photodetector comprised of two well-matched photodiodes, PD1 and PD2, is used. A polarizing beamsplitter divides the probe beam into a signal path, which is focused onto and reflected from the sample toward PD1, and a reference path to PD2. By rotating a half-wave plate to balance the intensities of PD1 and PD2 and differencing their photocurrents within the detector, we effectively suppress the common-mode noise inherent in the probe laser.

Display Full Size

Download

Figure 1. Thermal property mapping system for TIMs. (a) Optical configuration of the FDTR scanning system; (b) Schematic of thermal property scanning and sample configuration for TIMs characterization; (c) Phase images captured at different modulation frequencies. TIMs: thermal interface materials; FDTR: frequency-domain thermoreflectance; PBS: polarizing beamsplitter; ISO: isolator.

To ensure an optically flat and smooth transducer surface, we prepare the TIM sample in a sandwich fixture, as illustrated in Figure 1b. Scanning electron microscopy (SEM) and the corresponding energy dispersive spectroscopy images reveal that the TIM consists of a polydimethylsiloxane (PDMS) matrix incorporated with micrometer-sized Al₂O₃ particles and ZnO nanoparticles (Figure S1). A 110 nm gold transducer film was deposited by electron-beam evaporation on a double-side-polished fused silica slide purchased from MTI Corporation (Ra < 5 Å). The gold thickness is well above 50 nm to ensure optical opacity of the transducer layer. The low roughness of the polished fused silica slide provides a smooth reflective surface for FDTR measurements. The TIM was first applied onto the silicon substrate and leveled using a doctor blade to ensure an initially uniform coating. The transducer-coated silica slide was then placed onto the TIM surface and compressed using the fixture with four spring-loaded screws. The contact pressure was estimated to be approximately 36 psi based on the force measured by a force sensor and the contact area. The bondline thickness of the TIM layer was determined to be approximately 12.8 μm from micro-CT cross-sectional images, which also confirmed that the thickness is uniform across the imaging area (Figure S2). The sample was mounted on a piezoelectric stage with a 200 μm XYZ travel range and a 1 nm stepping precision. Both the pump and probe beams are focused through the silica substrate onto the gold transducer surface. Thermal property images were generated by performing point-by-point scans across the mapping area. At each scanning location, the phase of the thermal response was recorded across eight modulation frequencies using a lock-in amplifier. The phase images captured at different modulation frequencies are shown in Figure 1c.

2.2 Local thermal property measurements

Our FDTR microscopy set-up builds the thermal property image in a pixel-by-pixel manner, and thermal properties are extracted from local frequency-sweep measurements. Different from the conventional geometry with the transducer on the top layer, the pump beam in our setup passes through the silica layer and the generated heat will propagate bidirectionally, as shown in Figure 2a. We therefore use the bidirectional thermal model^[42] to extract thermal properties from FDTR signals. For each layer, the temperature rises (θ) and heat fluxes (Q) at the top and bottom surfaces are related by a thermal transfer matrix M:

Display Full Size

Download

Figure 2. Single point FDTR measurement on the TIM. (a) Thermal model for bidirectional heat conduction into the TIM and the silica slide; (b) Calculated sensitivity of thermal phase ϕ to thermal properties of the TIM; (c) Frequency-sweep FDTR signal obtained on the TIM sample with the best fit value k_TIM = 16.01 W m^-1 K^-1 and G_Au/TIM = 55.11 MW m^-2 K^-1; (d) Confidence interval for the FDTR measurement on the TIM. FDTR: frequency-domain thermoreflectance; TIM: thermal interface material.

(1)$\left[\begin{array}{c}{\theta }_b\\ Q_b\end{array}\right]=M\left[\begin{array}{c}{\theta }_t\\ Q_t\end{array}\right]=\left[\begin{array}{cc}\cosh \left(qd\right)& -\frac{1}{kq}\sinh \left(qd\right)\\ -kq\sinh \left(qd\right)& \cosh \left(qd\right)\end{array}\right]\left[\begin{array}{c}{\theta }_t\\ Q_t\end{array}\right]$

where the subscripts b and t represent the bottom and top surfaces, respectively; d is the thickness of the sample layer, and $q^2=i\omega C/k+{\xi }^2$ denotes the squared thermal wavenumber, with C and k the volumetric specific heat and thermal conductivity respectively, and ξ is the Hankel transform variable in the radial direction. Note that the thermal conductivity is treated as an isotropic scalar property when computing the thermal wavenumber, because the TIM materials we considered in this work contain randomly dispersed Al₂O₃ microparticles and ZnO nanoparticles in a polymer matrix, without intentionally aligned fillers or directional thermal pathways. For anisotropic TIMs embedded with aligned flakes of two-dimensional materials, or carbon nanotube arrays, thermal anisotropy can be measured using measurements with varied laser spot sizes, offset pump-probe beams^[43], or elliptical beams^[44], and the anisotropic thermal models developed in previous works can be used for data reduction. For the heat conduction downwards across the TIM layer into the substrate, one can multiply the transfer matrices together:

(2)$\left[\begin{array}{c}{\theta }_b^{+}\\ Q_b^{+}\end{array}\right]=M_n{T}_{n,n-1}{M}_{n-1}\cdots T_{21}{M}_1=\left[\begin{array}{cc}{A}_{+}& B_{+}\\ C_{+}& D_{+}\end{array}\right]\left[\begin{array}{c}{\theta }_t\\ Q_t^{+}\end{array}\right]$

where the superscript + denotes the heat flux propagating in the +z direction (downward), M_n is the transfer matrix of the n-th layer defined in Eq. (1), and T_n,n-1 is the transfer matrix of the interface between the (n-1)-th and the n-th layer:

(3)$T_{n,n-1}=\left[\begin{array}{cc}1& -1/G_{n,n-1}\\ 0& 1\end{array}\right]$

where G_n,n-1 is the interfacial thermal conductance. Similarly, the upward heat conduction can also be written as follows:

(4)$\left[\begin{array}{c}{\theta }_b^{-}\\ Q_b^{-}\end{array}\right]=\left[\begin{array}{cc}{A}_{-}& B_{-}\\ C_{-}& D_{-}\end{array}\right]\left[\begin{array}{c}{\theta }_t\\ Q_t^{-}\end{array}\right]$

An adiabatic boundary condition is assumed at the bottom and top surfaces of the fixture^[42], therefore $Q_t^{\pm }=-D_{\pm }/C_{\pm }{\theta }_t$. The total amount of heat deposited into the transducer is $Q_t=Q_t^{+}+Q_t^{-}$, we can therefore obtain the Green’s function as:

(5)$G=\frac{{\theta }_t}{Q_t}=\frac{-D_{+}{D}_{-}}{D_{+}{C}_{-}+D_{-}{C}_{+}}$

which represents the thermal response to a unit point source. In FDTR measurements, both the pump and the probe beams have a Gaussian profile of intensity, and the complex thermal response function is calculated as:

(6)$H\left(\omega \right)=\frac{\beta }{2\pi }{\int }_0^{\infty }G\exp \left[\frac{-(w_{\text{pump}}^2+w_{\text{probe}}^2){\xi }^2}{8}\right]d\xi$

where w_pump and w_probe are the 1/e² radii of the pump and probe laser spots respectively, and β is a constant proportional to the laser power, thermoreflectance coefficient, and the gain of the photodetector. The phase angle of H(ω) is typically used to extract thermal properties.

(7)$\varphi \left(\omega \right)=\mathrm{atan}\left(\mathrm{Im}\right[H\left(\omega \right)\left]/\mathrm{Re}\right[H\left(\omega \right)\left]\right)$

Thermal conductivity and the interfacial thermal conductance are extracted by fitting the predicted phase to the experimental phase using nonlinear least squares regression.

Figure 2b shows the sensitivity analysis of the phase signal to different thermal properties. Sensitivity is defined as the logarithmic derivative S_x = ∂lnϕ/∂lnx^[45], where x denotes the fitting or input parameter in the thermal model. Parameters used to calculate the sensitivities are calibrated and listed in Table S1. The frequency-dependent sensitivities arise from the variation in sensing depth with the modulation frequency of the pump laser. While the actual sensing depth is simultaneously determined by the laser spot size and the heating frequency^[46,47], our experiments are based on a fixed spot size, therefore we can use the simplified penetration length $d_p=\sqrt{k/\left(\pi Cf\right)}$ to qualitatively understand the frequency-modulated depth over which the sample is thermally perturbed by the laser beams. Under very low frequency ~1 kHz, the penetration length can reach 46 μm, resulting in low sensitivity to interface conductance, and the phase signal is mainly determined by thermal diffusion in the sample layer. To achieve simultaneous extraction of thermal conductivity and the interface conductance, we select 0.1-1 MHz to measure the thermal phase signals, where sensitivities to thermal conductivity and interface conductance are acceptable for data reduction. The sensitivity analysis also shows that the phase signal has relatively low sensitivity to the gold film thickness and thermal conductivity compared with the TIM thermal conductivity and interfacial thermal conductance. Therefore, uncertainties in the gold transducer parameters are not expected to significantly affect the extracted TIM thermal properties.

A representative FDTR signal of the TIM sample is shown in Figure 2c. Multi-parameter fitting via least-squares minimization yields k_TIM = 16.01 W m^-1 K^-1 and G_Au/TIM = 55.11 MW m^-‍2 K^-1. An uncertainty analysis of the extracted thermal properties is performed by implementing the multivariate error propagation scheme based on the Jacobian matrix^[48]. The confidence interval has an elliptical shape as shown in Figure 2d, the relative uncertainty of the thermal conductivity of the TIM k_TIM is 7.73%, and that of the interfacial thermal conductance G_Au/TIM is 17.95%. These results are consistent with reported values^[49] and manufacturer specifications.

3. Results and Discussion

3.1 Thermal property mapping

Figure 3a,b present the mapping of thermal conductivity and interfacial conductance, respectively. The results reveal significant heterogeneity at the microscale with a randomly scattered distribution of high-thermal-conductivity regions. The area-averaged thermal conductivity calculated from our mapping is 12.97 W m^-1 K^-1, which is comparable to the macroscopic value of 13.64 W m^-1 K^-1 measured via ASTM D5470^[49]. Local peak values reach nearly 20 W m^-1 K^-1, while in voids or defect areas within the sample, the thermal conductivity drops to approximately 0.1 W m^-1 K^-1. The spatial distribution of interfacial thermal conductance G_TIM exhibits similar microscale heterogeneity. Regions with high thermal conductivity generally correspond to high interfacial conductance. This positive correlation occurs because local agglomerations of conductive fillers not only enhance local heat transfer but also create more effective contact points at the interface^[50].

Display Full Size

Download

Figure 3. Thermal conductivity mapping of the TIM and its correlation with microstructure. (a) Spatial distribution of thermal conductivity and (b) interfacial thermal conductance; (c) The corresponding micro-CT slice and (d) threshold-segmented 3D reconstruction of the contact area, with the dimensionless intensity values obtained by linearly normalizing the CT grayscale signals from 0 to 1 for visualization. TIM: thermal interface material; CT: computed tomography.

To achieve structural visualization of the same area used for FDTR mapping under the same loading, we used micro-CT imaging that is able to penetrate across the transducer layers without disassembling the sample fixture. The co-registration of the micro-CT image and the FDTR image is based on aligning the edges of the square-shaped transducer films. Alignment details are presented in Section S5, with an alignment uncertainty of approximately 0.6 μm. As shown in Figure 3c, the thermally conductive filler particles correspond to the areas with high intensity values, because the higher density of particles attenuates X-rays more strongly than the polymer matrix^[51]. Figure 3d shows the 3D reconstruction image obtained by applying a threshold segmentation to Figure 3c. The blue regions represent aggregated fillers, and the white regions indicate areas with lower densities containing voids or low concentrations of oxide fillers.

To further demonstrate the structural-property correspondence, we selected three representative regions for further comparison: H1 with lots of scattered areas of high local thermal conductivity, L1 covering the assembly void defect, and L2 with relatively low local thermal conductivity. The CT images of region H1 indicate a filler-rich, spatially connected structure, suggesting percolated local heat-conduction pathways. In contrast, region L1 contains low-density voids with very low local thermal conductivity. These voids are introduced due to trapped air during the fixture assembly, as indicated by the low-density area detected by the CT imaging, see Figure 3c,d. The agreement between the thermal conductivity map and CT image demonstrates the capability of detecting defects in TIMs using FDTR microscopy. The L2 region exhibits lower thermal conductivity than the average, manifested as lower grayscale intensity in the CT slice and a filler-deficient structure in the reconstructed 3D image. These areas are characterized by a sub-average filler concentration where the isolated particles fail to establish percolating networks. Interestingly, the spatial mapping of G_TIM also reveals localized high-conductance domains that do not consistently overlap with high-k regions and exhibit smaller cluster sizes. We attribute this discrepancy to the intimate mechanical contact between individual Al₂O₃ particles and the gold transducer under compression, which locally enhances G_TIM even in regions with moderate bulk thermal conductivity^[50].

3.2 Thermal property evolution under varying contact pressures

To investigate the effect of contact pressure on the thermal properties of TIMs, FDTR imaging was performed at different pressures during both loading and unloading processes. The contact pressure was increased from 10 psi to 40 psi and then released to 0 psi. The measurement commenced at 10 psi because a minimum compressive load is required to flatten the thermal grease and establish reliable contact between the TIM and the transducer. The FDTR measurements at each pressure were conducted using a 10× objective with pump and probe beam radii of approximately 2.5 μm along the x-direction and 2.7 μm along the y-direction, over a modulation frequency range of 0.1-1 MHz.

Figure 4a,b show the area-averaged thermal conductivity and interfacial thermal conductance as a function of contact pressure, respectively. During loading and unloading, the averaged thermal conductivity remains nearly constant from 10 to 40 psi, and decreases when the pressure is further unloaded to zero. On the other hand, the interface conductance shows weak pressure dependence, possibly because the fluidity of the thermal grease ensured sufficient contact with the transducer layer.

Display Full Size

Download

Figure 4. Thermal property evolution of the TIM under varying contact pressures. (a) Thermal conductivity and (b) interfacial thermal conductance as a function of contact pressure during loading and unloading processes; (c-h) Thermal conductivity maps at different contact pressures during loading and unloading processes. TIM: thermal interface material.

Figure 4c,d,e,f,g,h further visualize the evolution of the local thermal-conductivity maps during loading and unloading. When the loading is increased from 10 to 20 psi, the scattered pattern of high-thermal conductivity areas is re-arranged, as indicated in the squared areas in Figure 4c,d. Interestingly, further increasing the loading does not rearrange the high-thermal conductivity patterns, instead, these features appear to shift slightly, which may be associated with pressure-induced creep or lateral displacement of the thermal grease, see squared areas in Figure 4d,e. Similar flowing behavior of high-thermal conductivity region is also observed in the unloading cycle (Figure 4f,g), and the high-thermal conductivity areas are re-arranged again when the pressure is further reduced from 20 psi to 10 psi (Figure 4g,h). Such behavior is consistent with the weak pressure-dependent thermal conductivity: the translation of thermally conductive particles does not alter the heat conductive pathways from 20 psi to 40 psi. These imaging results demonstrate that FDTR can be used as an in situ thermal characterization tool for TIMs under different contact pressure.

3.3 Statistical analysis of thermal conductivity

We conducted a statistical analysis of the thermal conductivity distribution within a representative 200 μm × 100 μm area, spanning from 80 μm to 180 μm along the y-axis of Figure 3a. We have excluded the pore-containing region such that the thermal conductivity histogram reflects its intrinsic statistical distribution. As shown in Figure 5a, the probability distribution exhibits an asymmetric profile that deviates from the standard Gaussian form. The measured average thermal conductivity is 13.85 W m^-1 K^-1, which is in agreement with reported values^[49] and manufacturer specifications. To further quantify the spatial heterogeneity of the thermal conductivity across the sampled domain, a statistical analysis of the thermal conductivity distribution within the sub-regions is presented in Figure 5b. Results indicate that the H1 region exhibits a higher average conductivity than the global mean due to the locally dense filler distribution, while the L2 region shows diminished values owing to sparse filler density; both regions display a positively skewed distribution, featuring a distinct long tail.

Display Full Size

Download

Figure 5. (a) Overall probability density distribution of thermal conductivity in TIMs and (b) statistical distributions of thermal conductivity across representative sub-regions and the corresponding statistical fitting using Eq. (15); (c) Schematic diagram of the Gaussian convolution model used to determine the effective thermal conductivity; (d) PDFs of the effective thermal conductivity calculated under different stochastic effects. TIMs: thermal interface materials; PDFs: probability density functions.

To explain the origin of the thermal conductivity distribution, we developed a statistical model. The measurement process of FDTR is modeled as a convolution of the ground truth thermal conductivity with a multiplicative superposition of two Gaussian expressions, as illustrated in Figure 5c. The measured effective thermal conductivity, k, is expressed as the integral of the local intrinsic thermal conductivity distribution, k_loc(r), weighted by the effective beam intensity distribution, P(r), over the entire plane:

(8)$k={\iint }_{R^2}{k}_{\mathrm{loc}}\left(r\right)\cdot P(r-r_0)dr$

where r₀ denotes the center of the laser spot. The normalized intensity distribution of the heating and probe beams, P(r), is given by $P\left(r\right)=\frac{2}{\pi w_0^2}\exp (-\frac{2r^2}{w_0^2})$. This effective radius w₀ is defined as the root-mean-square of the pump and probe beam radii, given by $\sqrt{(w_{\text{pump}}^2+w_{\text{probe}}^2)/2}$, where w_pump and w_probe are the 1/e² radii of the pump and probe laser spot respectively. SEM characterization reveals that the TIM contains Al₂O₃ fillers with an average diameter of ~3 μm, while the ZnO nanoparticles are approximately 100 nm in size. No significant heterogeneity in ZnO nanoparticle distribution was observed in SEM images (Figure S1). Given that the dimension of the ZnO nanoparticles is significantly smaller than the characteristic measurement spot size, the PDMS/ZnO mixture is treated as a homogeneous effective matrix with a thermal conductivity of k_m = 1 W m^-1 K^-1[52]. The micron-sized Al₂O₃ particles are modeled as a discrete dispersed phase with a thermal conductivity of k_p. Consequently, the local intrinsic thermal conductivity distribution is formulated as $k_{\mathrm{loc}}\left(r\right)=k_m+(k_p-k_m)\sum _i\Pi \left(\frac{r-r_i}{R_i^p}\right)$, where r_i and $R_i^p$ represent the center coordinates and radius of the i-th particle, respectively, $\Pi \left(x\right)$ is the rectangular function, $\Pi \left(x\right)=1$ for all x ∈ [-1,1] and $\Pi \left(x\right)=0$ otherwise. Due to the non-close-packed arrangement and comparable scales of particle size and beam diameter, we assume that the measured conductivity is primarily influenced by the nearest-neighbor particle^[53]. Under this approximation, the local conductivity distribution simplifies to $k_{\mathrm{loc}}\left(r\right)=k_m+(k_p-k_m)\Pi \left(\frac{r-r_p}{R_p}\right)$, where r_p and R_p represent the center coordinates and radius of the nearest-neighbor particle, respectively. Evaluating the integral in Eq. (8) in polar coordinates, we obtain:

(9)$k=k_m+(k_p-k_m)H({\rho }_0,R_p)$

(10)$H({\rho }_0,R_p)=\frac{4}{w_0^2}\exp (-\frac{2{\rho }_0^2}{w_0^2}){\int }_0^{R_p}\rho \exp (-\frac{2{\rho }^2}{w_0^2})I_0\left(\frac{4\rho {\rho }_0}{w_0^2}\right)d\rho$

where ρ = |r - r_p| and ρ₀ = |r₀ - r_p| are distances relative to the particle center, and H(ρ₀,R_p) represents the geometric overlap factor between the particle and the laser spot. The rigorous evaluation of H(ρ₀,R_p) involves a mathematically intractable integral of Bessel functions I₀. To further simplify the numerical behavior, the Gaussian Aperture Approximation, a method widely adopted in fields such as Particle Image Velocimetry, is employed to model the hard-edge particle as an equivalent soft-edge Gaussian distribution^[54,55]. This approximation effectively captures the primary features of the particle-beam interaction while simplifying the governing equations into an analytical solution as follows.

(11)$H({\rho }_0,R_p)\approx [1-\exp (-\frac{2R_p^2}{w_0^2})]\exp (-\frac{2{\rho }_0^2}{w_0^2+R_p^2})$

We first consider the stochastic spatial distribution of particles for a given size. Assuming the particle locations follow a homogeneous Poisson point process, the probability density function (PDF) of the nearest-neighbor distance ρ₀ in a 2D homogeneous Poisson point process^[56] follows a Rayleigh distribution and is derived as:

(12)$f_P\left({\rho }_0\right)=2\pi \lambda {\rho }_0\exp (-\pi \lambda {\rho }_0{}^2)$

where λ is spatial number density of the particles. Applying the probability density transformation rule $f_K(k\mid R_p)=f_P\left({\rho }_0\right)\left|\frac{d{\rho }_0}{dk({\rho }_0,R_p)}\right|$, we derive the conditional PDF of the thermal conductivity for a fixed particle size as follows.

(13)$f_K(k\mid R_p)=\pi \lambda \frac{w_0^2+R_p^2}{2(k-k_m)}{\left\{\frac{(k_p-k_m)[1-\exp (-\frac{2R_p^2}{w_0^2})]}{k-k_m}\right\}}^{-\pi \lambda \left(\frac{w_0^2+R_p^2}{2}\right)}$

Eq.(13) maps the stochastic spatial distribution of the fillers into the observed fluctuations in the local thermal conductivity distribution. This distribution arises from the combined effects of the nearest-neighbor distance distribution and the relative size between the laser spot and the thermally conductive particles. As described by Eq. (12), the nearest-neighbor distance follows a Rayleigh distribution, where the probability density first increases rapidly and then gradually decreases to zero. Since the fillers are not dimensionless points but possess finite physical radii, the effective local thermal conductivity measured by the laser is fundamentally governed by the geometric overlap area between a given particle and the laser spot. Consequently, for a fixed particle radius, the nonlinear mathematical relationship between the nearest-neighbor distance and this overlap area transforms the initial Rayleigh distribution of distances into a power-law distribution for local thermal conductivity, where the probability density gradually increases from the theoretical minimum thermal conductivity, k_m, to the theoretical maximum $k_{max}=k_m+(k_p-k_m)[1-\exp (-\frac{2R_p^2}{w_0^2})]$. In practice, the particle radii exhibit a log-normal distribution^[57]:

(14)$f_R\left(R_p\right)=\frac{1}{R_p{\sigma }_r\sqrt{2\pi }}\exp [-\frac{{(\ln R_p-{\mu }_r)}^2}{2{\sigma }_r^2}]$

where μ_r is the mean of lnR_p, and ${\sigma }_r^2$ is the variance of lnR_p. By applying the law of total probability for continuous random variables, the final PDF of the measured effective thermal conductivity is determined by integrating the conditional PDF over the particle size distribution.

(15)$f_K\left(k\right)={\int }_0^{\infty }{f}_K(k\mid R_p)f_R\left(R_p\right)d{R}_p$

Eq.(15) represents a weighted superposition of the power-law responses associated with each specific filler size, where the contribution of each particle radius is governed by its corresponding log-normal probability. This formulation elucidates that the resulting thermal conductivity distribution is essentially governed by the convolution of the intrinsic size distribution of fillers and the laser-particle distance decay effect. As shown in Figure 5a,b, we fit the distribution data using the proposed model, which reasonably captures the asymmetric nature of the observed thermal conductivity distribution. Further details of the model fitting are provided in Supplementary materials.

Given the experimental parameters (k_m = 1 W m^-1 K^-1, k_p = 30 W m^-1 K^-1, w₀ = 1.6 μm, μ_r = 0.1 μm, σ_r = 0.2 μm, λ = 0.4 μm^-2), numerical integration is performed with the upper limit of integration truncated at 10 μm. The final PDF of the measured effective thermal conductivity is highly consistent with the asymmetric probability distribution observed experimentally. Notably, by decoupling the spatial and size-related contributions within this statistical framework, we elucidate that the spatial stochasticity introduces a characteristic non-linear decay effect, whereas the randomness in particle size primarily dictates the log-normal skewness of the profile. The PDFs of the effective thermal conductivity calculated under varied stochastic conditions are compared in Figure 5d. From the comparison of the three models, it is evident that the log-normal-like characteristic of the experimentally measured thermal conductivity originates from the log-normal distribution of particle sizes within the finite measurement volume. Meanwhile, the spatial randomness introduces a non-linear decay effect, ultimately shaping the result into a non-Gaussian distribution with distinct heavy-tail characteristics. Similar long-tailed asymmetrical thermal conductivity distribution is also observed in an elastic thermal pad, see Section S6.4. Our model provides an explicit causal link between the particle size distribution and the asymmetric thermal conductivity histogram.

4. Conclusion

In summary, we demonstrated FDTR thermal-property imaging of TIMs at the microscale, enabling simultaneous mapping of local thermal conductivity and interfacial thermal conductance in a sandwich-configured sample fixture. The thermal-property images reveal pronounced spatial heterogeneity in a commercial Al₂O₃/ZnO-loaded thermal grease. Co-located micro-CT imaging associates these features with the underlying filler architecture: high-conductivity regions correspond to filler aggregation, whereas low-conductivity regions correlate with pores or filler-deficient structures. Pressure-dependent measurements further demonstrate that while the averaged thermal properties remain relatively stable over a moderate pressure range, the thermal-property maps capture local redistribution of conductive features during loading and unloading. We develop a statistical model showing the link between the particle size distribution and the statistical behavior of local thermal conductivity. The combined FDTR–micro-CT approach presented in this work provides a visualization tool for characterizing the structure-property correspondence at the microscale, which can facilitate future development of advanced thermal management materials.

Supplementary materials

The supplementary materials for this article is available at: Supplementary materials.

Authors contribution

Yao Y: Writing-original draft, visualization, validation, software, methodology, investigation, formal analysis, data curation.

Yang H: Writing-review & editing, validation, methodology, investigation.

Yang R: Conceptualization, writing-review & editing, project administration, funding acquisition.

Qian X: Conceptualization, writing-review & editing, methodology, formal analysis, supervision.

Conflicts of interest

Ronggui Yang is an Editorial Board Member of Thermo-X. The other authors declare no conflicts of interest.

Ethical approval

Not applicable.

Consent to participate

Not applicable.

Consent for publication

Not applicable.

Availability of data and materials

Data supporting the findings of this study are available from the corresponding author upon reasonable request.

Funding

This work was supported in part by the National Key R & D Project from the Ministry of Science and Technology of China (Grant No. 2022YFA1203100).

Copyright

© The Author(s) 2026.