Understanding the Thermal Time Constants of GaN HEMTs through Model Order Reduction Technique

: This paper described a comparison between a numerical Finite Element Analysis (FEA) and an analytical approach in order to extract the thermal time constants and the thermal resistances of simple but realistic structures. Understanding the complex contribution of multidimensional thermal spreading, the effect of multiple layers, and the correlation with the heat source length is mandatory due to the severe mismatch of thermal expansion in different epitaxial layers and high operating temperatures. This is especially true on GaN HEMT (High Electron Mobility Transistor) with the continuous decrease of the gate length and the increase of the power density. Moreover, in this paper, we extracted the time constants with a Model Order Reduction (MOR) technique based on the Ritz vector approach with inputs coming from the FE matrices. It was found that the time constants obtained by an analytical solution and a model order extraction from FEA were exactly the same. This result validated the idea that our MOR technique provides the real time constants and resistances for our device structures and in this case uniﬁed the analytical world with the numerical one.


Introduction
In recent years, the gallium nitride (GaN) technology has gained increasing interest in high-frequency power amplifier and high-voltage power conversion applications [1]. The high-power density present in GaN HEMTs technology leads to a high temperature in the transistor channel. This self-heating phenomenon degrades the performance and the reliability of the device. Knowing the channel temperature under operating conditions, the location of the hotspot and transient temperature response have become mandatory in the RF electronic industry.
Thermal issues of GaN HEMTs have been studied from various aspects. A large number of analytical models for thermal resistance have been developed. However, they are in many cases restricted to simple geometries, uniform boundary conditions, and linear thermal conductivity. Albers [2] and Yovanovich [3] used Fourier series-based formulation to calculate the surface temperature solution for a 3D rectangular structure with an arbitrary number of layers. Geer et al. [4] developed a semi-analytical approach that considers general boundary and interface conditions as well as internal heat sources. Masana solved the surface temperature for a 2D structure [5], and in [6], Darwish proposed a closed-form analytical model based on a prolate spheroidal and ellipsoid coordinate system to determine the thermal resistance of multi-finger AlGaN-GaN devices. These models have focused on the issue of steady-state temperature and thermal resistance; however, few works have been published on the understanding of the thermal time constants. Bagnall [7,8] solved the heat equation for a single layer with the heat flux evenly distributed at the top and a constant temperature at the bottom. He also provided an Bagnall [7,8] solved the heat equation for a single layer with the heat flux evenly distributed at the top and a constant temperature at the bottom. He also provided an analytical formula for thermal time constants. In all these approaches, the Fourier's Law is only considered. It is known that the Fourier's approach may be limited for nanoscale devices [9,10], but our focus of approach was to address and improve thermal models for electronics industry beyond a simple RC cell for real devices.
In this paper, a comparison between the thermal time constants and thermal resistances obtained both using the analytical model [7] and numerical model (ANSYS software) coupled with MOR [11] was carried out. This work can be considered as an extended study of [7]. However, in our article, the time constants were extracted through model order reduction from a numerical Partial Differential Equation. The excellent part is that this model order reduction technique was based on a numerical approach, thereby yielding the real physical time constants. The way of extracting the thermal time constants was found to be quite different in [7], because the authors used a fitting procedure associating with a modeling equation and a minimization of an error function.
The motive of this work was essentially to prove that the real physical time constants can be determined by using our MOR technique.
The flow presented in Figure 1 represents the process used to obtain the paired values (Ri, τ ). Part 2 presents the process of extracting the thermal time constants as well as the thermal resistance with our MOR technique.
In part 3, a one-dimensional analysis was performed for single-layer and two-layer structures. A comparison between the analytical approach and the numerical one is detailed.
In part 4, we consider a two-dimensional approach that takes into account the length of the heat source.
In part 5, some simulation results are compared with the experimental data in order to validate the approach proposed in this article.

Extraction of Real Thermal Time Constants with Ritz Model Order Reduction Technique
Nonlinear Model order Reduction is a very complex task. The assumption for the following development was that thermal properties are constant, either isotropic or orthotropic, but not temperature dependent. In that case, the formulation of the Finite Element problem is linear.

