Application of Scattering Parameters to DPL Time-Lag Parameter Estimation at Nanoscale in Modern Integration Circuit Structures

: This paper presents the methodology of material parameters’ estimation for the dual-phase-lag (DPL) model at the nanoscale in modern integration circuit (IC) structures. The analyses and measurements performed were used in the unique dedicated micro-electro-mechanical system (MEMS) test structure. The electric and thermal domain of this structure was analysed. Finally, the silicon dioxide (SiO 2 ) temperature time-lag estimation procedure is presented based on the scattering parameters measured by a vector network analyser for the considered MEMS structure together with the 2-omega method. The proposed methodology has the ability to estimate the time-lag parameter with high accuracy and is also suitable for the temperature time-lag estimation for other manufacturing process technologies of ICs and other insulation materials used for integrated circuits such as silicon nitride (Si 3 N 4 ), titanium nitride (TiN), and hafnium dioxide (HfO 2 ).


Introduction
The modelling and analysis of thermal and electromagnetic phenomena in integrated circuits (ICs) have great significance for the reliable design of modern nanoscale electronic structures and micro-electro-mechanical systems (MEMSs). It is also highly recommended to use this approach during various design-related activities. The most important of them are as follows: continuous decrease of semiconductor technology nodes, design cost reduction of microelectronic systems, and application of modern transistors in ICs (e.g., a fin field-effect transistor-FinFET [1,2], a gate-all-around field-effect transistor-GAAFET [3], and a vertical-slit field-effect transistor-VeSFET [4]).
Therefore, the presented work focuses on the following issues: • The effective modelling and simulation methods of electromagnetic (EM) and thermal phenomena of integration circuits, micro-electro-mechanical systems (MEMSs), and 3D integrated systems (see [1,5]).

•
The design of an equivalent electric circuit extractor and electrical simulation tools, e.g., [2,6,7]. • The verification of new simulation tools dedicated for the design of experimental application specific integrated circuits (ASICs). • Estimation of IC material parameters for the nanometric technologies.
The first and last issues will be presented in the paper and are performed on scattering parameters measurements of a designed MEMS test structure with a printed circuit board (PCB) test board.
for x ∈ R n , n ∈ N, t ∈ R + ∪ {0} (1) where q is heat flux vector; k is material thermal conductivity; T is distributed temperature; and c v and q V are volumetric heat capacity and value of internally generated heat, respectively. The FK equation has been commonly used to describe heat flow for over two centuries. However, the application of the FK equation has some significant limitations at the nanoscale [10].
The main limitation is related to the assumption of the infinite heat flow velocity in the semiconductor structure, which is very important in the case of nanostructures with a Knudsen number of Kn < 0.1 (where the characteristic length of IC is ca. ten times longer than the heat carrier mean free path 41.8 nm at 300 K [8]). Furthermore, the simultaneous change of temperature gradient and heat flux is a non-physical behaviour [9][10][11][12][13]; these phenomena have been presented in the case of nanosized structures in [8]. The second problem is associated with the radiative heat transfer and the Casimir effect (i.e., photons tunnelling) [14,15]; these phenomena occur in MEMSs [16] at small distances between surfaces (or particles). Both phenomena cause more dynamic and intense heat transfer than the classical Fourier-Kirchhoff and radiation models.
The new and more general dual-phase-lag (DPL) model of heat flow in solids at the nanoscale was introduced by Tzou [17,18] and is described with the following equations: where τ T and τ q are temperature time-lag and heat flux time-lag, respectively. The simplified version of the DPL heat transfer model can be obtained for one time-lag parameter τ q for τ q > 0 and τ T = 0. This equation was proposed by Cattaneo and Vernotte (CV) in [19][20][21][22][23][24][25][26] and is also used to eliminate the paradox of the heat-flux infinite velocity transfer. The comprehensive presentation of the DPL model aspects is included in the following published sources:  [22][23][24]30].
The DPL model can be used to describe heat transfer behaviour in various solid materials (e.g., metals, bulk crystals, dielectrics, and amorphous materials). For most of them, indicative dual-phase time-lag parameters values (τ T , τ q ) can be approximated using theoretical formulas [1]. Unfortunately, these parameters are strictly dependent on material density, fractal dimensions of material structure [30], and so on. Therefore, the values of the time-lag parameters should be empirically and independently characterised for each IC and MEMS technology process types. It is possible to perform empirical measurements of their values using a unique tool based on the femtosecond pulse laser with reflection mirrors for the following:
However, the test probes used for this purpose should be prepared not exactly like real structures in ICs and MEMSs. The first investigation of the dual-phase time-lag parameters estimation using dedicated MEMS structures with a measurement system implemented on a PCB was depicted in [39][40][41][42][43]. The authors of this paper propose the scattering wave parameters presented in [2] to extract the value of the dual-phase time-lag parameter τ T .

