Estimation of Thermal Resistance Field in Layered Materials by Analytical Asymptotic Method

: In this paper, the problem of the quantitative characterization of thermal resistance ﬁelds in a multilayer sample is addressed by using the classical front face ﬂash method as the thermal excitation and infrared thermography (IRT) as the monitoring sensor. In this challenging problem, the complete inverse processing of a multilayer analytical model is difﬁcult due to the lack of sensitivity of some parameters (layer thickness, depth of thermal resistance, etc.) and the expansive computational iterative processing. For these reasons, the proposed strategy is to use a simple multilayer problem where only one resistive layer is estimated. Moreover, to simplify the inverse processing often based on iterative methods, an asymptotic development method is proposed here. Regarding the thermal signal reconstruction (TSR) methods, the drawback of these methods is the inability to be quantitative. To overcome this problem, the method incorporates a calibration process originating from the complete analytical quadrupole solution to the thermal problem. This analytical knowledge allows self-calibration of the asymptotic method. From this calibration, the quantitative thermal resistance ﬁeld of a sample can be retrieved with a reasonable accuracy lower than 5%.


Introduction
In the domain of material and microelectronic device characterizations, many groups and studies have emerged with the goal to estimate both the thermal properties and resistance of these thin multilayers. One of the most common procedures is based on the 3ω technique [1], which has been extensively used since it offers absolute measurements of the heat flux and the temperature and is also well suited for characterizations at low temperatures. At high temperatures, contactless photothermal methods, such as the thermoreflectance [2,3] and infrared radiometry [4][5][6] techniques, have been implemented and allow measuring the relative change in temperature and heat flux. Nevertheless, calibration is a complex task within those experimental configurations. In this community, very few studies are based on field characterizations. The notion of fields introduces the requirement to propose new methods for which the computation time is reasonable, passing from only one time vector, to find the inverse of a matrix. On the other hand, within the research groups involved in quantitative infrared thermography, many authors have developed a technique for qualitative in-depth defect detection by using flash thermography and thermal signal reconstruction (TSR) [7,8] or the lock-in technique [9]. Few authors, e.g., Tranta et al. [10] and Roche et al. [11], have linked the TSR method with quantitative estimation of the defect depth. The same kind of studies have been performed by Ibarra et al. [12] with a technique based on pulsed phase transformation (PPT) and thermographic signal reconstruction similar to [13]. As depicted in [14] one alternative to these qualitative methods is the use of parametrized 1D simulations with the quadrupole method as presented in the book of Maillet et al. [15]. Some works were performed by coupling infrared thermography (IRT) and the quadrupole approach, such as the one of Bendada et al. [16][17][18]. In fact, with these analytical methods, it becomes easy to calculate a multilayer sample with several thermal properties as well as thermal resistance. Then, by using nonlinear regression methods as proposed by Muller et al. [14], the parameters of the model can be adjusted by least square minimization between the model output and the experimental data. This approach is difficult when the initial values are poorly known and could be time consuming when addressing a large amount of data, such as temperature fields, with IRT. In the publication of Feuillet et al. [19], artificial delaminations were manufactured by introducing PTFE (Teflon) inlays in carbon-fiber-reinforced plastics, and they attempted to reduce the inversion by using singular value decomposition (SVD) [20][21][22][23]. Nevertheless, in this challenging problem, the complete inverse processing of a multilayer analytical model is difficult due to the lack of sensitivity of some parameters (layer thickness, depth of thermal resistance, etc.) and the expansive computational iterative processing. For these reasons, the proposed strategy is to use a simple multilayer problem where only one resistive layer is estimated. Moreover, to simplify the inverse processing often based on iterative methods, an asymptotic development method is proposed. The main novelty comes from the calibrating method between the analytical calculation of a theoretical base and the estimated parameters of the asymptotic method. Then with the knowledge of this calibration curve, the estimated parameters based on the measured one can be turn to quantitative ones.
Finally, in this work, a well-known flash setup is presented. A method is developed to model the heat transfer within the multilayer sample by using the quadrupole formulation. Then, in the inverse processing section, the autoregressive asymptotic development method is presented as well as the calibration procedure. Finally, the results and discussion section presents the characterization of a standard sample as a validation of the method.