Finite Element Linear Formulation of the Heat Equation
The thermal model was obtained from the Finite Element formulation of the classical heat equation based on the Fourier phenomenological approach. The assumptions were Part 2 presents the process of extracting the thermal time constants as well as the thermal resistance with our MOR technique.
In part 3, a one-dimensional analysis was performed for single-layer and two-layer structures. A comparison between the analytical approach and the numerical one is detailed.
In part 4, we consider a two-dimensional approach that takes into account the length of the heat source.
In part 5, some simulation results are compared with the experimental data in order to validate the approach proposed in this article.

Extraction of Real Thermal Time Constants with Ritz Model Order Reduction Technique
Nonlinear Model order Reduction is a very complex task. The assumption for the following development was that thermal properties are constant, either isotropic or orthotropic, but not temperature dependent. In that case, the formulation of the Finite Element problem is linear.

Finite Element Linear Formulation of the Heat Equation
The thermal model was obtained from the Finite Element formulation of the classical heat equation based on the Fourier phenomenological approach. The assumptions were adiabatic conditions on lateral directions, uniform heat flux on the surface that represents the M.
where K the conductivity matrix and M the specific heat matrix are (n-by-n) symmetric and positive-definite matrices, T is (n-by-1) temperature vector at the mesh nodes, and F is the (n-by-1) load vector that accounts for both power generation and Dirichlet boundary conditions. n is the total number of mesh nodes. In our approach, ANSYS simulation software was solely used to extract the conductivity matrix K, the specific heat matrix M, and the load vector F using the excellent tricks proposed in [12]. Then, a reduced thermal equivalent model (τ i ; Ri) was obtained using the MOR method based on the Ritz vector approach detailed in [11,13].
The idea here was not to discuss the various MOR techniques that can be found in the literature, but to present our approach. The MOR has been extensively studied in the past years; however, to the best of our knowledge, this is the first time that it has been used to demonstrate the extraction of thermal time constants with MOR, which provides the real thermal time constants obtained by the analytical solution of the heat equation.
The main advantages of the Ritz vector method are that important response modes are not neglected and the accuracy with fewer vectors is improved compared with the use of eigenvectors. Moreover, using this method, the steady state temperature is always exacted. Indeed, the first Ritz vector is obtained by solving the steady state KT = F. The other vectors originate from orthogonal projections. When the time is infinite, their contribution decreases and tends to zero. With a pure eigenvalue decomposition, there are several major issues: - The number n of numerical equations is very large (few 100,000), so it is very difficult to compute all the eigenvalues and the eigenvectors; - The second difficulty involves selecting the main eigenvalues; - The steady state is not reached if all the eigenvalues are not considered.
The procedure for generating orthogonal Ritz vectors leads to an m-by-n projection matrix as shown in [11], where m is the number of time constants retained for the final solution. m depends on the precision proposed for the reduced model (see, for example the influence of m in results presented Figure 5). In the generation of the transformation matrix Φ m = (ϕ 1 . . . ϕ m ) of the m Ritz vectors, the relationships are verified: where T is the vector of the original temperatures (n) and p the vector of temperatures in the « reduced » system (m << n). During the projection stage, Φ m matrix is M normalized. From this transformation, the heat Equation can be written in terms of reduced coordinates with K* a (m-by-m) matrix The next stage consists in diagonalizing this system to decouple the equations:

of 13
The solution in the initial basis can be easily computed by: Thus: The selection of the output node for Tmax for example is very easy because one only needs to select a specific coordinate in T.

1-D, Single-Layer Simulation
The selection of the output node for Tmax for example is very easy because o needs to select a specific coordinate in T.  Where ρ is the density, Cp is the specific heat, K is the thermal conductivity, a the thickness of the structure. The values are reported in Table 1.   where ρ is the density, C p is the specific heat, K is the thermal conductivity, and L is the thickness of the structure. The values are reported in Table 1. Table 1. SiC Thermal properties and layer width used in the simulation.   These results highlight on one hand that more than one time constant is needed to represent the transient thermal behavior of a simple 1-D structure, and on the other hand they exhibit the difficulty and complexity of the transient approach. Furthermore, Bagnall established the analytical time constant formula for the structure presented in [7].