Test Structure
The DPL model was verified using a dedicated test structure containing a 104 nm thin silicon oxide layer (SiO 2 : ε r = 3.9 [44,45], we assume a constant value up-to 13.5 GHz based on [46]; tan δ = 0.0002 [47]) placed between two platinum resistors (the top resistor: Pt: W top = 2.75 µm, thickness h = 150 nm, length L = 460 µm, resistance R top = 155 Ω at 20 • C [39,48,49]; the bottom resistor: Pt, W bottom = 4.96 µm, thickness h = 150 nm, length L = 460 µm, resistance R bottom = 262.3 Ω at 20 • C [39,48,49]). Generally, each resistor could be used either as a heat source line or a temperature sensor. The bottom resistor is treated as a heater, and the upper resistor is a temperature sensor. Each resistor has a four-wire connection for independent measurements of voltage and current. These elements are placed on a 100 nm wide silicon dioxide (SiO 2 ) layer stacked on a 400 µm thick silicon layer. The test structure presented in Figure 1c was manufactured in MEMS technology by the Institute of Electron Technology (ITE) in Warsaw, Poland [39][40][41][42][43][44][45]. The bottom and upper resistor dependence on temperature and their temperature rise dependence on heating power were measured and characterised in [39,40]. The entire structure is bonded to a metal copper foil on a PCB test board (Figure 1a

Thermal Model
The 2D DPL thermal model of the MEMS test structure is described with the application of the harmonic temperature models proposed in [1]:

Thermal Model
The 2D DPL thermal model of the MEMS test structure is described with the application of the harmonic temperature models proposed in [1]: where T = T(x, y, z) = T r + jT i ∈ Z is harmonic temperature field distribution for angular frequency ω = 2πf, associated with the temperature in the following way: Moreover, T m and ϕ are amplitude and phase-shift for a given angular frequency, respectively. Please note that the maximal complex values of the harmonic temperature field are taken into consideration; thus, the coefficient √ 2 is introduced in Equations (4)-(7). The dissipated heat power generation per unit volume in the bottom resistor (g 0 > 0) is described in Appendix A. The heat power generation outside the bottom resistor is equal to zero g 0 = 0 W/m 3 . The following boundary conditions were applied: • The upper side of the structure-the heat-free convection is encoded using Neumann boundary condition with a heat transfer coefficient equal to 658.763 W/(m 2 K).

•
The left and right side of the structure-the symmetry is applied using Neumann boundary conditions. • The Dirichlet boundary condition for the bottom side of the structure T r = T ambient , The material parameters are included in Table 1. In addition, calculated temperature distributions for 10 Hz, 1 kHz, 100 kHz, and 1 GHz are presented in Table 2.
As can be observed in Figure 2, the changes of temperature amplitude |T(x,y,f )| in SiO 2 are located around the bottom heating resistor at f = 100 MHz and higher frequencies. Therefore, SiO 2 thermal characterisation should be performed for f ≥ 100 MHz. Consequently, the simulations were carried out for the frequency in the range f = 200-900 MHz.
Energies 2021, 14, 4425 6 of 14 want to propose a methodology of the DPL time lag estimation for other insulation materials applied in the integrated circuits, such as silicon nitride (Si3N4), titanium nitride (TiN), hafnium dioxide (HfO2), and others. Thus, numerous simulations were carried out for the temperature time-lag-delay parameter in the range τT = 50-120 ns to find the estimated value. The experiments and analyses carried out by Tomasz Raszkowski in his PhD thesis [50,51], concerning the transformation of the DPL equation into the Grünwald-Letnikov fractional space-derivative dual-phase-lag equation [16,46,47], indicate a lower influence of the heat flux time lag parameter τq for the whole thermal system response than for the temperature time lag parameter τT.
The analysis of the influence of time lag differences τT-τq [1] (pp. 140- Figure 2, pp. 143- Figures 5 and 6, pp. 142- Table 2) as well as the total thermal delay of complex IC thermal system [1] (pp. 151, see Table 6) show that the heat-flux time lag parameter is less critical than the temperature-time lag for final thermal system response. Therefore, we adopted its value equal to τq ≈ 3 ns, as proposed in the literature. Simulations were carried out and are presented in Figure 3. It can be noticed that the most significant change in the phase of temperature depending on the time-lags is recognised at the frequency range near 300 MHz and 500 MHz. However, the most significant amplitude of temperature changes on the upper resistor is observed close to the lower frequency f ≈ 300 MHz (see Values of a temperature time-lag and a heat flux time-lag for bulk SiO 2 material have been proposed in the literature equal to τ T ≈ 60 ns and τ q ≈ 3 ns, respectively (see [1,30]). Nevertheless, we also have to determine this value for the ITE process technologies. We want to propose a methodology of the DPL time lag estimation for other insulation materials applied in the integrated circuits, such as silicon nitride (Si 3 N 4 ), titanium nitride (TiN), hafnium dioxide (HfO 2 ), and others. Thus, numerous simulations were carried out for the temperature time-lag-delay parameter in the range τ T = 50-120 ns to find the estimated value.
The experiments and analyses carried out by Tomasz Raszkowski in his PhD thesis [50,51], concerning the transformation of the DPL equation into the Grünwald-Letnikov fractional space-derivative dual-phase-lag equation [16,46,47], indicate a lower influence of the heat flux time lag parameter τ q for the whole thermal system response than for the temperature time lag parameter τ T .
The analysis of the influence of time lag differences τ T -τ q [1] (pp. 140- Figure 2, pp. 143- Figures 5 and 6, pp. 142- Table 2) as well as the total thermal delay of complex IC thermal system [1] (pp. 151, see Table 6) show that the heat-flux time lag parameter is less critical than the temperature-time lag for final thermal system response. Therefore, we adopted its value equal to τ q ≈ 3 ns, as proposed in the literature. Simulations were carried out and are presented in Figure 3. It can be noticed that the most significant change in the phase of temperature depending on the time-lags is recognised at the frequency range near 300 MHz and 500 MHz. However, the most significant amplitude of temperature changes on the upper resistor is observed close to the lower frequency f ≈ 300 MHz (see Figure 3a). Therefore, the highest sensitivity of time-lag estimation is at f ≈ 300 MHz. Finally, f ≈ 300 MHz frequency neighbourhood, in the range of τ T = 50-120 ns and τ q ≈ 3 ns, was selected as the best parameter value set for further investigation (Figure 4b).
Energies 2021, 14, x FOR PEER REVIEW 7 of 15 Figure 3a). Therefore, the highest sensitivity of time-lag estimation is at f ≈ 300 MHz. Finally, f ≈ 300 MHz frequency neighbourhood, in the range of τT = 50-120 ns and τq ≈ 3 ns, was selected as the best parameter value set for further investigation (Figure 4b).
(a) (b)    Next, based on the collected simulation results, the linear regression was used to model the relationship between the temperature time-lag τT and the averaged phase-shift of the upper resistor temperature arg T(x,y,f) is used for the analysed MEMS test structure, which is presented in Figure 1c (τT  (50 ns, 120 ns) and τq = 3 ns) with R 2 = 0.999774.
 T (arg(T upper (f))) Next, based on the collected simulation results, the linear regression was used to model the relationship between the temperature time-lag τ T and the averaged phase-shift of the upper resistor temperature arg T(x,y,f ). The function is used for the analysed MEMS test structure, which is presented in Figure 1c (τ T ∈ (50 ns, 120 ns) and τ q = 3 ns) with R 2 = 0.999774.

Electrical Analysis and Final Measurement Procedure
The electrical behaviour or the MEMS test structure (Figure 1c) was determined using the Vector Network Analyzer PNA-X Microwave Agilent N5242A after application of the de-embedding procedure [2]. The electrical reference planes were set to the test board SMA (SubMiniature version A) connectors using Agilent N4691B Electronic Calibration Module. All non-utilised connectors were terminated. The complex scattering parameters (s 11 , s 12 , s 21 , s 22 ) have to be measured at 32,000 spot frequency points in the range 10 kHz-15 GHz using 50 Ω characteristic impedance of the measurement equipment.
The analysis of the electric behaviour of the MEMS test structure (Figure 1d) and the PCB test board is based on the previous study about transmission lines presented in [2]. Therefore, the following conclusions were derived for the chosen frequency close to f ≈ 300 MHz: • The upper and lower resistor lengths are very short compared with the wavelength at f ≈ 300 MHz. Therefore, the distributed character of both resistors shall be taken into account for f > c/(20 L ε r −1/2 ) ≈ 64 GHz, where c is the speed of the light in the vacuum, and the relative electric permittivity of SiO 2 is equal to ε r ≈ 3.9. Consequently, only a de-embedded procedure for a PCB is required for the PCB transmission lines' elimination [2].

•
The investigation presented in [2] shows that the eddy-current losses (and skin effect) in a conductor and dielectric losses can be neglected with a 1.4% error up to 10 GHz and an error lower than 0.3% up to 2 GHz. • The extracted parasitic parameters in Table 2 show that inductance and capacitance can be neglected for f ≈ 300 MHz. • It should also be noted the self-heating can be neglected for a small current flowing through the platinum resistor (e.g., I bottom ≤ 5 mA [44]).

•
As can be noticed, the application of the network analyser allows for the direct measurement of the power wave transmitted through the MEMS test structure s 12 = s 21 , which represents the insertion loss of the analysed network. To determine the temperature time-lag delay, we used the s 21 measurements for a frequency close to f = 150 MHz (n = 21 spot frequency points in the range 145 MHz ≤ f ≤ 155 MHz) and determined the averaged value and standard deviation of this parameter arg(s 21 ) = 62.3219 (±0.676051) • at f = 150 MHz, with the sum of error squares (SS) = 9.01973 and degrees of freedom (DOFs) = 21 (see Figure 5). Finally, the temperature time-lag delay can be calculated from Equation (8): τ T = 61.9791(±1.61694) ns at f = 300 MHz (9)

Results and Discussions
As can be observed in Figure 2, the changes of temperature amplitude |T(x,y,f )| in SiO 2 are localised around the bottom heating resistor at f = 100 MHz and higher frequencies. Therefore, SiO 2 thermal characterisation should be performed for f ≥ 100 MHz.
Based on the simulation, it can be noticed that the best scopes for which we obtain the highest sensitivity of the measurement method are the frequency range near 300 MHz to 500 MHz (see Figure 4). However, the most significant amplitude of temperature changes on the upper resistor is observed close to the lower frequency f ≈ 300 MHz (Figure 3a).
Next, based on the collected simulation results, the linear regression was used to model the relationship between the temperature time-lag τ T and the averaged phase-shift of the upper resistor temperature arg T(x,y,f ) in the form of a function: τ T | f = 300MHz = 2.391746 (±0.064775)·arg T upper − 87.07922(±4.52158) [ns, • ] (10) for the analysed MEMS test structure, which is presented in Figure 1c (τ T ∈ (50 ns, 120 ns) and τ q = 3 ns) with R 2 = 0.999774.

Results and Discussions
As can be observed in Figure 2, the changes of temperature amplitude |T(x,y,f)| in SiO2 are localised around the bottom heating resistor at f = 100 MHz and higher frequencies. Therefore, SiO2 thermal characterisation should be performed for f ≥ 100 MHz.
Based on the simulation, it can be noticed that the best scopes for which we obtain the highest sensitivity of the measurement method are the frequency range near 300 MHz to 500 MHz (see Figure 4). However, the most significant amplitude of temperature To determine the temperature time-lag delay, we used the s 21 measurements for a frequency close to f = 150 MHz and determined the averaged value and standard deviation of this parameter arg(s 21 ) = 62.3219 (±0.676051) • at f = 150 MHz (see Figure 4). Finally, the temperature time-lag delay can be calculated from Equations (8) and (10)  τ T = 61.9791(±1.61694) ns at f = 300 MHz (11) for I top ≤ 5 mA and I bottom ≤ 5 mA [44].

Conclusions
This paper includes electromagnetic and thermal analyses of the MEMS test structure developed to validate the DPL model based on the methodology and tools presented in [2] for electric and thermal domains [1], respectively. The detailed analysis allows determining the best measurement frequency as well as electric and thermal model simplifications. Moreover, the new approach is associated with conducting measurements of electromagnetic scattering parameters together with simplification of the EM and thermal models of the MEMS test structure.
As a final result, we used a network analyser to determine the temperature time-lag delay in SiO 2 τ T = 61.9791 (±1.61694) ns for assumed heat flux time-lag delay τ q = 3ns, which is a big challenge in the case of ICs and MEMSs. The final time-lag delay values are applicable for the devices manufactured in the technology process developed by the Institute of Electron Technology in Warsaw, Poland.
The proposed method is suitable for determining the time-lag delay parameters of other materials and other technological processes of ICs, which can be very important during thermal and electro-thermal investigations of advanced electronic nano structures. The ability of high accuracy estimation of temperature distribution in the newest electronic devices is very significant for the real thermal model analysis. As a result, it is crucial for the reliability of whole nanoscale appliances in modern semiconductor structures.  Acknowledgments: The authors would like to express their special thanks to M. Janicki and J. Topilko for sharing papers [39,40] and [43][44][45].

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Analysis of Dissipated Power and Internal Heat Generation Source in MEMS Test Structure
With a limited error, we can assume that the bottom resistor of the considered MEMS test structure is a lamped resistive element for the electric domain (L << λ/20, see Section 4.2 and Table 2). The power dissipated in this resistor is equal to p(t) = u bottom (t)·i bottom (t), where i bottom (t) = A·sin(2πf E t) is a current flowing through this resistor; A is the amplitude of current, e.g., A = (5/ √ 2)mA; u bottom (t) = |Z bottom |·sin(2πf E t + φ) is a voltage drop on the resistor, φ = argZ bottom ≈ 0.12 • at 150 MHz, |Z bottom | ≈ R bottom at 150 MHz; and f E is a frequency in the electric domain (e.g., 150 MHz).
In this appendix, we will use the Fourier transform of x(t) in the limit sense defined according to the following equation [1]: where δ(t) is the Dirac delta function; ω 0 = 2πf E is pulsation, T = 1/f E ; and X k is a set of complex Fourier coefficients given by The series of expansion coefficients of the periodic function x(t) into the complex series {Xn} will be defined as x(t) ↔ {X k } for k = −∞, . . . ,0, . . . ,+∞. The current flowing through the resistor has two complex coefficients associated with two alternating currents (ACs) ±f E in the electric frequency domain for k = ±1. where V bottom = L·W bottom ·t bottom ; |Z bottom |≈R bottom and Re(e ±j2φ ) = cos(2φ) ≈ 1 for φ ≈ 0.12 • at 150 MHz; the expression g bottom,k / √ 2 is used as the last left term in Equation (4). It should be underlined that the heat generated per unit in the bottom resistor is associated with the direct current (DC, for k = 0) and two alternating currents ±2f E in the thermal frequency domain for k = ±2. This is why the harmonic temperature field distribution (4) is performed at a frequency of, e.g., 300 MHz, while scattering parameters measurements (s 21 ) are determined at half the frequency of 150 MHz. It should also be noted that detailed information on the operational principles of the vector network analyser can be found in [52]. The application of ±f E in the electric domain and ±2f E in the thermal domain was presented for the first time in [53].