Experimental Setup and Materials
The setup illustrated in Figure 1a and presented in Figure 1b is the classical front face flash system associated with an IR camera. The flash lamp comes from Uniblitz and has an energy E of 1600 J. Each lamp is synchronized with the IR camera by using an analogical TTL link to perform a pretrigger mode before the flash. The IR camera is an MCT long-wave (λ = 9-11 µm) FLIR SC7600 with a matrix sensor of 320 × 256 pixels and a pitch of 30 × 30 µm 2 . With the used lenses and the sample distance, the resulting spatial resolution is approximately 300 × 300 µm 2 per pixel.
For this study, a homemade reference sample was machined to calibrate the depth and thickness of the thermal resistance, represented by an air square hole, as illustrated in Figure 2. The sample is composed of polycarbonate layer of thickness e PC and the air thickness of the thermal resistance is noted e Rth .The geometrical and thermal properties are reported in Figure 2 and extracted from reference [24] for aluminum and [25] for the polycarbonate. In this simple case, only a bilayer sample was designed to demonstrate the capacity of the proposed method. However, the method can be extended to multilayer stack. It is important to note that the minimum tolerance guaranteed by the supplier is approximately +50 µm of the desired thickness (see Figure 2).

Direct Problem
As illustrated in figure 3, the direct problem of the heat transfer within a multilayer sample can be solved by using the quadrupole method [15]. It is important to note that according to the multilayer used for characterization, the complete stack can be overlaid onto a substrate with different thermal properties to also control the boundary condition imposed at the rear face of

