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.

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 as6 

f(k,r,t)t+vrf(k,r,t)+Fkf(k,r,t)=(f(k,r,t)t)coll,
(1)

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 Gummel11 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 Stratton13 and Blotekjaer14 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 structure26 and full-wave EM.27 FKT has also demonstrated large-signal RF simulations28 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.

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.

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 as20 

Jn=μn[nEC+kTnnnkTnlnγn+λnfntdknTn].
(2)

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. μ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 γn and λn account for the Fermi–Dirac statistics and are given as

γn=nNCexp(ηn),
(3)
λn=F1/2(ηn)F1/2(ηn).
(4)

Fj(x) is the Blakemore Fermi–Dirac integral of order-j. The argument ηn is the reduced Fermi energy and is defined as ηn=(FnEC)/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 TngTn+(1g)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

Wnt+Sn=JnECq+dWndt|coll.,
(5)

where Wn is the electron energy density and Sn 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 JnECq 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. Sn is given as

Sn=5rnλn2[kTnqJn+fnhfκ^nTn]
(6)

and

κ^n=k2qnμnTn.
(7)

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

κn=32k2λnqnμnTn.
(8)
TABLE I.

Parameters in CHT.

ParameterDefault valueMeaning
rn 0.6 Parameter in energy flux 
fnhf Parameter in energy flux 
g Parameter to select carrier temperature [Tn → gTn + (1 − g TL)] 
fntd Thermal diffusion coefficient 
nmin 103 cm−3 Parameter to improve numeric stability 
ParameterDefault valueMeaning
rn 0.6 Parameter in energy flux 
fnhf Parameter in energy flux 
g Parameter to select carrier temperature [Tn → gTn + (1 − g TL)] 
fntd Thermal diffusion coefficient 
nmin 103 cm−3 Parameter to improve numeric stability 

The collision term in the energy balance is given as

dWndtcoll.=ζnWnWn0τe.
(9)

Wn0 is the equilibrium electron energy density, and τe is the energy relaxation time for electrons. The parameter ζn is adjusted to improve the stability of the numerical algorithm and is given as

ζn=1+nminn[n0nmin]max[0,(TLTn)/100K],
(10)

with nmin and n0 are adjustable density parameters to aid convergence. Since TnTL for our device, n0 does not impact ζn. The default value of nmin is set to 103cm3, which leads to ζn1 for all test cases reported here. While τ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 τe’s energy dependence. In this study, we vary τe in the fs to ps range to analyze its impact on IV curves and electron heating. Additionally, we elucidate that τe’s impact on IV 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,1EC,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,

Jn=μnijkTn(niB(uiuj)njB(ujui)),
(11)

where i and j are the major mesh points, μn is the electron mobility (assumed constant between the mesh points), ij is the length of the mesh edge, B(x)=x/(ex1) is the Bernoulli function, uiuj=qEij/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.

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

Jn=μn(2mn)3/23π23[(kTn)3/2F3/2(ηn)EC+(kTn)5/2F3/2(ηn)],
(12)

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

Snkin=μnq(2mn)3/23π23[(kTn)5/2F5/2(ηn)EC+(kTn)7/2F7/2(ηn)].
(13)

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

Jn=μnL(2mn)3/23π23[B(ξn)(kT0)5/2F3/2(ηn,0)B(ξn)(kT1)5/2F3/2(ηn,1)]=q(J01J10),
(14)
ξn=Δ(EC)Δ(FnEC)[lnF3/2(ηn,1)lnF3/2(ηn,0)],
(15)

where the Bernoulli function B(ξn)=ξn/[exp(ξn)1] and Δ(x)x1x0. The same procedure is applied to (13) to obtain the Scharfetter–Gummel discretization of the electron kinetic energy flux,