1-D, Single-Layer Simulation
where α is the thermal diffusivity of the layer defined by: Table 2 shows an excellent agreement between the analytical method [7] used to extract the thermal time constants with their fractional weight and our numerical method based on the Ritz vector MOR theory applied to Equation (1). This investigation proves that results determined by MOR for any structure are likely to provide the real physical values for thermal time constants and thermal resistances.    These results highlight on one hand that more than one time constant is needed to represent the transient thermal behavior of a simple 1-D structure, and on the other hand they exhibit the difficulty and complexity of the transient approach. Furthermore, Bagnall established the analytical time constant formula for the structure presented in [7].
where α is the thermal diffusivity of the layer defined by: Table 2 shows an excellent agreement between the analytical method [7] used to extract the thermal time constants with their fractional weight and our numerical method based on the Ritz vector MOR theory applied to Equation (1). This investigation proves that results determined by MOR for any structure are likely to provide the real physical values for thermal time constants and thermal resistances.  These results highlight on one hand that more than one time constant is needed to represent the transient thermal behavior of a simple 1-D structure, and on the other hand they exhibit the difficulty and complexity of the transient approach. Furthermore, Bagnall established the analytical time constant formula for the structure presented in [7].
where α is the thermal diffusivity of the layer defined by: Table 2 shows an excellent agreement between the analytical method [7] used to extract the thermal time constants with their fractional weight and our numerical method based on the Ritz vector MOR theory applied to Equation (1). This investigation proves that results determined by MOR for any structure are likely to provide the real physical values for thermal time constants and thermal resistances.  Figure 5 shows the influence of the number of the Ritz values. In the first example, we can observe that with five values, a very good accuracy was obtained from 1 ns to the steady state, while this one was reached in all cases.   Figure 5 shows the influence of the number of the Ritz values. In the first example, we can observe that with five values, a very good accuracy was obtained from 1 ns to the steady state, while this one was reached in all cases.

One-Dimension, Two-Layers Simulation
The GaN HEMTs structure may require to take into consideration more than the substrate to determine the transient response of the device. Thus, the GaN epitaxy layer was added. Figure 6 shows the structure described where t1 = 70 µm and t2 = 1.7 µm. The extraction of the different thermal time constants was performed using MOR. Parameter values for these studies are reported in Table 3. The normalized temperature is plotted Figure 7 for both structures (SiC and SiC+GaN), whereas the thermal time constant spectrum is displayed in Figure 8. In Figure 7, curves exhibit a very similar behavior with a very small error. In Figure 8, the thermal time constants have more or less the same value and the same fractional weight. This small variation is due to the small influence of the couple (GaN thickness, GaN thermal conductivity) with respect to the substrate layer.

One-Dimension, Two-Layers Simulation
The GaN HEMTs structure may require to take into consideration more than the substrate to determine the transient response of the device. Thus, the GaN epitaxy layer was added. Figure 6 shows the structure described where t 1 = 70 µm and t 2 = 1.7 µm. The extraction of the different thermal time constants was performed using MOR. Parameter values for these studies are reported in Table 3. The normalized temperature is plotted Figure 7 for both structures (SiC and SiC+GaN), whereas the thermal time constant spectrum is displayed in Figure 8. In Figure 7, curves exhibit a very similar behavior with a very small error. In Figure 8, the thermal time constants have more or less the same value and the same fractional weight. This small variation is due to the small influence of the couple (GaN thickness, GaN thermal conductivity) with respect to the substrate layer.            Figure 6 can be considered as one "equivalent" layer structure where: By making use of Equation (9): In general, for GaN HEMTs transistors, the width of the substrate layer is negligible compared with the width of the GaN layer (L 1 >> L 2 ).
Thus, using the Taylor expansion of this function, Equation (13) gives: If α 1 ·L 1 α 2 ·L 2 : The thermal time constant can be expressed as: Table 4 shows that Equations (6) and (7) found using the Taylor expansion demonstated good agreement with MOR results.