Direct Problem
As illustrated in Figure 3, the direct problem of the heat transfer within a multilayer sample can be solved by using the quadrupole method [15]. It is important to note that according to the multilayer used for characterization, the complete stack can be overlaid onto a substrate with different thermal properties to also control the boundary condition imposed at the rear face of the sample. Usually, a foam bulk is used to mimic an adiabatic condition, whereas a metallic plate is used to mimic the imposed temperature. Using the quadrupole formulation [15], the thermal response of any assembly (see Figure 3) can be generalized according to the Equation (1). The quadrupole method enables taking into account many layers as well as thermal resistances. A MATLAB homemade application (https://fr.mathworks. com/matlabcentral/fileexchange/74199-1dt-multilayer-thermalquadrupole-solver-and-builder) was developed to generate any assembly and to realize a sensitivity study for any parameters of the multilayer. The model is then expressed as: In the Equation (1), M i is the matrix relative to layer i and R i is the matrix associated with the thermal resistance between layers i and i + 1. They are both expressed in Equation (2). Moreover, θ e (K) and θ s (K) represent the Laplace transforms of the front and rear face temperatures, respectively. Ψ e (p)(W.m −2 ) is the Laplace transform of the heat flux at the front face. Usually, the heat flux waveforms can be Dirac, heaviside or a pulse with a duration linked with the flash lamp discharge time. Finally, h(W.m −2 .K −1 ) represents the convective heat losses.
In Equation (2) Since the heat flux amplitude is unknown, the data are normalized according to any time step included in the semi-infinite behavior of the temperature relaxation inside the first cross layer according to the following formulation: T n (t) = T(t ≥ t re f ) T(t re f ) . To calculate the inverse Laplace transform of the temperature, the Stehfest algorithm [26] is used. The dimensionless temperatures as a function of time and for different polycarbonate layer thicknesses and thermal resistances are plotted in Figure 4. The thermal properties of the polycarbonate layer and the aluminum layer are given in Figure 2.

Autoregressive Asymptotic Method
The asymptotic method consists of the development of the amplitude and time delay between a reference temperature signal and a measured or generated one. The main assumption is that this pseudothermal contrast is sensitive only to the intensity variations A (due to the presence of the thermal resistance) and a delay τ(s) proportional to both the depth of the thermal resistance and its thickness. If both parameters are estimated only the amplitude will be used to quantitatively estimated the thermal resistance. From the asymptotic Equation (3), and to avoid the influence of noisy data, an autoregressive inverse algorithm [27] is applied for a long time, as described by the Equation (4). The inverse processing is based on the Moore-Penrose inversion method [28,29]. The advantage of using this autoregressive method is the ability to estimate the parameters as a function of time to create a strong link with the maximum sensitivity corresponding to the early "separation" between the reference signal and the signal to process. With the above remarks under consideration, the inverse of the system can be written as follows where A(−) is the amplitude of the asymptotic development, τ(s) is the time delay, t(s) is the time, and the indices m and re f correspond to the measured temperature and the reference one. From Equation (4), the pseudoinverse X is obtained by using the singular value decomposition (SVD) methods [20][21][22][23].
The principle of SVD is to decompose a matrix I of size m by n into a product of three matrices, where U is an orthogonal matrix of size m by m, V is an orthogonal matrix of size n by n, and S is a "diagonal" matrix of size m by n. The diagonal elements of S are called singular values. They are positive and ranked from largest to smallest. If only the first p singular values are different from zero (p < n), the SVD can then be further reduced. Thus, the pseudoinversion of Equation (4) gives: From the reference sample and the temperature response calculated and plotted in Figure 4, the autoregressive method is applied to all the data, with the Rth = 0 plot as the reference, corresponding to an assembly without any thermal resistance. The estimated parameters A(−) is reported in Figure 5 as a function of the generated thermal resistance Rth . The time required to calculate all the time steps of the autoregressive method is approximately 6 s with a 2.3 GHz quad-core i5 MacBook Pro with 16 GB of RAM. The data can be fitted by using polynomial regression for the thermal resistance (Figure 5b). These results can be qualified by using a calibrating curve to link the nonquantitative estimation given by the asymptotic development method to the absolute value of he thermal resistance.

Results and Discussion
In this part, the results on the sample described in Section 2 are measured under the following parameters: (i) flash duration of 1 ms, (ii) camera frequency acquisition of 200 Hz, and (iii) 4000 images corresponding to a final time of 20 s and a pixel area of 300 × 300 µm 2 for the spatial resolution. The obtained results are plotted in Figure 6: From these measurements, 5 regions of interest (ROIs) appeared. For each ROI, the corresponding temperature evolution of one pixel can be extracted: 4 of them in the center of the thermal resistance square hole and the fifth in the center of the image. All these temperatures plot are normalized (baseline subtraction and reference division T * n (t) = , where t 0 (s) is the time of the flash and t re f (s) any time during the -1/2 slope) and superposed onto the "healthy" analytical model (calculated for the bilayer of Figure 2) assuming perfect contact between all the layers. This normalization also avoid the non-uniformity of the global temperature coming from the non-homogeneity of the flash lamp. The results represented in Figure 7 show a -1/2 slope for a short time and the influence of the thermal resistance as well as the different depth effects as a function of time. In the Figure 7, it can be clearly observed that for the deepest ROI the temperature evolution quit the -1/2 slope further. For higher thermal resistance, the behavior is closer to adiabatic case with variable inflexion tendencies. Moreover, it is very interesting to note that the ROI extracted in the center of the image seems to be affected by a small thermal resistance (no exact superposition with the supposed healthy model can be achieved). From this, we can estimate in Figure 8 the corresponding parameters. First, in Figure 8a, the obtained amplitude is plotted. Here, one can see that the estimated values are in good agreement with the ones estimated using the analytical solution represented in Figure 5a. Then, from the calibrated curves calculated and represented in Figure 5b, the absolute and quantitative fields of the amplitude value can be represented, as shown in Figure 8b. Here, the time required to calculate the parameters with the inverse process (4) is approximately 7 s with a 2.3 GHz quad-core i5 MacBook Pro with 16 GB of RAM. This computational time includes the inversion of a 51,012 × 4000 matrix representative of the number of pixels in space (218 × 234) and the number of time steps of the thermograms.
In this validation case, the geometry and thermal properties of the studied sample are well known (see Figure 2). Then, from the results in Figure 8, the agreement with the value is good, and the respective theoretical and estimated values are reported in Table 1.
From the data Table 1 one can see that the error is always lower than 5%, with the maximum uncertainty for ROI number 3 corresponding to a hole of 500 µm. One can also note that a threshold has been performed before the identification in order to treat only the pertinent points belonging to the patches. Indeed, a certain number of points in the ROIs borders are not pertinent because of the strong gradient present here due to the 3D effect of the diffusion not taken into account in the identification process (see Figure 8). For example, for ROI3, 1500 points over 1900 have been used for the identification.   Figure 2 or the results in Figure 8; **, Systematic maximum machining uncertainty of +50 µm given by the supplier; see Section 2; ***, Error 100 × 1 − Rth est Rth theo calculated with the average value Rth theo of the Rth theo gap.

Conclusions
In this paper, a methodology for quantitative thermal resistance imaging estimation is proposed. This method is based on asymptotic development and quantitative calibration by using an analytical quadrupole formulation of the problem. The main advantages of this method are as follows: (i) the possibility to address any type of multilayer (high numbers of layers, etc.), (ii) a fast inverse process that can be adapted for online monitoring (inverse processing in less than 7 s), and (iii) the ability to overcome the knowledge of the initial excitation heat flux as well as the absolute measured temperature. The main drawbacks of the proposed methods are (i) the 1D character of the model and (ii) the sensitivity to only one thermal resistive layer. Nevertheless, the complete method applied on a heterogeneous bi-layer sample reveals a very good result with a complete estimation of the thermal resistance fields in fast time (complete calculation time lower than 7 s) and accuracy with maximum thermal resistance errors around 5%.