Snkin=μnqL(2mn)3/23π23[B(ξE)(kT0)7/2F5/2(ηn,0)B(ξE)(kT1)7/2F5/2(ηn,1)]=S01kinS10kin,
(16)
ξE=Δ(EC)Δ(FnEC)[lnF5/2(ηn,1)lnF5/2(ηn,0)].
(17)

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 (FnEC)0 and temperature Tn,0 to a final mesh point with Fermi level (FnEC)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 

H=Wn(Tf)Wn(Ti)(FnEC)[n(Tf)n(Ti)].
(18)

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 J01 and J10 from (14), respectively. Therefore, the electronic heat flux flowing from point 0 to point 1 is given by

S01heat=S01kin(Tn,0)S01kin(Tn,1)(FnEC)0[J01(Tn,0)J01(Tn,1)],
(19)

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

S10heat=S10kin(Tn,1)S10kin(Tn,0)(FnEC)1[J10(Tn,1)J10(Tn,0)],
(20)

and the net electronic heat flux along the edge is Snheat=S01heatS10heat. 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).

The GaN HEMT structure simulated in FKT and CHT is illustrated in Fig. 1. The n+ GaN source/drain contacts are doped with 1019cm3 donors and complete dopant ionization is assumed. The GaN region is semi-insulating containing 9×1015cm3 fully ionized acceptor atoms compensated with 1016cm3 donor-like traps located at an energy level 0.6 eV below EC and capture cross section of 1019m2. 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.

FIG. 1.

Cross section of a GaN device simulated in FKT and CHT. The 2DEG exists in GaN at the interface between AlGaN (barrier) and GaN (channel). Figure is not drawn to scale.

FIG. 1.

Cross section of a GaN device simulated in FKT and CHT. The 2DEG exists in GaN at the interface between AlGaN (barrier) and GaN (channel). Figure is not drawn to scale.

Close modal
FIG. 2.

Mesh structure of the GaN HEMT device, shown in Fig. 1, used for the simulations in both CHT and FKT. The mesh is refined near the 2DEG in the channel, and near the material interfaces, to obtain accurate results.

FIG. 2.

Mesh structure of the GaN HEMT device, shown in Fig. 1, used for the simulations in both CHT and FKT. The mesh is refined near the 2DEG in the channel, and near the material interfaces, to obtain accurate results.

Close modal

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.

FIG. 3.

Comparison of the equilibrium band diagram produced by CHT and FKT along the z-direction (vertical) at x=0.9μm. EC denotes the conduction band edge and EF refers to the Fermi level at thermal equilibrium.

FIG. 3.

Comparison of the equilibrium band diagram produced by CHT and FKT along the z-direction (vertical) at x=0.9μm. EC denotes the conduction band edge and EF refers to the Fermi level at thermal equilibrium.

Close modal
FIG. 4.

Comparison of the equilibrium electron density produced by CHT and FKT along the z-direction (vertical) at x=0.9μm.

FIG. 4.

Comparison of the equilibrium electron density produced by CHT and FKT along the z-direction (vertical) at x=0.9μm.

Close modal
FIG. 5.

Comparison of the equilibrium electron density produced by CHT and FKT along the x-direction (channel) at 0.6 nm below the AlGaN/GaN interface.

FIG. 5.

Comparison of the equilibrium electron density produced by CHT and FKT along the x-direction (channel) at 0.6 nm below the AlGaN/GaN interface.

Close modal
TABLE II.

Material parameters that are common in both CHT and FKT.

MaterialDielectric constantElectron affinity (eV)Bandgap (eV)Electron mobility (cm2/V s)Effective mass (m0)a
AlGaN 8.9 3.55 4.14 … … 
GaN 8.9 4.1 3.457 440 0.2 
SiC 9.66 … … … … 
MaterialDielectric constantElectron affinity (eV)Bandgap (eV)Electron mobility (cm2/V s)Effective mass (m0)a
AlGaN 8.9 3.55 4.14 … … 
GaN 8.9 4.1 3.457 440 0.2 
SiC 9.66 … … … … 
a

m0 is the free electron mass

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.

FIG. 6.

