Various simulations of a GaN HEMT are used to study the behaviors of two different energy-transport models: the Fermi kinetics transport model and a hydrodynamics transport model as it is implemented in the device simulator Sentaurus from Synopsys. The electron transport and heat flow equations of the respective solvers are described in detail. The differences in the description of electron flux and the discretization methods are highlighted. Next, the transport models are applied to the same simulated device structure using identical meshes, boundary conditions, and material parameters. Static simulations show the numerical convergence of Fermi kinetics to be consistently quadratic or faster, whereas the hydrodynamic model is often sub-quadratic. Further comparisons of large-signal transient simulations reveal the hydrodynamic model produces certain anomalous electron ensemble behaviors within the transistor structure. The fundamentally different electron dynamics produced by the two models suggest an underlying cause for their different numerical convergence characteristics.

## I. INTRODUCTION

Recent decades have seen increased interest in GaN devices for power and radio frequency (RF) electronics.^{1,2} Numerical simulation is a necessary step for the design and analysis of these devices. Robust simulation models rely on producing predictive results. Accurate simulation of these devices must capture the important physical processes involved and the carrier transport mechanisms. While it is possible to obtain reasonable results using ad hoc solutions despite incorrectly modeling the physics, these results are not predictive and do not aid in device design.^{3–5} Furthermore, these numerical methods may be problematic for large-signal transient simulations, which are particularly important for RF power applications. This paper seeks to explore some numerical behaviors of established charge transport solvers for RF GaN HEMTs.

Charge transport in a semiconductor is described by the Boltzmann transport equation (BTE). The BTE is an integrodifferential equation that describes the evolution of the non-equilibrium distribution function of electrons under the influence of an external electric field. The BTE is given as^{6}

where $f$ is the distribution function, $F$ is the external force, $v$ is the group velocity, $k$ is the wavevector in the reciprocal space, $r$ is the position vector in the real space, $t$ is the time, and the term on the right-hand side represents the rate of change of $f$ due to the collisions between the carriers.

Monte Carlo solvers are a useful statistical method to solve the BTE. The Monte Carlo technique generally employed is the ensemble Monte Carlo (EMC) method.^{7} The accurate results produced by this method come at the cost of immense computational load. To reduce the computational load encountered by the EMC, the hybrid cellular automaton/Monte Carlo (CA/MC) algorithm was developed. The CA/MC scheme attempts to speed-up the simulation by pre-calculating the scattering rates and storing them in a table. It was shown that the CA/MC method required about 0.2–0.5 CPU-second per iteration, up to 25 times faster than the traditional EMC.^{8} However, this speed improvement offered by the hybrid scheme is still not comparable to the shorter simulation time taken by deterministic solvers.

Deterministic solvers are an alternative approach that does not rely on random sampling to solve the BTE. While the deterministic methods are generally much faster than the Monte Carlo methods, they rely on physical approximations to the BTE to obtain the solutions.^{9} Deterministic solvers convert the BTE into an equivalent infinite series of differential equations in phase space by applying the method of moments over the reciprocal space. Truncating the series after a finite number of moments typically results in one solution variable more than the available equations. Therefore, a closure relation is required to make the system of equations tractable.^{10}

Scharfetter and Gummel^{11} developed one of the first deterministic Boltzmann solvers. They proposed a robust discretization of the drift-diffusion (DD) model with physical approximations. However, as the device dimensions were scaled, the DD model failed to account for the large electric fields and the hot-carrier effects that arise, e.g., in FET channels. These effects are typically accounted for in the DD model by empirically modifying the mobility.^{12}

The DD model was extended by Stratton^{13} and Blotekjaer^{14} to include a balance equation for the average energy of the carriers. These models are generally referred to as energy-transport models, and the hydrodynamic model is a prime example.^{10} Hydrodynamic models typically use the parabolic band approximation to describe the kinetic energy vs momentum dispersion relation. This approximation fails to hold when the electric fields in the device reach large magnitudes. To account for this, the model relies on empirical field-dependent and temperature-dependent mobility models. The commonly used closure relation in hydrodynamic models is the description of heat flow using Fourier’s law and an approximate thermal conductivity.

