*2.1. Finite Di*ff*erence Method Approximation*

Thermal analysis was carried out for the real microelectromechanical system (MEMS) nanostructure, presented in Figure 1, which was manufactured at the Polish Institute of Electron Technology [14,15]. This test structure includes two parallel platinum resistors, each 10 μm long. One resistor is treated as a heater, whereas the second resistor plays the role of a temperature sensor. The cross-sectional area of each resistor is 100 nm wide and 100 nm high. The distance between these resistors is 100 nm. The resistors are placed on a 100 nm wide silicon dioxide layer and both are stacked on a 500 μm thick silicon layer. A detailed description of the investigated structure, as well as the measurement process are found in [14,15].

**Figure 1.** Cross-sectional area of the microelectromechanical system (MEMS) test structure where the Dirichlet boundary conditions are marked with a blue line, the Neumann boundary conditions are marked with a red dashed line, and the internal heat source is marked using red color.

In the investigated cross-sectional area of the MEMS test structure, the Dirichlet boundary conditions are at the bottom, while the Neumann boundary conditions are used on the left, right, and the top parts of the structure and their environment. The boundary conditions can be described by the following equations:

$$\begin{aligned} T\_k(t) &= 0, t \in \mathbb{R}\_+ \cup \{0\}, \\ k &\in \{1, 2, 3, \dots, n\_x\} \end{aligned} \tag{3}$$

$$\begin{array}{c} \begin{aligned} q\_k(t) &= 0, t \in \mathbb{R}\_+ \cup \{0\}, \\ k \in \left\{ n\_x + 1, 2 \cdot n\_x + 1, \dots, \left( n\_y - 1 \right) \cdot n\_x + 1 \right\} \cup \left\{ n\_{\overline{x}, 2} \cdot n\_{\overline{x}, \dots, \left( n\_y - 1 \right) \cdot n\_x \right\} \cup \\ &\cup \left\{ \left( n\_y - 1 \right) \cdot n\_x + 1, \left( n\_y - 1 \right) \cdot n\_x + 2, \dots, n\_y \cdot n\_x \right\} \end{aligned} \end{array} \tag{4}$$

where *n*x and *<sup>n</sup>*y reflects the number of discretization nodes in both axes OX and OY, respectively.

Thermal simulations for the investigated MEMS structure were performed for its two-dimensional (2D) cross-sectional area, in the middle of the resistor. The thermal analysis of the structure was carried out using the DPL model. In order to make the analysis easier, the DPL model described by the system Equation (2) was equivalently transformed to the following 2D form [16]:

$$\begin{cases} \mathbf{c}\_{\upsilon} \cdot \boldsymbol{\tau}\_{q} \cdot \frac{\partial^{2}T(\mathbf{x},\mathbf{y},t)}{\partial t^{2}} + \boldsymbol{c}\_{\upsilon} \cdot \frac{\partial T(\mathbf{x},\mathbf{y},t)}{\partial t} - \mathbf{k} \cdot \boldsymbol{\tau}\_{T} \cdot \frac{\partial \Delta T(\mathbf{x},\mathbf{y},t)}{\partial t} - \mathbf{k} \cdot \Delta T(\mathbf{x},\mathbf{y},t) = \\ \quad = q\_{V}(\mathbf{x},\mathbf{y},t) + \boldsymbol{\tau}\_{q} \cdot \frac{\partial q\_{V}(\mathbf{x},\mathbf{y},t)}{\partial t} \text{ for } \mathbf{x},\mathbf{y} \in \mathbf{R}, t \in \mathbf{R}\_{+} \cup \{0\} \end{cases} \tag{5}$$

As presented in Appendix A, the term τ*q* ∂*qv*(*t*) ∂*t* can be neglected with an error of about −1.91%, for the simulation time *t* ≥ 3<sup>τ</sup>T and |*err*| ≤ 0.58% for *t* ≥ 10<sup>τ</sup>*T*. The *qv<sup>t</sup>* + τ*q* can also be used instead of *qv*(*t*) + τ*q* ∂*qv*(*t*) ∂*t* for *qv*(*t*) = *c*1 · *<sup>H</sup>*(*t*), where *H*(*t*) is a Heaviside function, *c*1 is a constant where |*<sup>c</sup>*1| < +<sup>∞</sup>. It is worthwhile highlighting that the presented values were obtained considering the temperature and heat flux time lags for platinum. The Laplacian Δ*T* was approximated using the finite difference method (FDM) approach for a rectangular mesh with a constant distance between nodes in both dimensions. The considered approximation is presented below:

$$\begin{aligned} \Delta T(\mathbf{x}, y, t) &\approx \frac{T(\mathbf{x} + \Delta \mathbf{x}, y, t) + T(\mathbf{x}, y + \Delta \mathbf{x}, t) - 4 \cdot T(\mathbf{x}, y, t) + T(\mathbf{x} - \Delta \mathbf{x}, y, t) + T(\mathbf{x}, y - \Delta \mathbf{x}, t)}{\left(\Lambda \mathbf{x}\right)^2}, \\ \text{dla } &\mathbf{x}, y \in R, t \in R\_+ \cup \{0\} \end{aligned} \tag{6}$$

Then, the DPL Equation (5) becomes an ordinary differential equation (ODE) of a time variable. Such a constructed system of equations, including each investigated mesh node, was solved based on the backward differentiation formulas (BDF) approach [17–20], with initial conditions *<sup>T</sup>*(*<sup>x</sup>*,*y*,*<sup>t</sup>*) = 0 and ∂*T*(*<sup>x</sup>*,*y*,*<sup>t</sup>*)/∂*<sup>t</sup>* = 0 for *t* ≤ 0, which means zero initial conditions for Equation (16).