Drain current ID vs drain voltage VDS for different gate voltages VGS simulated by CHT and FKT with no electron heating, i.e., drift-diffusion.

FIG. 6.

Drain current ID vs drain voltage VDS for different gate voltages VGS simulated by CHT and FKT with no electron heating, i.e., drift-diffusion.

Close modal

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 (τ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 IDVDS 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.

FIG. 7.

Drain current ID vs drain voltage VDS for different gate voltages VGS simulated by CHT and FKT including hot electrons, τe=0.2 ps.

FIG. 7.

Drain current ID vs drain voltage VDS for different gate voltages VGS simulated by CHT and FKT including hot electrons, τe=0.2 ps.

Close modal
FIG. 8.

Comparison of electron temperatures produced by CHT and FKT along the x-direction (channel) at 6 nm below the AlGaN/GaN interface for VGS=0 V, VDS=10 V, and τe=0.2 ps for both simulators.

FIG. 8.

Comparison of electron temperatures produced by CHT and FKT along the x-direction (channel) at 6 nm below the AlGaN/GaN interface for VGS=0 V, VDS=10 V, and τe=0.2 ps for both simulators.

Close modal

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 τ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 nth 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=4 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 <1010.

FIG. 9.

Calculated Newton solver residuals of CHT (left) and FKT (right) hot electron simulations with τe=0.2 ps at VGS=4 V. The residuals are functions of the Newton iteration and the quiescent drain voltage.

FIG. 9.

Calculated Newton solver residuals of CHT (left) and FKT (right) hot electron simulations with τe=0.2 ps at VGS=4 V. The residuals are functions of the Newton iteration and the quiescent drain voltage.

Close modal

A quantitative investigation of the convergence properties of the transport solvers is next presented. In particular, we analyze the rate of convergence (ROC), α, of the numerical solvers. The ROC is calculated by37 

α(VGS,VDS)log|en+1(VGS,VDS)/en(VGS,VDS)|log|en(VGS,VDS)/en1(VGS,VDS)|.
(21)

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 (α=2) for all the drain biases. It should be noted that there exists some FKT ROC values which fall slightly below α=2. These values, on the order of α=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.

FIG. 10.

Estimated ROC of CHT (left) and FKT (right) hot electron simulations with τe=0.2 ps and VGS ranging from 4 to 0 V in 1 V steps. The black dashed line indicates the region of ROC values representative of quadratic convergence.

FIG. 10.

Estimated ROC of CHT (left) and FKT (right) hot electron simulations with τe=0.2 ps and VGS ranging from 4 to 0 V in 1 V steps. The black dashed line indicates the region of ROC values representative of quadratic convergence.

Close modal

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 τ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 τe to 5.0 fs, while keeping the τ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.

FIG. 11.

Comparison of electron temperatures produced by CHT and FKT along the x-direction (channel) at 6 nm below the AlGaN/GaN interface at VGS=0 V. (i) Temperature profile produced by FKT on increasing VDS to 50 V at τe=0.2 ps. (ii) Temperature profile produced by CHT for VDS=10 V and τe=0.2 ps. (iii) Temperature profile produced by CHT for VDS=10 V and τe=5.0 fs.

FIG. 11.

Comparison of electron temperatures produced by CHT and FKT along the x-direction (channel) at 6 nm below the AlGaN/GaN interface at VGS=0 V. (i) Temperature profile produced by FKT on increasing VDS to 50 V at τe=0.2 ps. (ii) Temperature profile produced by CHT for VDS=10 V and τe=0.2 ps. (iii) Temperature profile produced by CHT for VDS=10 V and τe=5.0 fs.

Close modal
FIG. 12.

Estimated ROC of CHT (left) hot electron simulations with τe=0.005 ps and σPZ=9×1012cm2. Estimated ROC of FKT (right) hot electron simulations with τe=0.2 ps, σPZ=9×1012cm2, and drain voltage up to 50 V. The symbols correspond to VGS, as indicated in the legends.

FIG. 12.

Estimated ROC of CHT (left) hot electron simulations with τe=0.005 ps and σPZ=9×1012cm2. Estimated ROC of FKT (right) hot electron simulations with τe=0.2 ps, σPZ=9×1012cm2, and drain voltage up to 50 V. The symbols correspond to VGS, as indicated in the legends.

Close modal

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 τe far below 0.2 ps, it might seem natural to increase τe far above 0.2 ps in FKT to produce much higher Tn. However, when τ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 (σ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 σPZ, with τ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 σPZ, while CHT shows a super-linear ROC that decreases as σPZ is increased. Furthermore, the CHT simulations failed to converge for all (VGS, VDS), when σPZ15×1012cm2, which is indicated as “N/A” in Table III. Here, mean{α(VGS,VDS)} is the average ROC of all (VGS, VDS) at which the simulations converged, where VGS[4,0] V and VDS[0,10] V, as shown in Figs. 10 and 12.

TABLE III.

Summary of CHT device and FKT ROC values. “N/A” refers to simulations that failed to converge.

mean{α(VGS, VDS)}
σPZ (×1012 cm−2)CHTFKT
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{α(VGS, VDS)}
σPZ (×1012 cm−2)CHTFKT
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 

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.

FIG. 13.

Electron temperatures Tn across the transistor for VGS=0 when VDS is increased from 0 to 1 V in 1 ns as computed with the CHT model (left) and the FKT model (right). The AlGaN barrier/GaN channel interface is located at z = 0, and the position of the gate in the x-direction is represented by the dotted lines. Multimedia view: https://doi.org/10.1063/5.0118104.1, https://doi.org/10.1063/5.0118104.2

FIG. 13.

Electron temperatures Tn across the transistor for VGS=0 when VDS is increased from 0 to 1 V in 1 ns as computed with the CHT model (left) and the FKT model (right). The AlGaN barrier/GaN channel interface is located at z = 0, and the position of the gate in the x-direction is represented by the dotted lines. Multimedia view: https://doi.org/10.1063/5.0118104.1, https://doi.org/10.1063/5.0118104.2

Close modal
FIG. 14.

Electron temperatures Tn across the transistor computed by with the CHT model (left) and the FKT model (right) over the time span 0<t<0.15 ns as the drain voltage VDS changed at a rate of 1 V/ns. The AlGaN barrier/GaN channel interface is located at z = 0, and the position of the gate in the x-direction is represented by the dotted lines. Multimedia view: https://doi.org/10.1063/5.0118104.3, https://doi.org/10.1063/5.0118104.4

FIG. 14.

Electron temperatures Tn across the transistor computed by with the CHT model (left) and the FKT model (right) over the time span 0<t<0.15 ns as the drain voltage VDS changed at a rate of 1 V/ns. The AlGaN barrier/GaN channel interface is located at z = 0, and the position of the gate in the x-direction is represented by the dotted lines. Multimedia view: https://doi.org/10.1063/5.0118104.3, https://doi.org/10.1063/5.0118104.4

Close modal

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.

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

The authors have no conflicts to disclose.

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).

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

1.
K.
Hoo Teo
,
Y.
Zhang
,
N.
Chowdhury
,
S.
Rakheja
,
R.
Ma
,
Q.
Xie
,
E.
Yagyu
,
K.
Yamanaka
,
K.
Li
, and
T.
Palacios
, “
Emerging GaN technologies for power, RF, digital, and quantum computing applications: Recent advances and prospects
,”
J. Appl. Phys.
130
,
160902
(
2021
).
2.
B. J.
Baliga
, “
Gallium nitride devices for power electronic applications
,”
Semicond. Sci. Technol.
28
,
074011
(
2013
).
3.
M.
Stettler
,
M.
Alam
, and
M.
Lundstrom
, “
A critical examination of the assumptions underlying macroscopic transport equations for silicon devices
,”
IEEE Trans. Electron Devices
40
,
733
740
(
1993
).
4.
J.
Bude
, “MOSFET modeling into the ballistic regime,” in 2000 International Conference on Simulation Semiconductor Processes and Devices (Cat. No.00TH8502) (IEEE, Piscataway, NJ, 2000), pp. 23–26.
5.
T.
Grasser
,
H.
Kosina
, and
S.
Selberherr
, “Investigation of spurious velocity overshoot using Monte Carlo data,” in Simulation of Semiconductor Processes and Devices 2001, edited by D. Tsoukalas and C. Tsamis (Springer, Vienna, 2001), pp. 54–57.
6.
K.
Hess
,
Advanced Theory of Semiconductor Devices
(
Prentice Hall
,
1988
).
7.
M. V.
Fischetti
and
S. E.
Laux
, “
Monte Carlo analysis of electron transport in small semiconductor devices including band-structure and space-charge effects
,”
Phys. Rev. B
38
,
9721
9745
(
1988
).
8.
M.
Saraniti
and
S.
Goodnick
, “
Hybrid fullband cellular automaton/Monte Carlo approach for fast simulation of charge transport in semiconductors
,”
IEEE Trans. Electron Devices
47
,
1909
1916
(
2000
).
9.
S.
Selberherr
,
Analysis and Simulation of Semiconductor Devices
(
Springer Science & Business Media
,
1984
).
10.
T.
Grasser
,
T.-W.
Tang
,
H.
Kosina
, and
S.
Selberherr
, “
A review of hydrodynamic and energy-transport models for semiconductor device simulation
,”
Proc. IEEE
91
,
251
274
(
2003
).
11.
D.
Scharfetter
and
H.
Gummel
, “
Large-signal analysis of a silicon read diode oscillator
,”
IEEE Trans. Electron Devices
16
,
64
77
(
1969
).
12.
O.
Penzin
,
L.
Smith
,
A.
Erlebach
,
M.
Choi
, and
K.-H.
Lee
, “
Kinetic velocity model to account for ballistic effects in the drift-diffusion transport approach
,”
IEEE Trans. Electron Devices
64
,
4599
4606
(
2017
).
13.
R.
Stratton
, “
Diffusion of hot and cold electrons in semiconductor barriers
,”
Phys. Rev.
126
,
2002
2014
(
1962
).
14.
K.
Blotekjaer
, “
Transport equations for electrons in two-valley semiconductors
,”
IEEE Trans. Electron Devices
17
,
38
47
(
1970
).
15.
S.
Strauss
,
A.
Erlebach
,
T.
Cilento
,
D.
Marcon
,
S.
Stoffels
, and
B.
Bakeroot
, “TCAD methodology for simulation of GaN-HEMT power devices,” in 2014 IEEE 26th International Symposium on Power Semiconductor Devices & IC’s (ISPSD) (IEEE, Piscataway, NJ, 2014), pp. 257–260.
16.
R.
Kraus
and
A.
Castellazzi
, “
A physics-based compact model of SiC power MOSFETs
,”
IEEE Trans. Power Electron.
31
,
5863
5870
(
2016
).
17.
J.
Ajayan
,
T.
Ravichandran
,
P.
Mohankumar
,
P.
Prajoon
,
J.
Charles Pravin
, and
D.
Nirmal
, “
Investigation of DC-RF and breakdown behaviour in Lg = 20 nm novel asymmetric GaAs MHEMTs for future submillimetre wave applications
,”
Int. J. Electron. Commun.
84
,
387
393
(
2018
).
18.
P.
Zheng
,
D.
Connelly
,
F.
Ding
, and
T.-J. K.
Liu
, “
FinFET evolution toward stacked-nanowire FET for CMOS technology scaling
,”
IEEE Trans. Electron Devices
62
,
3945
3950
(
2015
).
19.
Y.-C.
Wu
and
Y.-R.
Jhan
,
3D TCAD Simulation for CMOS Nanoeletronic Devices
(
Springer
,
2017
).
20.
Synopsys
,
Sentaurus Device User Guide
(
Synopsys, Inc.
,
Mountain View, CA
,
2020
).
21.
G.
Baccarani
, “HFIELDS: A highly flexible 2-D semiconductor-device analysis program,” in Proceedings of the NASECODE IV Conference (Boole Press, 1985), pp. 3–12.
22.
J.
Litsios
,
S.
Müller
, and
W.
Fichtner
, “Mixed-mode multi-dimensional device and circuit simulation,” in Simulation of Semiconductor Devices and Processes (Springer, 1993), pp. 129–132.
23.
G.
Baccarani
,
M.
Rudan
,
M.
Lorenzini
,
W.
Fichtner
,
J.
Litsios
,
A.
Schenk
,
P.
van Staa
,
L.
Kaeser
,
A.
Kampmann
,
A.
Marmiroli
,
C.
Sala
, and
E.
Ravanelli
, “Device simulation for smart integrated systems (DESSIS),” in Proceedings of Third International Conference on Electronics, Circuits, and Systems (IEEE, Piscataway, NJ, 1996), Vol. 2, pp. 752–755.
24.
C. K.
Sarkar
,
Technology Computer Aided Design: Simulation for VLSI MOSFET
(
CRC Press
,
2018
).
25.
M.
Grupen
, “
An alternative treatment of heat flow for charge transport in semiconductor devices
,”
J. Appl. Phys.
106
,
123702
(
2009
).
26.
M.
Grupen
, “
Energy transport model with full band structure for GaAs electronic devices
,”
J. Comput. Electron.
10
,
271
290
(
2011
).
27.
M.
Grupen
, “
Three-dimensional full-wave electromagnetics and nonlinear hot electron transport with electronic band structure for high-speed semiconductor device simulation
,”
IEEE Trans. Microwave Theory Tech.
62
,
2868
2877
(
2014
).
28.
N. C.
Miller
,
J. D.
Albrecht
, and
M.
Grupen
, “Large-signal RF GaN HEMT simulation using Fermi kinetics transport,” in 2016 74th Annual Device Research Conference (DRC) (IEEE, Newark, DE, 2016), pp. 1–2, 1.
29.
M.
Grupen
, “
Reproducing GaN HEMT kink effect by simulating field-enhanced barrier defect ionization
,”
IEEE Trans. Electron Devices
66
,
3777
3783
(
2019
).
30.
N. C.
Miller
,
J. D.
Albrecht
, and
M.
Grupen
, “
Delaunay-Voronoi surface integration: A full-wave electromagnetics discretization for electronic device simulation
,”
Int. J. Numer. Model.
29
,
817
830
(
2016
).
31.
N. C.
Miller
,
Large-Signal RF Simulation and Characterization of Electronic Devices Using Fermi Kinetics Transport
(
Michigan State University
,
2017
).
32.
N. C.
Miller
,
M.
Grupen
,
K.
Beckwith
,
D.
Smithe
, and
J. D.
Albrecht
, “
Computational study of Fermi kinetics transport applied to large-signal RF device simulations
,”
J. Comput. Electron.
17
,
1658
1675
(
2018
). 2.
33.
R.
Stratton
, “
Diffusion of hot and cold electrons in semiconductor barriers
,”
Phys. Rev.
126
,
2002
(
1962
).
34.
C.
Kittel
and
H.
Kroemer
,
Thermal Physics
, 2nd ed. (
Freeman
,
1980
).
35.
H.
Ye
,
G. W.
Wicks
, and
P. M.
Fauchet
, “
Hot electron relaxation time in GaN
,”
Appl. Phys. Lett.
74
,
711
713
(
1999
).
36.
J.-Z.
Zhang
, “
Energy relaxation of hot electrons in Si-doped GaN
,”
J. Appl. Phys.
115
,
203704
(
2014
).
37.
J. R.
Senning
,
Computing and Estimating the Rate of Convergence
(
Gordon College
,
Wenham
,
2007
).