Synopsys Sentaurus is an industry-leading suite of semiconductor device and process simulation tools. Sentaurus has been widely used to simulate devices for various applications including power electronics,^{15,16} RF electronics,^{17} and CMOS technology.^{18,19} Sentaurus provides a wide range of simulation packages to handle the vast application space and support a streamlined simulation workflow. The range of packages available includes Sentaurus Process for process simulation, Sentaurus Device Monte Carlo, Sentaurus Device Electromagnetic Wave Solver, and, the package under study in this work, Sentaurus Device (SDevice).^{20} The history of SDevice can be traced back to the development of DESSIS in 1992. DESSIS was developed by merging HFIELDS from the University of Bologna,^{21} SIMUL from ETH Zurich,^{22} BONSIM from Bosch, and in collaboration with ST Microelectronics.^{23} In 1993, Integrated Systems Engineering (ISE) was founded and assumed responsibility for the simulator. ISE merged with Synopsys in 2004, and DESSIS became the basis of SDevice.^{24} SDevice supports several carrier transport models that can be used depending on the device under investigation. The hydrodynamic model, implemented in SDevice, will be investigated in this work.

FKT is also an energy transport based BTE solver that employs the method of moments but uses an alternative treatment of the electronic heat flow. It substitutes the idea of thermal conductivity with the heat capacity of an ideal Fermi gas to be used as the closure relation.^{25} It has been demonstrated that FKT, unlike the hydrodynamic model, could incorporate the electronic band structure^{26} and full-wave EM.^{27} FKT has also demonstrated large-signal RF simulations^{28} and defect dynamics simulations.^{29} More details on the solver are reported in the literature.^{30–32}

This paper examines the physical applicability of both the FKT and the hydrodynamic model employed by Sentaurus, hereafter referred to as commercial hydrodynamic transport (CHT) solver. First, the system of equations and the discretization scheme employed by the solvers are presented in Sec. II. Section III presents the GaN HEMT structure and the simulation parameters used in the respective simulators. Next, the thermal equilibrium simulations are performed and the results are compared. The output characteristics and the electron temperatures of the HEMT structure, produced by the non-equilibrium simulations, are compared. Further, the differences in the simulation frameworks are investigated by analyzing the numerical convergence behavior after varying the simulation setup to ensure similar simulation conditions for the purpose of comparison. Finally, the effect of varying the polarization sheet charge density at the AlGaN/GaN interface on the convergence characteristics and the transient simulation behavior of the simulators are investigated.

## II. COMPUTATIONAL BACKGROUND

The semiconductor device simulators considered here are both of the energy-transport variety. As such, they have a great deal in common. For example, they both tessellate a problem domain with a suitable mesh comprised of triangles for 2D simulations or tetrahedra for 3D simulations. This mesh conforms to the simulated device’s geometry and is refined appropriately in regions where solution variables change rapidly, such as material interfaces and metallurgical junctions. Solution variables assigned to each mesh point include an electrostatic potential to account for voltage drops and the corresponding electric fields throughout the device. To each mesh point is also assigned a Fermi–Dirac distribution for an electron gas occupying the semiconductor’s conduction band to account for currents flowing through the device under the influence of the internal fields. Integrals of the distribution function over the density of conduction band states can provide the particle and energy densities of mobile electrons throughout the device. These densities can be treated as additional solution variables at each mesh point or they can be equivalently associated with alternative solution variables such as the Fermi distribution’s chemical potential (Fermi level) and temperature.

Both device simulators calculate the solution variables by solving Poisson’s equation, the electron continuity equation, and the conservation of electron energy at each mesh point using the box-integration method.^{9} The electrostatic potentials are approximated as piecewise linear along the edges connecting neighboring mesh points, and the electron fluxes are approximated as piecewise constant between mesh points. The electron particle and kinetic energy fluxes can be derived from the 1st and 3rd moments of the BTE, respectively.^{6} Both simulators represent these fluxes in a discrete form on the device mesh using the powerful Scharfetter–Gummel method.^{9,11} Because the mobile electrons are treated as ideal Fermi gases, both simulators also account for the electronic heat that flows between neighboring electron distributions.

Although the simulators use these common methodologies, there are some differences in how they are implemented that lead to significant differences in the simulation results and numerical convergence characteristics. Principal differences include the precise forms of the electron fluxes obtained from moments of the BTE, the way Scharfetter–Gummel discretization is applied to these fluxes, and the way electronic heat flow is treated. The following details these differences along with the associated simulation results and numerical behaviors. For simplicity, effective mass and constant mobility approximations will be used.

### A. CHT simulator