2-D, One-Layer Simulation
In GaN HEMT RF power devices, the gate length is typically between 0.1 and 1 µm. The structure used for the 2-D simulation is shown in Figure 9. With the technique used previously, a 2-D simulation using ANSYS software was performed. The K and M matrices and load vector were extracted using our MOR process. Then, the different thermal time constants and thermal resistances were determined. Figure 10 shows the comparison between the normalized time domain thermal re- With the technique used previously, a 2-D simulation using ANSYS software was performed. The K and M matrices and load vector were extracted using our MOR process. Then, the different thermal time constants and thermal resistances were determined. Figure 10 shows the comparison between the normalized time domain thermal response of an evenly distributed heat flux applied both on the full length (400 µm) and a length lg = 0.1 µm. The 2-D solution for lg = 0.1 µm exhibited a faster response than for lg = 400 µm. This behavior reveals the existence of more thermal time constants and a broadening of the time constant spectrum. With the technique used previously, a 2-D simulation using ANSYS software was performed. The K and M matrices and load vector were extracted using our MOR process. Then, the different thermal time constants and thermal resistances were determined. Figure 10 shows the comparison between the normalized time domain thermal response of an evenly distributed heat flux applied both on the full length (400 µm) and a length lg = 0.1 µm. The 2-D solution for lg = 0.1 µm exhibited a faster response than for lg = 400 µm. This behavior reveals the existence of more thermal time constants and a broadening of the time constant spectrum.    From the previous observations, complementary 2-D studies were performed for a heating region length "lg" varying from 0.1 µm to 1 µm. Figure 12 shows that all the extracted thermal time constants were independent of the gate length except for the shortest one τ21. Regarding the thermal resistances ( Figure  13), the only non-constant resistance was the one corresponding to the shortest thermal time constant (τ21) i.e., the one corresponding to the close vicinity of the heat source. In addition, another interesting result relies on the fact that the shortest thermal capacitance was linear in function of "lg" while the thermal resistance was linear with log (lg) [14].  From the previous observations, complementary 2-D studies were performed for a heating region length "lg" varying from 0.1 µm to 1 µm. Figure 12 shows that all the extracted thermal time constants were independent of the gate length except for the shortest one τ 21 . Regarding the thermal resistances (Figure 13), the only non-constant resistance was the one corresponding to the shortest thermal time constant (τ 21 ) i.e., the one corresponding to the close vicinity of the heat source. In addition, another interesting result relies on the fact that the shortest thermal capacitance was linear in function of "lg" while the thermal resistance was linear with log (lg) [14]. This observation on the trend of the thermal resistances is fully consistent with Darwish and al. [6], where it is explained that only a part of the resistance depends on log (lg). heating region length "lg" varying from 0.1 µm to 1 µm. Figure 12 shows that all the extracted thermal time constants were independe the gate length except for the shortest one τ21. Regarding the thermal resistances (Fi 13), the only non-constant resistance was the one corresponding to the shortest the time constant (τ21) i.e., the one corresponding to the close vicinity of the heat sourc addition, another interesting result relies on the fact that the shortest thermal capacit was linear in function of "lg" while the thermal resistance was linear with log (lg) This observation on the trend of the thermal resistances is fully consistent with Dar and al. [6], where it is explained that only a part of the resistance depends on log (lg) Then, the temperature transient response for this structure can be modeled usi Foster-type circuit with 20 constant-value RC blocks with one variable-paramete block.

2-D, Two-Layer Simulation
As previously mentioned, GaN HEMTs require more than the substrate to be con ered in determining the transient response of the device; thus, approximating the structure of GaN transistors, a GaN epitaxy layer was added with a gate length betw 0.1 and 1 µm. Figure 14 shows the structure used for the 2-D two-layer simulations the first simulation, lg = 0.1 µm, L2 = 1.7 µm, and L2 = 70 µm. Then, the temperature transient response for this structure can be modeled using a Foster-type circuit with 20 constant-value RC blocks with one variable-parameter RC block.

