Investigation of Heat Diffusion at Nanoscale Based on Thermal Analysis of Real Test Structure

Abstract

This paper presents an analysis related to thermal simulation of the test structure dedicated to heat-diffusion investigation at the nanoscale. The test structure consists of thin platinum resistors mounted on wafer made of silicon dioxide. A bottom part of the structure contains the silicon layer. Simulations were carried out based on the thermal simulator prepared by the authors. Simulation results were compared with real measurement outputs yielded for the mentioned test structure. The authors also propose the Grünwald–Letnikov fractional space-derivative Dual-Phase-Lag heat transfer model as a more accurate model than the classical Fourier–Kirchhoff (F–K) heat transfer model. The approximation schema of proposed model is also proposed. The accuracy and computational properties of both numerical algorithms are presented in detail.

1. Introduction

1.1. State of the Art

Nowadays, a heat transfer problem in the case of modern electronic structure is one of the most important research areas in the high-tech industry. The main reason for this fact is a significant downsizing of integrated circuits that is implemented in all modern electronic devices. Moreover, it also connects with a meaningful increase of an operation frequency of the mentioned devices. Both of these facts cause a rapid growth of a heat density generated inside such small structures. Consequently, the increased internal heat generation results in a huge increase of a temperature in critical parts of a device during its operation. It should also be highlighted that assurance of a proper cooling condition in the case of nanosized electronic structures is a non-trivial issue. Thus, all of these factors may influence an unstable operation and shorten the life cycle of an entire structure. Therefore, a thermal analysis seems to be one of the most important steps in the process of designing and developing modern electronic structures.

For a long time, heat conduction problems have been solved based on Fourier’s theory [1,2]. This approach uses the Fourier’s law and resulting Fourier–Kirchhoff (F–K) model for solids. They can be described by the following system of equations:

$$\{\begin{array}{cc}\begin{array}{c}q\left(x,t\right)=-k\cdot \nabla T(x,t)\\ c_v\frac{\partial T\left(x,t\right)}{\partial t}+\nabla \circ q\left(x,t\right)=q_V\left(x,t\right)\end{array}& x\in R^n,n\in N,t\in R_{+}\cup \left\{0\right\}\end{array}$$

(1)

where q is a heat flux density vector; k means a thermal conductivity of investigated material; T is a function regarding temperature rise above the ambient temperature; c_v means a volumetric heat capacity being a product of a specific heat of a material for a constant pressure (c_p) and its density (ρ); and q_V reflects the value of internally generated heat.

The classical F–K equation can also be directly derived from the classical thermodynamics for a positive-definite entropy-production rate, for which the thermodynamic state transition in a heat transfer process is extremely slow (quasi-stationary process). Therefore, the process time (t) should be longer than a system’s relaxation times.

Despite the numerous advantages of the Fourier–Kirchhoff approach, the research has shown that the application of the mentioned model may not reflect real phenomena [3,4,5]. One of its biggest disadvantages is an assumption considering the infinite speed of heat propagation. Apart from that, an instantaneous change of a temperature gradient or a heat flux should also be emphasized as a non-physical behavior postulated by the investigated theory. None of them, especially in the case of electronic structure with physical dimensions in nanometer-scale, has been empirically confirmed [6,7]. Thus, a thermal model including improvements of the F–K approach should be used instead of the classical theory.

The thermodynamic state is nonequilibrium in a fast transient process where time is comparable with the system relaxation times (e.g., an Umklapp phonon–phonon scattering process relaxation time in semiconductors [8,9]). The extended irreversible thermodynamic (EIT) can be applied to describe this kind of heat transfer system, assuming the second order the Taylor series expansion of entropy in EIT [10,11]. As a consequence, in mid-1990s, the new Dual-Phase-Lag (DPL) model was introduced by Tzou [12,13].

This model is an advanced version of an approach based on the Fourier–Kirchhoff theory and includes so-called time lags describing the time needed to change the temperature gradient, as well as the heat flux density. Each change is reflected by a separate lag value: a temperature time lag (τ_T) and a heat flux time lag (τ_q), respectively. The mathematical description of analyzed model can be presented in the form of the following system of equations:

$$\begin{array}{cc}\{\begin{array}{c}\nabla \circ q\left(x,t\right)=-c_v\frac{\partial T\left(x,t\right)}{\partial t}+q_V\left(x,t\right)\\ k{\tau }_T\frac{\partial \nabla T\left(x,t\right)}{\partial t}+{\tau }_q\frac{\partial q\left(x,t\right)}{\partial t}=-k\cdot \nabla T\left(x,t\right)-q\left(x,t\right)\end{array}& x\in R^n,n\in N,t\in R_{+}\cup \left\{0\right\}\end{array}$$

(2)

The other general approach can be derived by using the General Equation for Non-Equilibrium Reversible-Irreversible Coupling (GENERIC) equation proposed in 1997 by M. Grmela and H. C. Öttinger [14,15,16]. This approach can be useful to obtain a model for the Monte Carlo simulation of modern nanostructures [17] and also for deterministic approaches [18]. Apart from that, an approach known as Guyer-Krumhansl-type heat conduction [19,20] can be also used. Moreover, the research presented by Pop et al. [21], related to heat generation and transportation problems in nanosized transistors, is also worth considering.

Let us consider, more precisely, one of the most common heat transfer models for electronic structures—the DPL model (

Table 1. Selected Heat Transfer Model transformed into time delayed Partial Differential Equation.

Heat Transfer Model	Derived Time Delayed PDE
(a)
DPL equation [12,13] for τT > τq > 0, Maxwell–Cattaneo–Vernotte equation [26,27,28] for τq > 0, τT = 0 (more details in Vermeersch and De Mey as well as Kovács and Ván papers [29,30]), and F–K Equations (1)–(2) for τq = τT = 0: $\{\begin{array}{c}{c}_v\frac{\partial T\left(x,t\right)}{\partial t}=-\nabla \circ q\left(x,t\right)\\ q\left(x,t\right)+{\tau }_q\frac{\partial q\left(x,t\right)}{\partial t}=-k\cdot \nabla T\left(x,t\right)-k{\tau }_T\frac{\partial \nabla T\left(x,t\right)}{\partial t}\end{array}$	$c_v\frac{\partial T\left(x,t+{\tau }_q\right)}{\partial t}=k\cdot \Delta T\left(x,t+{\tau }_T\right)$ ⇓ $c_v\frac{\partial T\left(x,t+{\tau }_q-{\tau }_T\right)}{\partial t}=k\cdot \Delta T\left(x,t\right)$
(b)
Ballistic-conductive heat transfer model for τq > 0, τQ > 0, k12∙k12 ≤ 0 (more details [30]) $\{\begin{array}{c}\begin{array}{l}{c}_v\frac{\partial T\left(x,t\right)}{\partial t}=-\nabla \circ q\left(x,t\right)\\ {\tau }_q\frac{\partial q\left(x,t\right)}{\partial t}+q\left(x,t\right)=-k\cdot \nabla T\left(x,t\right)-k_{21}\cdot \nabla Q\left(x,t\right)\end{array}\\ {\tau }_Q\frac{\partial Q\left(x,t\right)}{\partial t}+Q\left(x,t\right)=k_{12}\cdot \nabla \circ q\left(x,t\right)\end{array}$	$\begin{array}{l}{c}_v\frac{\partial T\left(x,t+{\tau }_q\right)}{\partial t}=\\ =k\cdot \Delta T\left(x,t\right)-c_v{k}_{12}{k}_{21}\frac{\partial \Delta T\left(x,t-{\tau }_Q\right)}{\partial t}\end{array}$
(c)
Proposed DPL model [12,13], with fractional order of the temperature function space derivative based on Grünwald–Letnikov theory for τq > 0, τT > 0, 2 < αx < 2.5 (more details in [31]): $\begin{array}{l}{c}_v\frac{\partial T\left(x,t\right)}{\partial t}+{\tau }_q{c}_v\frac{{\partial }^{2}T\left(x,t\right)}{\partial t^2}=\\ =k\cdot {}_{GL}{\Delta }^{{\alpha }_x}T\left(x,t\right)+k{\tau }_T\frac{\partial {}_{GL}\Delta {}^{{\alpha }_x}T\left(x,t\right)}{\partial t}\end{array}$	$c_v\frac{\partial T\left(x,t+{\tau }_q\right)}{\partial t}=k\cdot {}_{GL}\Delta {}^{{\alpha }_x}T\left(x,t+{\tau }_T\right)$ ⇓ $c_v\frac{\partial T\left(x,t+{\tau }_q-{\tau }_T\right)}{\partial t}=k\cdot {}_{GL}\Delta {}^{{\alpha }_x}T\left(x,t\right)$

a—left column), the ballistic-conductive heat transfer model (Table 1b—left column), and the model proposed in this paper, the DPL model with fractional order of the temperature function space derivative based on the Grünwald–Letnikov theory (Table 1c—left column). Suppose also an isotropic medium parameter properties and idempotent equations parameters k > 0, c_v > 0, k₂₁, k₁₂, τ_T, τ_q, τ_Q, etc. To simplify the analysis, the Taylor-series expansions mapping into a function with time delay is considered, as presented in Equation (3). The terms of higher orders are neglected.

$$f\left(x,t\right)+\tau \frac{\partial f\left(x,t+\tau \right)}{\partial t}\to f\left(x,t+\tau \right)$$

(3)

Then, time delayed PDEs are obtained and presented in the right column of Table 1. In the first case (Table 1a—right column), the state is appointed by using a gradient of a heat flux (−k∙ΔT) delayed by τ_T-τ_q. It should be emphasized that results produced by the DPL model are consistent with many experiments and measurements at nanoscale (for more, see [9]).

In the case of the ballistic-conductive heat transfer model (Table 1b), the rate of value temperature T(x, t + τ_q) change at t+ τ_q and x is dependent on a gradient of heat flux (−k∙ΔT) at the current time (t), increased by a difference of growth rate of the averaged value of temperature (T) over the infinitesimal neighborhood of point x and the value of temperature (T) in this point, both from the past time (t − τ_Q):

$$DT\left(x,t+{\tau }_q\right)\propto k\cdot \Delta T\left(x,t\right)\cdot Dt+\left|c_v{k}_{12}{k}_{21}\right|\cdot D\Delta T\left(x,t-{\tau }_Q\right)$$

(4)

where D(∙) is the difference operator (“change in”) corresponding to changes for Dt→0 and also for k₁₂∙k₂₁ ≤ 0. Therefore, the temperature T(x, t + τ_q) dynamic is additionally intensified by the dynamics of the rate of heat flux gradient (∂ΔT(x, t − τ_Q)/∂t) in the past time (t − τ_Q).

The approach proposed in this work is based on the temperature T(x, t + τ_q) dynamic control, using the fractional order of Laplace operator (_GLΔ^αx) in the DPL model (Table 1c—right column). The integral or differential behavior of this operator and time–space discretization scheme is obtained by the value of parameter α_x. The application of the fractional derivative can be interpreted as an inhomogeneous space for the heat transfer in relation to the local energy distribution (interpreted as temperature at nanoscale) in infinitesimal neighborhood of considered point. Therefore, the infinitesimal distance in Laplacian is modulated in relation to the temperature distribution around the considered point at nanoscale. The physical interpretation of the fractional derivatives is also presented in [22,23,24,25].

Hence, the proposed model (and proposed time–space scheme) is a link between experimentally confirmed DPL mesoscopic model with the ballistic heat transport model with dynamic temperature changes’ intensification useful for quasi 1-D nanostructures and for radiative heat transport without phonon collisions (e.g., metal nanowires and ultra-thin metal–oxide films).

One of the biggest advantages of the DPL approach is a universality of its application. It can be used in the case of parabolic partial differential equations, as well as for hyperbolic ones. It means that the Dual-Phase-Lag model can successfully replace the classical model based on parabolic Fourier–Kirchhoff differential equation. Of course, the DPL approach also has some disadvantages. For example, it is characterized by a bigger complexity than F–K model, what results in longer time needed for computational simulation results, especially in the case of complex structures. However, this issue will not be investigated in this paper.

The crucial achievement described in this paper is the application of the Grünwald–Letnikov fractional space-derivative in Dual-Phase-Lag heat transfer model to more accurate modeling of delays and modulation of Laplacian in relation to the temperature distribution around the considered point. This approach seems to be beneficial to represent real heat transfer behavior at the nanoscale. The proposed solution is more accurate than the F–K heat transfer model. The computation time reduction is also obtained in comparison to the DPL model. Mentioned features will be demonstrated by using a special test nanostructure developed in microelectromechanical system (MEMS) technology.

The structure of the paper is as follows. Primarily, a description of investigated real test structure is prepared. Then, analyses related to mathematical modeling of the temperature distribution inside the test structure are included. The structure’s discretization method is considered, as well. Key components of authors’ approach, i.e., DPL model approximation scheme and heat transfer enhancement (including Grünwald–Letnikov fractional derivative and its application), are described in Section 2.3 and Section 2.4. Next, simulation results are demonstrated and compared to real measurements described in [32,33]. Finally, results are discussed, and the research is summarized.

1.2. MEMS Test Structure Description

The MEMS test nanostructures were manufactured in the Polish Institute of Electron Technology. The structure consists of two parallel platinum resistors, which have lengths equal to 10 µm. The cross-sectional area of each resistor is a 100 nm square. Moreover, the distance between both resistors is also equal to 100 nm. They are placed on a thin layer of silicon dioxide, with a thickness of 100 nm. All layers are stacked on a 0.5 mm thick silicon wafer. A more detailed structure description can be found in Janicki et al. [33].

The test structure was bonded to a metal-core PCB characterized by high thermal conductivity. Moreover, the structure was connected to a biasing circuit, using high-frequency coaxial cables that were mounted to a tiny Hirose U.FL connector attached to the silicon die. Real photos of the analyzed structure, as well as the control circuit, are included in Janicki et al. [32,33].

During the measurement process, one of platinum resistors played a role of a heater, while the second one served as a temperature sensor. A detailed description of the measurement process of the test structure and obtained results is presented in [32,33].

2. Mathematical Description of Proposed Methodology

2.1. General Description

Thermal simulation was performed for the two-dimensional cross-sectional area of the investigated structure in the middle of the resistors’ length. To obtain the temperature simulation results, the DPL model, Equation (2), was used. In order to make the analysis more effective, the system of Equation (2) was transformed to an equivalent form, presented below, for two-dimensional space [34]:

$$\begin{array}{c}{c}_v\cdot {\tau }_q\cdot \frac{{\partial }^{2}T\left(x,y,t\right)}{\partial t^2}+c_v\cdot \frac{\partial T\left(x,y,t\right)}{\partial t}-k\cdot {\tau }_T\cdot \frac{\partial \Delta T\left(x,y,t\right)}{\partial t}-k\cdot \Delta T\left(x,y,t\right)=q_V\left(x,y,t\right)\\ x,y\in R,t\in R_{+}\cup \left\{0\right\}\end{array}$$

(5)

The Laplacian of a temperature function ΔT was approximated by using the Finite Difference Method (FDM) according to the following formulas:

$$\Delta T\left(x,y,t\right)=\frac{{\partial }^{2}T\left(x,y,t\right)}{\partial x^2}+\frac{{\partial }^{2}T\left(x,y,t\right)}{\partial y^2},x,y\in R,t\in R_{+}\cup \left\{0\right\}$$

(6)

$$\frac{{\partial }^{2}T\left(x,y,t\right)}{\partial x^2}\approx \frac{T\left(x+\Delta x,y,t\right)-2\cdot T\left(x,y,t\right)+T\left(x-\Delta x,y,t\right)}{{\left(\Delta x\right)}^2},x,y\in R,t\in R_{+}\cup \left\{0\right\}$$

(7)

$$\frac{{\partial }^{2}T\left(x,y,t\right)}{\partial y^2}\approx \frac{T\left(x,y+\Delta y,t\right)-2\cdot T\left(x,y,t\right)+T\left(x,y-\Delta y,t\right)}{{\left(\Delta y\right)}^2},x,y\in R,t\in R_{+}\cup \left\{0\right\}$$

(8)

Thus, considering the same difference between nodes in both axes, i.e., Δx = Δy, Laplacian ΔT can be approximated in the following way:

$$\begin{array}{l}\Delta T\left(x,y,t\right)\approx \frac{T\left(x+\Delta x,y,t\right)+T\left(x,y+\Delta x,t\right)-4\cdot T\left(x,y,t\right)+T\left(x-\Delta x,y,t\right)+T\left(x,y-\Delta x,t\right)}{{\left(\Delta x\right)}^2},\\ \mathrm{dla}x,y\in R,t\in R_{+}\cup \left\{0\right\}\end{array}$$

(9)

On the basis of this methodology, the authors’ method for structure discretization and FDM matrices generation was proposed. Moreover, taking into consideration the proposed approximation, the DPL equation, in the form of Equation (5), has become an ordinary differential equation of a time variable. Finally, the prepared matrix system of equations was solved for different points of time, using a class of Backward Differentiation Formulas (BDF) [35,36,37,38].

In the following subsections, the structure discretization, as well as the proposed discretization scheme for DPL model describing the temperature distribution in the cross-sectional area of investigated test structure, is characterized in detail.

2.2. Structure Cross-Sectional Area Discretization

Primarily, the investigated cross-sectional area of the structure, as presented in Figure 1, was discretized by using two-dimensional discretization mesh characterized by the following formulas [34]:

$$q_k\left(t\right)=q\left(x,y,t\right),x=i\cdot \Delta x,y=j\cdot \Delta y$$

(10)

$$T_k\left(t\right)=T\left(x,y,t\right),x=i\cdot \Delta x,y=j\cdot \Delta y$$

(11)

i ∊ {1, 2,…, n_x}, j ∊ {1, 2,…, n_y}, k ∊ {1, 2,…, n_x ∙n_y}

where n_x and n_y describe a number of discretization nodes in the x-axis and y-axis, respectively. On the other hand, the product n_x∙n_y reflects the entire number of nodes used to discretize the structure’s cross-section.

Nodes are numbered from the left to the right side, along the x-axis. After reaching the last point in a single row, the next part of the structure, being the nearest row from the top of the current row, was taken into consideration and numbered in the same way. Thus, node no. 1 was placed in the left bottom corner of the structure, while the node with the highest possible number, equal to n_x∙n_y, was located in the top right corner. The graphical representation of used discretization mesh for analyzed cross-section of the test structure, as well as the way mesh nodes were numbered, is demonstrated in Figure 2. It is also worth highlighted that the distance between nodes in both dimensions was set to 10 nm.

2.3. Dual-Phase-Lag Approximation Scheme for Test Structure

In order to obtain the solution of temperature distribution problem inside the investigated cross-sectional area of the structure, the Dirichlet boundary conditions were considered at the bottom, while left, right and the top parts have been modeled using Neumann boundary conditions. Described situation can be characterized by the following equations:

$$\begin{array}{c}{T}_k\left(t\right)=0,t\in R_{+}\cup \left\{0\right\},\\ k\in \left\{1,2,3,\dots ,n_x\right\}\end{array}$$

(12)

$$\begin{array}{c}{q}_k\left(t\right)=0,t\in R_{+}\cup \left\{0\right\},\\ 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_x,2\cdot n_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{array}$$

(13)

Considering imposed boundary conditions and the FDM assumptions, the DPL model can be explained by the following system:

$$\{\begin{array}{cc}\begin{array}{c}M\cdot \ddot{T}\left(t\right)+D\cdot \dot{T}\left(t\right)+K\cdot T\left(t\right)=b\cdot u\left(t\right)\\ y\left(t\right)=c^T\cdot T\left(t\right)\end{array}& t\in R_{+}\cup \left\{0\right\}\end{array}$$

(14)

where index T means a transposition, while $\mathit{T}$, $\dot{\mathit{T}}$, and $\ddot{\mathit{T}}$ are n_x∙n_y-element vectors reflecting the temperature function and its first and the second time derivatives, respectively. Moreover, taking into account the way of nodes numbering, matrices I_diag, M_FDM, M, D, K, and c^T; vectors c_v, k, τ_q, τ_T, b, and y; and the function u(t) may be reflected as indicated below:

$$\begin{array}{l}{I}_{diag},M_{FDM},M,D,K,c^T\in R^{n_x\cdot n_y\times n_x\cdot n_y},\\ c_v,k,{\mathsf{\tau }}_q,{\mathsf{\tau }}_{T}b,y\in R^{n_x\cdot n_y\times 1},\\ u\in R\end{array}$$

$$M_{FDM}=\left[\begin{array}{ccccccccccc}-3& 1& 0& \ddots & 1& 0& 0& \ddots & 0& 0& 0\\ 1& -4& 1& \ddots & 0& 1& 0& \ddots & 0& 0& 0\\ 0& 1& -4& \ddots & 0& 0& 1& \ddots & 0& 0& 0\\ \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots \\ 1& 0& 0& \ddots & -4& 1& 0& \ddots & 1& 0& 0\\ 0& 1& 0& \ddots & 1& -4& 1& \ddots & 0& 1& 0\\ 0& 0& 1& \ddots & 0& 1& -4& \ddots & 0& 0& 1\\ \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots \\ 0& 0& 0& \ddots & 1& 0& 0& \ddots & -3& 1& 0\\ 0& 0& 0& \ddots & 0& 1& 0& \ddots & 1& -3& 1\\ 0& 0& 0& \ddots & 0& 0& 1& \ddots & 0& 1& -2\end{array}\right]$$

(15)

$$M=\mathrm{diag}\left({\mathsf{\tau }}_q\right)\circ \mathrm{diag}\left(c_v\right)\circ I_{diag}$$

(16)

$$D=\mathrm{diag}\left(c_v\right)\circ I_{diag}-\frac{1}{{\left(\Delta x\right)}^2}\cdot \mathrm{repmat}\left({\mathsf{\tau }}_T\right)\circ \mathrm{repmat}\left(k\right)\circ M_{FDM}$$

(17)

$$K=-\frac{1}{{\left(\Delta x\right)}^2}\mathrm{repmat}\left(k\right)\circ M_{FDM}$$

(18)

$$b=a\cdot {\left[\begin{array}{ccccccccc}0& \cdots & 0& 1& \cdots & 1& 0& \cdots & 0\end{array}\right]}^T,a\in R_{+}$$

(19)

$$c=I_{diag}$$

(20)

where I_diag is an identity matrix, operator ◦ reflects matrices multiplication resulting in a Hadamard product, and matrix function diag (∙) creates a diagonal matrix from a vector, while Matlab repmat function replicates a given vector and composes a matrix of required dimensions. It is also worth highlighting that non-zero values in vector b are observed only in the case of mesh nodes characterizing a heating area.

The system in Equation (14) of differential equations of the second order was equivalently transformed into the following system of the linear first-order equations [39]:

$$\{\begin{array}{cc}\begin{array}{c}E\cdot \dot{\overline{T}}\left(t\right)=A\cdot \overline{T}\left(t\right)+B\cdot u\left(t\right)\\ y\left(t\right)=C^T\cdot \overline{T}\left(t\right)\end{array}& t\in R_{+}\cup \left\{0\right\}\end{array}$$

(21)

where $\overline{T}$ and $\dot{\overline{T}}$ are 2∙n_x∙n_y-element vectors. The first part of vector $\overline{T}$ consists of elements of vector T, while its second part includes elements of $\dot{\mathit{T}}$. On the other hand, the first and the second parts of the vector $\dot{\overline{T}}$ consists of elements of vectors $\dot{\mathit{T}}$ and $\ddot{\mathit{T}}$, respectively. Moreover, matrices E, A, B, and C^T can be represented in the following way [39,40]:

$$E=\left[\begin{array}{cc}{I}_{diag}& \Theta \\ \Theta & M\end{array}\right],A=\left[\begin{array}{cc}\Theta & I_{diag}\\ -K& -D\end{array}\right],B=\left[\begin{array}{c}{\Theta }_1\\ b\end{array}\right],C=\left[\begin{array}{cc}c& \Theta \\ \Theta & c\end{array}\right]$$

(22)

where Θ and Θ₁ are null matrices. Moreover, we get the following:

$E,A,C^T\in R^{2\cdot n_x\cdot n_y\times 2\cdot n_x\cdot n_y},B\in R^{2\cdot n_x\cdot n_y\times 1},I,\Theta \in R^{n_x\cdot n_y\times n_x\cdot n_y},{\Theta }_1\in R^{n_x\cdot n_y\times 1}$

Thus, the full authors’ approximation scheme for the DPL model in two-dimensional Euclidean space is described by the system in Equation (21), including explanations formulated in Equations (15)–(22) and boundary conditions in Equations (12) and (13). The final solution of temperature-distribution changes over time in the analyzed cross-section of the test structure, based on the prepared approximation scheme, is determined by using the BDF method of a variable order between 1 and 5.

2.4. Heat Transfer Enhancement

The considered test structure was analyzed, including the surrounding air environment. Furthermore, platinum resistors’ separation distance d = 100 nm in the benchmark structure is comparable to the surrounded air mean free path length Λ ≈ 65 nm. Therefore, the heat flow can be approximated by the following equation [41]:

$$q_{air}=\frac{k_{air}}{\langle d+a\Lambda \rangle }\left(T_{Pt,A}-T_{Pt,B}\right)$$

(23)

where a is used to describe the interaction between the gas molecules and the solid walls [42], typically a = 1, and <∙> stands for ensemble average. The final equivalent air conductance between the mentioned resistors was estimated by using the following simplified formula (compare with value in

Table 2. Considered material parameters’ values [46].

Layer	Material	$\mathit{k}\left[\frac{\mathbf{W}}{\mathbf{m}\cdot \mathbf{K}}\right]$	$\mathit{\rho }\left[\frac{\mathbf{kg}}{{\mathbf{m}}^3}\right]$	${\mathit{c}}_{\mathit{p}}\left[\frac{\mathbf{J}}{\mathbf{kg}\cdot \mathbf{K}}\right]$
1 (wafer)	Silicon (Si)	148	2330	712
2 (oxide)	Silicon dioxide (SiO2)	1.38	2220	745
3 (heater)	Platinum (Pt)	71.6	21,450	133
4 (thermometer)	Platinum (Pt)	71.6	21,450	133
5 (ambient)	Air	0.0263 *	1.1614	1.007

* See Section 2.4.

):

$$k_{air,eqv,PtA-PtB}=k_{air}\frac{d}{\langle d+a\Lambda \rangle }\approx 0.961k_{air}\approx 0.961\cdot 0.0263\mathrm{W}/\left(\mathrm{mK}\right)$$

(24)

The photon tunneling and the radiative heat transfer between platinum resistor parallel surfaces were also analyzed. The photon tunneling was neglected due to the small amount of heat flux transport ${\left(S_z\right)}_{\max }^{evamescant}=2.2\cdot {10}^{-19}{\mathrm{W}/\mathrm{m}}^2$ (calculated for vacuum environment) in comparison to the main heat flux Pt-SiO₂ q_{Pt_SiO2},_cond≈19.3 MW/m² and the heat conduction through q_{Pt_Pt,cond} ≈ 0.917 q_Pt-SiO2,cond >> ${\left(S_z\right)}_{\max }^{evamescant}$ in the air for the steady state [43,44]:

$${\left(S_z\right)}_{\max }^{evamescant}\approx \frac{k_B^2\cdot \left(T_{PtA}+273.15\mathrm{K}\right)}{24\hslash b^2}\left[{\mathrm{W}/\mathrm{m}}^2\right]$$

(25)

where b is an inter-atomic distance b ≈ 39.12 nm for Pt, k_B means the Boltzmann constant, and $\hslash$ is the reduced Planck constant. Moreover, the radiative heat transfer (q_{SB rad}) was also neglected, due to the insignificant emitted energy from warmer (T_PtA) to colder platinum resistor (T_PtB) parallel surfaces q_Pt-SiO2, _cond / q_{SB rad} ≈ 8.3∙10²⁵:

$$\begin{array}{l}{q}_{SBrad}\approx 5.6693\frac{\mathrm{W}}{{\mathrm{m}}^2\left(\frac{\mathrm{K}}{100}\right)}\cdot \left(\sqrt{2}-1\right)\cdot \left\{{\left(\frac{T_{PtA}+273.15\mathrm{K}}{100}\right)}^4-{\left(\frac{T_{PtB}+273.15\mathrm{K}}{100}\right)}^4\right\}\\ \mathrm{and}72\le q_{SBrad}\le 146{\mathrm{W}/\mathrm{m}}^2\mathrm{for}{T}_{PtA}\approx 59\mathrm{K}\mathrm{and}4.7\mathrm{K}\le T_{PtB}\le 31\mathrm{K}\end{array}$$

(26)

Moreover, considering investigation presented above, some additional changes in the proposed model are needed. Thus, for the air area between resistors’ surfaces and for contact areas between platinum and silicon dioxide, a fractional order of the temperature-rise function space derivative has been employed, according to the theory described in [31,45]. This theory, based on the Grünwald–Letnikov definition of the fractional derivative, allows us to establish the following formula reflecting the temperature-rise function’s improvement for the central difference of the FDM scheme [31,45]:

$$\begin{array}{l}{}_{GL}{D}_{0,s}^{{\alpha }_x}T\left(s\right)=\frac{1}{{\left(\Delta s\right)}^{\alpha }}\cdot {\displaystyle \sum _{k=0}^{\mathrm{round}\left({\alpha }_x,0\right)}{\left(-1\right)}^k\frac{\Gamma \left({\alpha }_x+1\right)}{\Gamma \left(k+1\right)\cdot \Gamma \left({\alpha }_x-k+1\right)}T\left(s-k\cdot \Delta s+\frac{{\alpha }_x\cdot \Delta s}{2}\right)},\\ \mathrm{for}{\alpha }_x\in R_{+},\Delta s\to 0\end{array}$$

(27)

where Δs is the mesh nodes distance, α is investigated fractional order, Γ is the special gamma function, and round(m,n) is the function rounding the value m to n digits. Considering fractional order of derivative, we needed approximating investigated function values in points between nodes of a discretization mesh. The mentioned approximation depends on function values in neighboring mesh points. Thus, the right-hand-side part of Equation (27), being an investigated approximation, can be reflected by using the following equations [45]:

$$\begin{array}{l}\frac{1}{{\left(\Delta s\right)}^{{\alpha }_x}}\cdot {\displaystyle \sum _{k=0}^2{\left(-1\right)}^k\frac{\Gamma \left({\alpha }_x+1\right)}{\Gamma \left(k+1\right)\cdot \Gamma \left({\alpha }_x-k+1\right)}T\left(x-k\cdot \Delta x+\frac{{\alpha }_x\cdot \Delta x}{2},t\right)}\cdot 1=\\ =\frac{\left(\frac{{\alpha }_x}{2}-1\right)\cdot T\left(x+2\cdot \Delta x,y,t\right)+\left(2-\frac{{\alpha }_x}{2}-{\alpha }_x\cdot \left(\frac{{\alpha }_x}{2}-1\right)\right)\cdot T\left(x+\Delta x,y,t\right)}{{\left(\Delta x\right)}^{{\alpha }_x}}+\\ +\frac{\left(\frac{{\alpha }_x\cdot \left({\alpha }_x-1\right)}{2}\cdot \left(\frac{{\alpha }_x}{2}-1\right)-{\alpha }_x\cdot \left(2-\frac{{\alpha }_x}{2}\right)\right)\cdot T\left(x,y,t\right)+\left(\left(2-\frac{{\alpha }_x}{2}\right)\cdot \frac{{\alpha }_x\cdot \left({\alpha }_x-1\right)}{2}\right)\cdot T\left(x-\Delta x,y,t\right)}{{\left(\Delta x\right)}^{{\alpha }_x}}\\ +\frac{\left(\frac{{\alpha }_x}{2}-1\right)\cdot T\left(x,y+2\cdot \Delta x,t\right)+\left(2-\frac{{\alpha }_x}{2}-{\alpha }_x\cdot \left(\frac{{\alpha }_x}{2}-1\right)\right)\cdot T\left(x,y+\Delta x,t\right)}{{\left(\Delta x\right)}^{{\alpha }_x}}+\\ +\frac{\left(\frac{{\alpha }_x\cdot \left({\alpha }_x-1\right)}{2}\cdot \left(\frac{{\alpha }_x}{2}-1\right)-{\alpha }_x\cdot \left(2-\frac{{\alpha }_x}{2}\right)\right)\cdot T\left(x,y,t\right)+\left(\left(2-\frac{{\alpha }_x}{2}\right)\cdot \frac{{\alpha }_x\cdot \left({\alpha }_x-1\right)}{2}\right)\cdot T\left(x,y-\Delta x,t\right)}{{\left(\Delta x\right)}^{{\alpha }_x}}\\ \mathrm{for}{\alpha }_x\in \left(2,2.5\right),\Delta x\to 0\end{array}$$

(28)

The formula above replaces the classic approximation of the Laplacian described by Equation (9).

3. Thermal Simulation and Results Analysis

3.1. Material Characterization and Initial Simulation Results

A thermal simulation of the test structure was prepared by using MathWorks® Matlab environment and proposed author’s approximation scheme for the DPL model. The used computational node includes 4-core, 8-threads Intel^® Core™ i7 2.6 GHz (3.6 GHz in Turbo mode) CPU, 32 GB DDR4 memory supported by 265 GB of swap file. In order to obtain simulation results, parameters’ values presented in Table 2 were taken into consideration.

The most problematic issue is related to establishing parameters of the air layer, being an ambient of investigated test structure. In particular, the air between two platinum resistors is crucial in presented investigation. Moreover, the contact layer between platinum resistors and the oxide layer also needs special consideration. In order to emphasize the problem, a sample analysis of temperature rises over the time were prepared for different values of the DPL model parameters, as well as for different values of a thermal conductivity for the mentioned part of the air layer.

The first part of the simulation process does not include the investigation demonstrated in Section 2.4. Two sets of reference DPL model parameters were considered during the simulation process:

• τ_T = 60 ps, τ_q = 3 ps (all layers)
• τ_T = 2.6 ps, τ_q = 0.0916 ps (platinum resistors) and τ_T = 60 ps, τ_q = 3 ps (all remaining layers).

The results, as demonstrated in Figure 3 and Figure 4, were additionally normalized in order to make their analysis easier. The normalization was prepared in the following way:

$$T_k^{norm}\left(t\right)=\frac{T_k\left(t\right)}{\underset{t,k}{\max }\left\{T_k\left(t\right)\right\}}k\in \left\{1,2,\dots ,n_x\cdot n_y\right\},t\in R_{+}\cup \left\{0\right\}$$

(29)

As it can be seen, the average temperature rise inside the heater is less than 95% of the maximal recorded temperature rise, while the temperature rise in the platinum sensor is characterized by a slightly more than 61% of the highest observed temperature rise. Differences in temperature rise values yielded for analyzed sets of DPL model parameters are observed between 1 ps and 1 ns.

Time shifts between observed lines in Figure 3 and Figure 4 were plotted in Figure 5. In order to show differences between results, excluding and including air conductivity investigation presented in Section 2.4, a time shift analysis over the time, demonstrated in Figure 6, was carried out. In this figure, “new k_air” means the thermal conductivity of the air layer between platinum resistors calculated based on the analysis shown in Section 2.4. As a comparison, the simulation results obtained for both DPL time lags equal to zero for the air layer are also included.

Taking into consideration the sensitivity of yielded results to even small changes of the air layer material parameters, in the second part of the simulation, values, calculated considering the investigation described in Section 2.4, were employed only.

3.2. Final Simulation Results and Comparison to Real Measurements of the Test Structure

The analysis in the previous subsections shows that there is a need for calculation of the proper material parameters for the air layer, especially between platinum resistors of the test structure. Moreover, results of further research suggest using the fractional order of the space derivative of the temperature rise function for the investigated region, as well as this one between platinum resistors and the silicon dioxide, according to the theory described in Section 2.4. Taking into consideration these facts, the simulation of the temperature distribution in the cross-sectional area of the real test structure was carried out.

It was assumed that differences between mesh nodes are equal to 10 nm in both axes. Furthermore, values of DPL model time-lag parameters were calculated according to the theory presented in Zubert et al. [8]. Thus, for Platinum resistors, heat-flux time lag was set at approximately 550 ps, while the considered value of the temperature time lag was equal to 15.6 ns. In the case of other materials, investigated parameters were equal to 18 and 480 ns, respectively. Simulation results and their comparison to the real measured data (collected and described in [32,33]) are presented in Figure 7 and Figure 8, for the transient and steady state analyses, respectively. Moreover, in Figure 7, simulation results received by using the classical F–K model are also included, for comparison purposes.

In Figure 7, the results for the heated platinum resistor were marked by dashed and dotted lines, while those ones obtained for the temperature sensor were plotted with dashed lines. Moreover, outputs for the heater and thermometer were marked by the red and blue curves, respectively. On the other hand, measurement data (collected and described in Janicki et al. [32,33]) were plotted, using green lines, while results obtained by using the F–K model were marked by the black color.

The maximal temperature rise above the ambient temperature was observed at the top surface of the heat source, i.e., in the first platinum resistor (layer 3). Its value is nearly 60 K. On the other hand, the temperature rise recorded at the surface of the temperature sensor, i.e., the second platinum resistor (layer 4), is about 6 K, which states approximately 10% of the temperature rise observed in the heat source.

Both of the investigated temperature rises coincide almost exactly with measurements of the real structure (green curves in Figure 7). Parameter α_x allows for changing the result curves’ slopes, while DPL model parameters τ_q and τ_T cause the curves to shift over the time. The mentioned change is proportional to the parameters’ values.

In order to check the quality of the obtained results for used the discretization mesh and mesh nodes distances (10 nm), which were described in Section 2.2, the simulation curves’ fitting to the measured one was considered. Each curve was plotted based on 703 points logarithmically distributed over the time. Then, such metrics as coefficient of determination (R2), root mean squared error (RMSE), sum of squared estimate error (SSE), and mean squared error (MSE) were calculated. Moreover, to assess the goodness of recognition as a volatility over the time, a Pearson correlation coefficient (corr) was also considered. Determined metrics values are presented in

Table 3. Evaluation of goodness of fitting of simulation results to real data.

Temperature Distribution Simulation	MSE	RMSE	SSE	R2	Corr
Heater	8.5837∙10−4	0.0293	0.6034	0.9572	0.9973
Thermometer	9.5540∙10−5	0.0098	0.0472	0.9554	0.9748

.

As it can be seen, both the heater and thermometer curves are characterized by relatively small values of MSE, RMSE, and SSE values. Moreover, the coefficient of determination and correlation coefficient suggest a proper recognition of a shape of the measured curve by a simulated one. Generally, it can be stated that calculated metrics (MSE, RMSE, SSE, coefficient of determination, and correlation coefficient) for the simulation fitting to the real data confirm highly accurate simulation results. This situation clearly shows that the proposed approach based on the theory described in Section 2 allows for the production of outputs reflecting the real thermal phenomena observed at the nanoscale. Moreover, as it was shown in Figure 7, the classical F–K model should not be used in the case of electronic nanosized structures.

4. Conclusions

This paper includes the investigation of the heat transfer problems at the nanoscale. A new approach to the heat transfer modeling in modern nanosized structures was considered. It combines the DPL model and the Grünwald–Letnikov fractional derivative. A combination of these mathematical tools allows for the preparation of a complex approach using the Finite Difference Method to temperature distribution determination at the nanoscale, with a high level of accuracy, as is confirmed by the measurement of a real test structure.

An important novelty described in this paper is the use of a DPL model with a fractional order space derivative of a temperature function based on the Grünwald–Letnikov derivative operator. This operator, as well as the proposed time–space discretization schema, is a bridge between experimentally confirmed DPL mesoscopic model with the ballistic heat transport model with dynamic temperature changes’ intensification useful for quasi 1-D nanostructures and for radiative heat transport without phonon collisions. This solution allows for the consideration of such physical behaviors like time needed for a heat flux or temperature gradient changes. Thus, a modeling of a heat diffusion can be investigated by making realistic assumptions, which was not possible in the case of the F–K model use and real relaxation thermal properties of material at mesoscopic scale. The research has shown that the proposed GL DPL model is more realistic than the commonly used Fourier–Kirchhoff model.

The manuscript describes also proposed an approximation scheme of a modern DPL heat transfer model based on the Finite Difference Method approach prepared for the two-dimensional cross-section of the real test nanometric electronic structure manufactured at the Institute of Electron Technology in Warsaw. The investigation has shown that there is a possibility to effectively implement a prepared algorithm that allows for the determination of a temperature distribution inside real nanoscale electronic structures, based on proposed an approximation scheme.

Thermal simulation has provided results which coincide almost exactly with the real measurements. It means that prepared methodology is highly accurate and allows modeling of the heat transfer problems by using a modern approach based on the use of the Dual-Phase-Lag model. The considered thermal model is an appropriate methodology for heat-diffusion modeling, especially at the nanoscale.

In the future, the reduction of the DPL model order reduction methodology will be considered in order to save simulation time, decrease a computational power requirements [47], and make the simulation process more efficient.