For the sake of simplicity, we will limit ourselves to the consideration of only electrons in steady-state conditions. Additionally, the lattice temperature, $TL$ is set to 300 K to disentangle the effects of electron heating and lattice heating, resolve convergence issues, and enable a fair comparison of CHT against FKT. The electron current density, modeled using CHT, is given as^{20}

The first term in the above equation accounts for the spatial variation of electrostatic potential, while the remainder terms account for the gradient of concentration and the carrier temperature gradient. $\mu n$ is the electron mobility, $n$ is the electron concentration, and $Tn$ denotes the temperature of electrons. The thermal diffusion constant $fntd$ is set equal to zero, which corresponds to the Stratton model.^{33} The parameters $\gamma n$ and $\lambda n$ account for the Fermi–Dirac statistics and are given as

$Fj(x)$ is the Blakemore Fermi–Dirac integral of order-$j$. The argument $\eta n$ is the reduced Fermi energy and is defined as $\eta n=(Fn\u2212EC)/kTn$. The effective conduction band density of states in the semiconductor is specified by $NC$. CHT provides an option, in lieu of $Tn$, to use a blended temperature in (2) according to $Tn\u2192gTn+(1\u2212g)TL$, where $TL$ is the lattice temperature and $g$ is a dimensionless parameter between zero and unity. For all CHT simulations reported in this paper, $g$ is set equal to unity.

The energy balance equation is given as

where $Wn$ is the electron energy density and $S\u2192n$ is the energy flux. $Wn$ comprises a thermal and a convective term; however, CHT implements $Wn$ without the convective term. That is, $Wn=nwn=n(3kTn/2)$. The term $J\u2192n\u22c5\u2207\u2192ECq$ represents the energy delivered to the carriers by the field per unit time, while $dWndt|coll.$ is the net energy gained from the lattice via phonon interactions. $S\u2192n$ is given as

and

Here, we use the default parameter values for $rn$ and $fnhf$ as listed in Table I. These parameters allow us to change the relative contributions of the convective and diffusive terms in the energy flux. With the default parameters, the prefactor for the diffusive term becomes

Parameter . | Default value . | Meaning . |
---|---|---|

r_{n} | 0.6 | Parameter in energy flux |

$fnhf$ | 1 | Parameter in energy flux |

g | 1 | Parameter to select carrier temperature [T_{n} → gT_{n} + (1 − g T_{L})] |

$fntd$ | 0 | Thermal diffusion coefficient |

n_{min} | 10^{3} cm^{−3} | Parameter to improve numeric stability |

Parameter . | Default value . | Meaning . |
---|---|---|

r_{n} | 0.6 | Parameter in energy flux |

$fnhf$ | 1 | Parameter in energy flux |

g | 1 | Parameter to select carrier temperature [T_{n} → gT_{n} + (1 − g T_{L})] |

$fntd$ | 0 | Thermal diffusion coefficient |

n_{min} | 10^{3} cm^{−3} | Parameter to improve numeric stability |

The collision term in the energy balance is given as

$Wn0$ is the equilibrium electron energy density, and $\tau e$ is the *energy relaxation time* for electrons. The parameter $\zeta n$ is adjusted to improve the stability of the numerical algorithm and is given as

with $nmin$ and $n0$ are adjustable density parameters to aid convergence. Since $Tn\u2265TL$ for our device, $n0$ does not impact $\zeta n$. The default value of $nmin$ is set to $103cm\u22123$, which leads to $\zeta n\u22481$ for all test cases reported here. While $\tau e$ is energy-dependent and would lead to a modification of the diffusive energy flux, here we treat it as an energy-independent parameter. However, CHT provides empirical methods, such as spline approximation of energy relaxation time over energy, to account for $\tau e$’s energy dependence. In this study, we vary $\tau e$ in the fs to ps range to analyze its impact on *I*–*V* curves and electron heating. Additionally, we elucidate that $\tau e$’s impact on *I*–*V* convergence in CHT. Detailed results are discussed in Sec. III.

To discretize the electron fluxes, CHT uses the ingenious Scharfettel–Gummel method.^{9,11} In its original form, this method assumed the electric field $E=(EC,1\u2212EC,0)/(qL)$ and the electron flux $Jn$ were constant on a mesh edge of length $L$ connecting mesh points 0 and 1. With these assumptions, the phenomenological drift-diffusion equation was treated as a first-order differential equation for the electron density $n$. Solving this differential equation for variations in $n$ along the mesh edge then produced a discrete form for the electron flux, given in (11), that proved to be numerically robust,