2-D, Two-Layer Simulation
As previously mentioned, GaN HEMTs require more than the substrate to be considered in determining the transient response of the device; thus, approximating the real structure of GaN transistors, a GaN epitaxy layer was added with a gate length between 0.1 and 1 µm. Figure 14 shows the structure used for the 2-D two-layer simulations. For the first simulation, lg = 0.1 µm, L 2 = 1.7 µm, and L 2 = 70 µm.
As previously mentioned, GaN HEMTs require more than the substrate to be ered in determining the transient response of the device; thus, approximating t structure of GaN transistors, a GaN epitaxy layer was added with a gate length b 0.1 and 1 µm. Figure 14 shows the structure used for the 2-D two-layer simulatio the first simulation, lg = 0.1 µm, L2 = 1.7 µm, and L2 = 70 µm.  Figure 15 shows that the GaN layer epitaxy had a larger impact on the thermal ent response than for large gate length, and due to the to the small GaN thickness a pared with the substrate layer, the impact on the thermal time constant was neglig shown in Figure 16.  Figure 15 shows that the GaN layer epitaxy had a larger impact on the thermal transient response than for large gate length, and due to the to the small GaN thickness as compared with the substrate layer, the impact on the thermal time constant was negligible, as shown in Figure 16.

Calibration between Simulation Results and Measurements
These results presented in Figure 17 concern the steady state for an eight-finger GaN

Calibration between Simulation Results and Measurements
These results presented in Figure 17 concern the steady state for an eight-finger GaN HEMT with gate length of 0.15 µm and gate width of 50 µm. They allow us to validate the FEM approach in terms of thermal resistance as well as the geometry and mesh. Figure 16. Thermal time constant spectrum for 2-D lg = 0.1 µm structure: comparison betwe layer and 2-layer structure.

Calibration between Simulation Results and Measurements
These results presented in Figure 17 concern the steady state for an eight-fin HEMT with gate length of 0.15 um and gate width of 50 um. They allow us to val FEM approach in terms of thermal resistance as well as the geometry and mesh.  A thermoreflectance measurement of the temperature is presented in [15,16]. The maximum of temperature variation from ambient for this device was 22 • C during a dissipated power of 0.7 W, leading to a thermal resistance of 31.4 • C/W.
The Finite Element simulation presented in Figure 18 was performed by considering symmetries for 1 4 of the structure. Adiabatic lateral boundary conditions, a baseplate of 25 • C, and a total dissipated power of 4 W were considered for the simulation. Thus, the extracted thermal resistance obtained was (148.6 − 25)/4, i.e., 30.9 • C/W.

Electronics 2021, 10, x FOR PEER REVIEW
A thermoreflectance measurement of the temperature is presented in [15, maximum of temperature variation from ambient for this device was 22 °C duri sipated power of 0.7 W, leading to a thermal resistance of 31.4 °C/W.
The Finite Element simulation presented in Figure 18 was performed by con symmetries for ¼ of the structure. Adiabatic lateral boundary conditions, a bas 25 °C, and a total dissipated power of 4 W were considered for the simulation. T extracted thermal resistance obtained was (148.6 − 25)/4, i.e., 30.9 °C/W.

Conclusions
In this paper, we provided elements for the understanding of thermal time c in GaN HEMTs. The main goal of this work was not to propose the best physi method to understand the submicronic GaN HEMT device properties, but instea approach based on a well-known Fourier theory and Finite Element Method tool are commonly used in labs and the electronics industry.
In summary, we intended to show: that there are several thermal time constants in a HEMT, even with a 1-D, sing the influence of the GaN layer on the thermal time constants; -the evolution of the thermal time constants with the dimensions of the dev in particular with the gate length, which is one of the key points to address frequency applications.

Conclusions
In this paper, we provided elements for the understanding of thermal time constants in GaN HEMTs. The main goal of this work was not to propose the best physics-based method to understand the submicronic GaN HEMT device properties, but instead a new approach based on a well-known Fourier theory and Finite Element Method tools, which are commonly used in labs and the electronics industry.
In summary, we intended to show: that there are several thermal time constants in a HEMT, even with a 1-D, single layer; -the influence of the GaN layer on the thermal time constants;