where $i$ and $j$ are the major mesh points, $\mu n$ is the electron mobility (assumed constant between the mesh points), $\u2113ij$ is the length of the mesh edge, $B(x)=x/(ex\u22121)$ is the Bernoulli function, $ui\u2212uj=qE\u2113ij/kTn$, and $E$ is constant between the mesh points. For high electron densities, CHT may also add some additional terms to (11) to account for Fermi–Dirac degeneracy.

### B. Fermi kinetics transport

Because FKT, like CHT, assumes an energy-dependent electron distribution function at points in the simulation domain, electron fluxes are computed from moments of the BTE. This method of moments is a powerful technique for inferring information about the part of the distribution function odd in momentum-space, such as net particle flux, from the even part of the distribution function. For example, drift-diffusion flux, which is often presented phenomenologically, can be derived from the first moment of the BTE assuming a Maxwell–Boltzmann distribution function, constant electron temperature, and constant mobility.^{6} In general, FKT uses electron fluxes derived for a Fermi–Dirac distribution, non-constant temperature, non-parabolic bands, and an energy-dependent mobility computed from electron scattering rates.^{26} To simplify comparisons with CHT, the electron current density from the first moment assuming parabolic bands and constant mobility is given by

where $Fj(\eta n)$ is the integral over energy $E$ of the electron Fermi function times $Ej$ and $Fj\u2032(\eta n)$ is its derivative with respect to $\eta n$. This integral is proportional to the Fermi integrals in (4)$Fj(\eta n)=Fj(\eta n)/\Gamma (j+1)$. The electron internal/kinetic energy flux is likewise derived from the third moment of the BTE,

Like CHT, FKT uses the Scharfetter–Gummel discretization, but it treats (12) as a first-order differential equation for the quantity $(kTn)5/2F3/2$. Solving this differential equation produces a discretized electron current density on the mesh edge given by

where the Bernoulli function $B(\xi n)=\xi n/[exp(\xi n)\u22121]$ and $\Delta (x)\u2261x1\u2212x0$. The same procedure is applied to (13) to obtain the Scharfetter–Gummel discretization of the electron kinetic energy flux,

Although FKT particle and kinetic energy fluxes as well as the flux discretization schemes differ somewhat from the corresponding quantities in CHT, perhaps the most significant departure is the treatment of electronic heat flow. It is evident from (6) that CHT treats electronic heat flow with Fourier’s law using an electronic thermal conductivity related to the electrical conductivity. FKT, on the other hand, exploits the kinetics of ideal Fermi gases.^{25} Because electron ensembles at points inside the simulation are assumed to have an energy-dependent distribution, i.e., independent of the direction of degenerate electron momentum vectors $k$, the electron ensembles are symmetric in $k$-space. This is the defining feature of an ideal gas. When simulated electron fluxes move along a mesh edge from an initial mesh point with Fermi level $(Fn\u2212EC)0$ and temperature $Tn,0$ to a final mesh point with Fermi level $(Fn\u2212EC)1$ and temperature $Tn,1$, the electrons must thermalize in order to change their temperature from $Tn,0$ to $Tn,1$. This thermalization process is achieved through electron–electron scattering, and it is the source of electronic heat flow between the electron ensembles/gases. The amount of heat $H$ needed to change the temperature of an ideal Fermi gas from initial temperature $Ti$ to final temperature $Tf$ is precisely defined by the thermodynamic identity,^{34}

Determining the net electronic heat flux on the mesh edge requires the rate electrons change from $Tn,0$ to $Tn,1$ and the rate electrons change from $Tn,1$ to $Tn,0$. These rates are simply $J0\u21921$ and $J1\u21920$ from (14), respectively. Therefore, the electronic heat flux flowing from point 0 to point 1 is given by

the electronic heat flux flowing from point 1 to point 0 is given by

and the net electronic heat flux along the edge is $Snheat=S0\u21921heat\u2212S1\u21920heat$. Combining this heat flow with the kinetic energy flux gives the total electron energy flux along the edge connecting mesh points 0 and 1, $Sn=Snkin+Snheat$.

FKT uses the same Joule heat CHT uses in (5). FKT typically computes the collision operator that determines the rate hot electrons transfer energy to the semiconductor lattice by integrating the phonon scattering rates over electron energy. However, for the sake of the comparisons with CHT presented here, FKT was modified to use (9).

## III. COMPARISONS OF SIMULATION RESULTS

The GaN HEMT structure simulated in FKT and CHT is illustrated in Fig. 1. The n$+$ GaN source/drain contacts are doped with $1019cm\u22123$ donors and complete dopant ionization is assumed. The GaN region is semi-insulating containing $9\xd71015cm\u22123$ fully ionized acceptor atoms compensated with $1016cm\u22123$ donor-like traps located at an energy level 0.6 eV below $EC$ and capture cross section of $10\u221219m2$. The Schottky barrier height at the gate contact boundary condition is assumed to be 1.5 eV, while the source and drain contacts are ohmic boundary conditions. For both simulators, AlGaN and SiC are treated as insulators, which means that carrier transport equations are not solved in these regions. All the simulations for both the simulators were run on the same mesh, as shown in Fig. 2.

### A. Thermal equilibrium

The thermal equilibrium simulation does not depend on the transport model used. Therefore, one would expect identical equilibrium results from the two simulators, given the mesh, the material and simulation parameters for the GaN HEMT are kept alike. Keeping the aforementioned conditions identical, we calibrate the two simulators and verify that they produce similar results. Table II shows the material parameters that are kept uniform. Figures 3–5 show the z-direction band diagram, z-direction electron density, and the x-direction electron density plots, respectively. As expected, excellent agreement between the simulators is obtained.

Material . | Dielectric constant . | Electron affinity (eV) . | Bandgap (eV) . | Electron mobility (cm^{2}/V s)
. | Effective mass (m_{0})^{a}
. |
---|---|---|---|---|---|

AlGaN | 8.9 | 3.55 | 4.14 | … | … |

GaN | 8.9 | 4.1 | 3.457 | 440 | 0.2 |

SiC | 9.66 | … | … | … | … |

Material . | Dielectric constant . | Electron affinity (eV) . | Bandgap (eV) . | Electron mobility (cm^{2}/V s)
. | Effective mass (m_{0})^{a}
. |
---|---|---|---|---|---|

AlGaN | 8.9 | 3.55 | 4.14 | … | … |

GaN | 8.9 | 4.1 | 3.457 | 440 | 0.2 |

SiC | 9.66 | … | … | … | … |

^{a}

m_{0} is the free electron mass

### B. Non-equilibrium simulations

First, a drift-diffusion simulation was run using both FKT and CHT. Once again, the material parameters used for the simulations and shown in Table II are kept identical. The voltage step size is fixed at 0.1 V for both the simulators. As seen in Fig. 6, both the simulators give nearly identical drain current. The minor differences in the obtained result can be attributed to the difference in the representation of the electron flux in CHT, [Eq. (2)], and the differences in the Scharfetter–Gummel discretization in the respective simulators.

Next, an energy-transport simulation was run using the two simulators. The parabolic band approximation was used to model the dispersion relation with an electron effective mass of $0.2m0$ ($m0$ is the free electron mass). The energy relaxation time ($\tau e$) is $0.2$ ps.^{35,36} The CHT-specific hydrodynamic simulation parameters used are shown in Table I. For the FKT simulation, the voltage step size is kept fixed and unchanged at 0.1 V. On the other hand, it was observed that an adaptive step size was required for the convergence of the CHT simulation. The maximum allowed voltage step size is set to 0.1 V, and the minimum allowed (initial) step size is lowered to 2.5 mV.

The $ID\u2212VDS$ characteristic and the electron temperature profile in the channel produced by the simulators is shown in Figs. 7 and 8, respectively. It would seem unexpected that despite the higher electron temperature obtained by CHT, the drain current produced is also marginally higher than that produced by FKT. This can be explained by the fact that both the simulators are using a constant mobility. As a result, the scattering effects associated with high electron temperatures are not considered. Thus, in this case, higher electron temperatures signify higher carrier velocity with constant mobility, and consequently higher current.

As noted above, while the output characteristic obtained are quite similar, the electron temperature profile in the channel shows a marked difference between the two simulators. The electron temperatures produced by FKT are considerably lower than those from CHT, when $\tau e=0.2$ ps is chosen for both simulators. This is because of the significantly higher electronic heat flow into the ohmic contacts through the degenerately doped source/drain regions in FKT. More physically realistic simulations should include lattice heating, and FKT would transfer electron energy entering the ohmic mesh points directly into the semiconductor crystal, leading to lattice temperature hot spots at ohmic contact mesh points, particularly those on the drain contact. However, the simplified simulations presented here neglect lattice heating in order to highlight electron energy transport, and that electronic energy flowing into ohmic mesh points is removed instantly from the simulation domain.

We next analyze the convergence of the nonlinear solvers to explore other differences of the transport simulators. The outputs of the Newton solvers were stored to disk and post processed using custom Python codes. In particular, the Newton solver residuals were used to analyze the convergence of the transport simulators. These residuals are a function of Newton iteration and are computed at each quiescent gate and drain voltage. We define the Newton solver residual computed at the $n$th Newton iteration for a specific gate and drain quiescent bias as $en(VGS,VDS)$.

An example of Newton solver residuals is illustrated in Fig. 9. Here, residual values computed in the CHT device and FKT simulations are reported in the left and right of the figure, respectively. These residuals were computed at $VGS=\u22124$ V over the quiescent drain voltage sweep from $VDS=0$ to 10 V. The FKT simulations exhibit excellent convergence characteristics for all quiescent drain voltage biases. However, as can be seen in Fig. 9, there are some instances of slow convergence exhibited by the CHT device simulator. For example, there is one Newton solver instance which takes ten iterations to complete, and the final residual value is well above the typical value of $<10\u221210$.

A quantitative investigation of the convergence properties of the transport solvers is next presented. In particular, we analyze the rate of convergence (ROC), $\alpha $, of the numerical solvers. The ROC is calculated by^{37}

The calculated ROCs processed from the CHT and FKT Newton solver residuals are illustrated in Fig. 10. These results indicate that FKT exhibits at least quadratic convergence ($\alpha =2$) for all the drain biases. It should be noted that there exists some FKT ROC values which fall slightly below $\alpha =2$. These values, on the order of $\alpha =1.8$, are the result of numerical precision of the solver. For example, the set of Newton solver residuals which contain a value at the 7th iteration shed light on the impact of numerical precision of the solver. The slope of the residual curve changes between iterations 6 and 7 due to hitting the numerical noise floor in the solver. This change in slope affects the estimated ROC value and would not occur if there was not a numerical noise floor in the calculations. The ROC of CHT is also quadratic at several quiescent bias solutions. However, the ROC of the CHT simulation is generally super-linear.

As noted above, the two simulators exhibit different convergence characteristics, with FKT generally displaying quadratic convergence, compared to CHT showing super-linear convergence. In addition, the simulators also produced widely different electron temperatures in the channel. It would be interesting to determine if this difference in the convergence behavior can be attributed to the much higher electron temperatures obtained from the CHT simulations. This can be verified by running the simulations to control for the temperature difference. Therefore, we analyze the convergence characteristics by changing the simulation setup so as to obtain similar peak electron temperatures.

We investigate two methods to obtain similar channel electron temperatures from the two simulators. First, we try to reduce the electron temperatures produced by CHT to the approximated range of temperatures yielded by the unaltered FKT simulation. For that purpose, we reduce the $\tau e$ used in CHT to an artificially low value. Reducing the energy relaxation time allows the electrons to rapidly lose their energy to the lattice, thereby reducing the electron temperature. For the second method, we try to increase the drain bias in the FKT simulation, keeping the CHT setup unchanged. The exceedingly high fields produced should result in hot electrons with temperatures comparable to CHT.

Accordingly, we first perform the CHT simulation by reducing the $\tau e$ to $5.0$ fs, while keeping the $\tau e$ in FKT unchanged at $0.2$ ps. Figure 11 shows that the peak electron temperature in the channel is in the same range of magnitude as FKT. The convergence characteristic for the CHT simulation is shown in Fig. 12 (left). We see that the ROC is generally super-linear and quadratic or above for some of the gate biases. Therefore, reducing the electron temperature did not significantly change the CHT convergence rate, and it still shows a slower ROC compared to FKT.

In addition to simulations that lowered electron temperature in CHT to approximate the original FKT temperatures in Fig. 8, further FKT simulations also increased electron temperatures to approximate the original CHT temperatures in Fig. 8. Because $Tn$ was reduced in CHT simulations by decreasing $\tau e$ far below 0.2 ps, it might seem natural to increase $\tau e$ far above 0.2 ps in FKT to produce much higher $Tn$. However, when $\tau e$ increases excessively, the electronic heat flow to the ohmic contacts in FKT limits further increase in channel electron temperature. A more effective way to heat FKT electrons to temperatures comparable to the CHT temperatures in Fig. 8 is to simply drive the FKT simulations to higher drain voltages. Thus, the FKT simulation was run up to $VDS=50$ V at $VGS=0$ V. Figure 11 shows that the peak electron temperature obtained for this setup is comparable to that obtained from the original CHT simulation. We observe that FKT continues to show excellent quadratic convergence, as seen in Fig. 12 (right). Despite increasing the electron temperature by applying high fields, FKT continues to show quadratic convergence compared to the super-linear convergence of CHT.

Finally, changing the polarization sheet charge density ($\sigma PZ$) at the AlGaN/GaN interface changes the field at the interface. The large electric fields that arise place further burden on the solvers and affect the convergence behavior. Therefore, a series of simulations were run by increasing $\sigma PZ$, with $\tau e=0.2$ ps for both the simulators, to characterize the convergence behavior. As summarized in Table III, FKT displays greater than quadratic ROC on average for all values of $\sigma PZ$, while CHT shows a super-linear ROC that decreases as $\sigma PZ$ is increased. Furthermore, the CHT simulations failed to converge for all ($VGS$, $VDS$), when $\sigma PZ\u226515\xd71012cm\u22122$, which is indicated as “N/A” in Table III. Here, mean${\alpha (VGS,VDS)}$ is the average ROC of all ($VGS$, $VDS$) at which the simulations converged, where $VGS\u2208[\u22124,0]$ V and $VDS\u2208[0,10]$ V, as shown in Figs. 10 and 12.

. | mean{α(V_{GS}, V_{DS})}
. | |
---|---|---|

σ_{PZ} (×10^{12} cm^{−2})
. | CHT . | FKT . |

9 | 1.856 | 2.399 |

10 | 1.809 | 2.387 |

11 | 1.044 | 2.392 |

12 | 1.484 | 2.409 |

13 | 1.510 | 2.403 |

14 | 1.560 | 2.397 |

15 | N/A | 2.391 |

16 | N/A | 2.387 |

17 | N/A | 2.384 |

18 | N/A | 2.388 |

. | mean{α(V_{GS}, V_{DS})}
. | |
---|---|---|

σ_{PZ} (×10^{12} cm^{−2})
. | CHT . | FKT . |

9 | 1.856 | 2.399 |

10 | 1.809 | 2.387 |

11 | 1.044 | 2.392 |

12 | 1.484 | 2.409 |

13 | 1.510 | 2.403 |

14 | 1.560 | 2.397 |

15 | N/A | 2.391 |

16 | N/A | 2.387 |

17 | N/A | 2.384 |

18 | N/A | 2.388 |

### C. Time-dependent simulations

Because the above results suggest the different CHT and FKT models may have led to different numerical convergence behavior, transient simulations were performed to test for additional differences in behavior. Starting from the thermal equilibrium initial condition, the drain voltage $VDS$ was ramped from 0 to 1 V in 1 ns, and the loss of electron energy to the semiconductor crystal lattice was set to zero in order to highlight differences in real-space energy transport. The electron temperatures across the device computed by CHT and FKT are shown in animated Fig. 13 (Multimedia view) (left) and (right), respectively. As expected from the previous static results, the electrons heat up to higher temperatures in the CHT calculation. However, the peak CHT temperatures are relatively distant from the gate, appearing underneath the drain contact. The peak FKT temperatures, on the other hand, appear at the corner of the gate closer to the drain. Because most of the $VDS$ voltage drop occurs at this corner of the gate, it is unclear why the peak CHT electron temperature would not also occur in this location. In addition to locations of peak electron temperature, the transient simulations also produced other differences that become apparent when focusing on the shorter time scale $0<t<0.15$ ns. The animated Fig. 14 (Multimedia view) (left) and (right) show these dynamic electron temperatures across the device as produced by the CHT and FKT models, respectively. As in the case of the longer time scales, FKT shows electron temperatures increasing monotonically in time with the peak located at the drain-side corner of the gate. The CHT model again produces a peak electron temperature occurring away from the gate and underneath the ohmic drain contact, but it also shows electron temperatures in the GaN buffer region under the gate decreasing significantly below the 300 K ambient temperature. As these electron dynamics occur at constant volume in the absence of physical work, it is not immediately clear what physical mechanism could cause electron temperatures inside the transistor to fall below the ambient temperature. The fundamentally different physical behaviors revealed by these dynamically evolving transient simulations may help explain the different numerical convergence characteristics previously observed for the CHT and FKT static simulations.

## IV. CONCLUSION

Deterministic Boltzmann solvers make use of the moments of the BTE to model the charge transport in semiconductors. Hydrodynamic model and FKT are both deterministic solvers that belong to this class of energy-transport models. They differ in the way they treat the electronic heat flow. While hydrodynamic models utilize an electron thermal conductivity, FKT employs the heat capacity of an ideal Fermi gas.

In this work, we compared the hydrodynamic model implemented in CHT with the FKT model. A GaN HEMT structure was used to demonstrate the simulators. It was observed that both the simulators gave reasonably similar non-equilibrium output characteristics. However, the electron temperatures in the channel, close to the AlGaN/GaN interface, were widely different. The rate of convergence of the solvers were also different, with FKT showing quadratic convergence or higher, and CHT showing super-linear convergence.

Next, we investigated if the difference in the convergence characteristics can be attributed to the different electron temperatures obtained. We employed the two methods to ensure that the two simulators produce similar temperatures. First, we reduced the energy relaxation time used in CHT simulation, keeping the FKT setup unchanged. Second, we increased the drain bias for the FKT simulation, producing hot electrons with temperature comparable to the unchanged CHT simulation. We then compared the convergence behavior of FKT and CHT in each of the above mentioned cases. As before, FKT continued to show greater than quadratic rate of convergence, and CHT showed super-linear convergence. We can conclude that the different convergence behavior is due to the difference in the transport models and the description of the electronic heat flow rather than the difference in the electron temperatures produced. Additionally, simulations were run after changing the polarization sheet charge density at the AlGaN/GaN interface. The results show that FKT exhibited rapid quadratic convergence for all values chosen, while CHT failed to converge for large values of the polarization charge.

To further explore the possible underlying causes of different numerical convergence, additional large-signal transient simulations were performed, and the evolving electron temperatures inside the device, as produced by the CHT and FKT models, were examined. During the transients, FKT electron temperatures throughout the HEMT structure increased monotonically with drain voltage, and the peak temperature consistently appeared at the drain side of the gate where the lateral electric fields are largest and most of the drain voltage drop occurs. CHT, on the other hand, produced peak electron temperatures located away from the gate and in the degenerately doped drain region. Also, CHT electron temperatures did not increase monotonically throughout the device during the transient drain volt sweep. At early stages of the voltage sweep, CHT produced electron temperatures in the buffer region underneath the gate that dropped nearly 100 K below the ambient temperature. Because these electron dynamics occurred under constant volume conditions, there is no apparent physical mechanism that may have caused this electron cooling. Instead, it may reveal some fundamental differences between the CHT and FKT models, particularly their different treatments of electronic heat flow. These results also suggest that FKT could potentially provide a robust and effective means for simulating GaN HEMT structures for certain applications including RF power electronics.

## ACKNOWLEDGMENTS

This work was supported by AFOSR (Grant No. LRIR 21RYCOR073).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Ashwin Tunga:** Formal analysis (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). **Kexin Li:** Methodology (supporting); Software (equal); Validation (supporting); Writing – review & editing (equal). **Ethan White:** Validation (equal); Writing – review & editing (equal). **Nicholas C. Miller:** Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). **Matt Grupen:** Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). **John D. Albrecht:** Writing – review & editing (equal). **Shaloo Rakheja:** Conceptualization (equal); Funding acquisition (equal); Investigation (equal); Project administration (equal); Resources (equal); Supervision (equal); Validation (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

*2000 International Conference on Simulation Semiconductor Processes and Devices (Cat. No.00TH8502)*(IEEE, Piscataway, NJ, 2000), pp. 23–26.

*Simulation of Semiconductor Processes and Devices 2001*, edited by D. Tsoukalas and C. Tsamis (Springer, Vienna, 2001), pp. 54–57.

*2014 IEEE 26th International Symposium on Power Semiconductor Devices & IC’s (ISPSD)*(IEEE, Piscataway, NJ, 2014), pp. 257–260.

*Proceedings of the NASECODE IV Conference*(Boole Press, 1985), pp. 3–12.

*Simulation of Semiconductor Devices and Processes*(Springer, 1993), pp. 129–132.

*Proceedings of Third International Conference on Electronics, Circuits, and Systems*(IEEE, Piscataway, NJ, 1996), Vol. 2, pp. 752–755.

*2016 74th Annual Device Research Conference (DRC)*(IEEE, Newark, DE, 2016), pp. 1–2, 1.