## Abstract

A model of self-heating is incorporated into a cellular Monte Carlo (CMC) particle device simulator. This is done through the solution of an energy balance equation (EBE) for phonons, which self-consistently couples charge and heat transport in the simulation. First, several tests are performed to verify the applicability and accuracy of the proposed nonlinear iterative solver in the presence of convective boundary conditions, as compared to a finite element analysis (FEA) solver as well as using the Kirchhoff Transformation. Finally, a fully coupled electro-thermal characterization of a GaN/AlGaN high electron mobility transistor (HEMT) is performed, and the effects of nonideal interfaces and boundary conditions are studied.

## 1 Introduction

It is well known that electronic devices suffer from reliability issues and performance degradation due to self-heating effects [1,2]. Self-heating and thermal management are therefore still the subjects of active research, and thermal effects are often not included in the context of Monte Carlo (MC) simulations at all.

Monte Carlo techniques have been widely used in the study and characterization of semiconductor materials and devices in a number of seminal works [36]. While being less computationally expensive, drift-diffusion and hydrodynamic simulation approaches are based on only low order moments of the Boltzmann transport equation (BTE), and they become inaccurate as devices are pushed smaller in size and closer to their material limits. In comparison, MC methods give a full stochastic solution to the BTE while making as few assumptions as possible; this comes at the expense of being more costly in terms of either simulation time or hardware requirements, or both. Here, we adopt the cellular Monte Carlo (CMC) approach [7], which greatly reduces the time of simulation at the expense of increased memory requirements in comparison with traditional ensemble Monte Carlo methods.

An excellent and comprehensive thermodynamic treatment of self-heating and heat transport in devices is given by Wachutka [8]. Fushinobu et al. [9] developed an approach to include self-heating by assuming an energy decay path and using a relaxation time approximation that has also been used elsewhere [10,11]. On the other hand, Pilgrim et al. modeled heat generation through a net phonon emission approach and solved a thermal resistance matrix to obtain temperature maps [12]. While these approaches incorporate temperature-dependent scattering rates, the temperature dependence of material properties, like thermal conductivity, is often neglected.

A rather comprehensive model of electro-thermal GaN device simulation is that of Hao et al. [13,14] which couples electron and phonon MC simulations along the 2-dimensional electron gas (2DEG) channel where heat generation is most significant while using Fourier's Law in the rest of the domain. In Refs. [13] and [14], they extract parameters such as electron density, energy, and drift velocity from the electron MC to compute the energy exchanged between electrons and phonons. In contrast, here we do not carry out full phonon MC simulations solving the phonon BTE for heat transport but instead solve a steady-state energy balance equation (EBE) based on Fourier's law. However, in this work heat generation is computed directly from the principles of energy conservation: as electron–phonon scattering events occur, the energy lost by electrons is recorded and used to compute a heat generation rate in each computational cell in the simulation domain. Also, we compute the energy transfer rate from electrons to both optical and acoustic phonons, respectively, rather than assuming an energy decay path from electrons to optical phonons to acoustic phonons. We feel that this novel way of computing the heat generation is more physically accurate than methods based on a relaxation-time approximation. In addition, we use a full electronic bandstructure and phonon dispersion, including all optical phonon modes without assuming them to be dispersionless, i.e., we allow the optical modes to contribute to the thermal conductivity based on the material-specific phonon dispersion. However, Fourier's law becomes inaccurate at the nanoscale where ballistic heat transfer effects become significant, and hence this model is unable to include such effects.

Bonani and Ghione [15] took the temperature-dependent thermal conductivity into account in the context of device simulation by using the Kirchhoff transformation [16], which has also been used in the study of larger scale thermal spreading problems in semiconductor packages [17]. However, convective boundary conditions, which physically represent an imperfect heat sink, become nonlinear under the Kirchhoff transformation. In addition, these types of conditions can also occur at imperfect material interfaces, which can introduce a “thermal resistance” to heat conduction across the interface. Hence, material interfaces represent a potential bottleneck for heat dissipation, while at the same time becoming more significant in high electron mobility transistors (HEMTs) where often a GaN/AlN heterojunction is used to create a 2DEG conduction channel for electron transport.

In this work, we expand the capabilities of the CMC framework to include a coupled flux-based electro-thermal device solver using a two-step iterative process. This allows us to retain the speed advantages of the CMC, while obtaining a physically accurate picture of steady-state heat generation within a device. Energy-dependent carrier scattering rates for materials have been extensively studied and are well accepted. Here, heat generation is computed directly from the principles of energy conservation: as electron–phonon scattering events occur, the energy lost by electrons is used to compute a heat generation rate in each computational cell in the simulation.

This heat generation rate is used as the forcing function in the energy balance equations for the optical and acoustic phonon modes, respectively; meanwhile, the relaxation time approximation is retained for the energy exchange between the optical and acoustic phonon modes. The resulting temperature map for each phonon mode is then used to update the temperature-dependent thermal conductivity, and this inner iterative loop is repeated until convergence is reached everywhere in the simulation. The scattering rates in each cell are subsequently updated according to the new temperature (as scattering increases with temperature), a new heat generation rate is obtained, and the process is repeated until the heat generation rate stabilizes and the temperature maps reach convergence in the outer iterative loop.

The final result is a temperature map for both acoustic and optical modes and the full electro-thermal operating characteristics, e.g., an IV curve accounting for the effects of self-heating.

## 2 Theory

The details of the CMC approach and its advantages over traditional ensemble Monte Carlo methods have been detailed previously [7], as well as the extension of the CMC to subsequent applications, such as HEMT and radio frequency (RF) simulation. Therefore, this work focuses solely on the development of a thermal solver that we have coupled to the electrical CMC and validated through coupled nonlinear electro-thermal simulations.

The principal novelties in this work are the inclusion of temperature-dependent thermal conductivities and convective, or Robin, boundary conditions, as well as a novel approach to calculating the position-dependent heat generation rates within the device under study.

### 2.1 Energy Balance Equation.

The particle-flux-based electro-thermal simulator developed here solves an EBE for phonons, derived directly from the phonon BTE, in lieu of the heat transport equation. However, it will be shown that the final EBE has the same general form as the heat transport equation. The advantages are that the EBE is better suited to self-consistent coupling with the electron dynamics and allows for greater accuracy by giving a separate solution for each phonon mode (or groups of modes, i.e., all acoustic or optical). The EBE for each phonon mode μ is written
$∂Wμδt=−∇·Fμ+∂Wμδt|e−p+∂Wμδt|p−p$
(1)

where $Wμ(r,t)=(1/Ω)ΣkEμ(k)fμ(r,k,t)$ is the ensemble energy density within the volume Ω of the reciprocal space, $Fμ(r,t)=(1/Ω)Σkv(k)Eμ(k)fμ(r,k,t)$ is the energy flux, and the two partial derivatives on the right-hand side constitute the generation terms due to electron–phonon (ep) and phonon–phonon (pp) interactions, respectively.

Under steady-state conditions, the time derivative on the left-hand side is zero, and the sum of electron–phonon and phonon–phonon contributions gives the total heat generation rate, $Pμ$
$∇·(κμ(r,T)∇T)=−(∂Wμδt|e−p+∂Wμδt|p−p)=−Pμ$
(2)

where Fourier's law has been used to approximate the energy flux, $Fμ(r)=−κμ(r,T)∇T$. Note that this means that ballistic heat transfer is not explicity included in our model, as Fourier's law is diffusive in nature implying the presence of scattering.

It should be noted at this point that Eq. (2) can be directly solved with a nonlinear finite element analysis (FEA) approach. Although in principle a MC device simulator could be based on finite element methods, particle tracking could become more complicated when using nonrectangular grids. Hence, our CMC (and many academic device simulation codes) uses finite difference methods, where solving such an equation is nontrivial.

Furthermore, a solver for a linear elliptical partial differential equation (PDE), the Poisson equation, is already required in any device simulator. Hence, manipulating Eq. (2) into the form of a linear elliptical PDE allows for a rather straightforward implementation in any device simulator using finite difference methods.

The main issue in performing this manipulation is the dependence of $κμ,C(r,T)$ on both position and temperature. For the former, we assume that the thermal conductivity is constant within any particular cell, C, but allowed to vary when moving from one cell to another. This is analogous to the approach used for the dielectric constant in Poisson's equation. Hence, we will express this restricted position dependency with the cell index C rather than via a full functional dependence on the position vector $r$. For the latter case of the temperature dependence, we've taken two separate approaches to perform this manipulation: (1) using the Kirchhoff transformation and (2) an iterative approach which overcomes many of the restrictions imposed by the Kirchhoff transformation.

### 2.2 Kirchhoff Transformation.

The Kirchhoff transformation is defined as [1517]
$Θμ,C(T)=Tref+1κμ,C(Tref)∫TrefTκμ,C(τ)∂τ$
(3)
where $Θ$ denotes a new “apparent” temperature, which implicitly includes the temperature dependence of the thermal conductivity. Tref is a chosen reference temperature, often 300 K, at which the thermal conductivity value, $κμ,C(Tref)$ is independently known. We choose the temperature-dependent thermal conductivity to be given by a power law fit to experimental data as in Ref. [18]
$κμ(T)=κref(TTref)α$
(4)
Using the apparent temperature, $Θ$, we can now write the EBE as a linear elliptical PDE
$∇2Θμ,C=−Pμ(r)κμ,C(Tref)$
(5)

However, a few significant restrictions are imposed by the Kirchhoff Transformation. First, Eq. (5) requires that both the temperature and its derivative in the normal direction, i.e., the heat flux, be continuous across a boundary [15]. This can be overcome by using the same value for the exponent α in the power law fit across the boundary, but this introduces a further approximation if the materials have a differing temperature dependence. Moreover, while the approximation is reasonable for most semiconductor materials, metals and oxides display vastly different temperature-dependent behavior which is relevant at their respective semiconductor interfaces.

Second, functions representing Dirichlet and simple Neumann boundary conditions, corresponding to a constant temperature or heat flux, respectively, are linear under the Kirchhoff transformation. In contrast, convective boundary conditions become nonlinear under the Kirchhoff transformation [16,17]. Convective boundary conditions describe thermal contacts that have an inherent thermal resistance, Rth, or a finite heat transfer coefficient, H, where H =1/Rth. Clearly, as $Rth→0, H→∞$ and the perfect thermal contact of a Dirichlet boundary condition is recovered. In addition, convective conditions can also occur at material interfaces, which often introduce a resistance to heat passing from one material to the other. Convective boundary conditions are expressed as
$−κ∇T=H(TL−TR)$
(6)

which relates the heat flux approaching the interface, $−κ∇T$, with the heat transport coefficient and the temperature jump across the interface, $(TL−TR)$. Since the Kirchhoff transformation is nonlinear for any nonlinear functional form of the thermal conductivity, it is generally not true that $H(ΘL−ΘR)=H(TL−TR)$. Bagnall et al. [17] overcame this nonlinearity in large-scale thermal spreading problems by choosing the reference temperature and thermal conductivity in the Kirchhoff transformation such that the convective boundary condition is roughly linear, i.e., $H(ΘL−ΘR)≈H(TL−TR)$. This approach requires some a priori knowledge of the temperature distribution expected at the convective boundary in order to choose $Tref$ around which to linearize the transformation so as to reduce the error due to the nonlinearity. However, this poses problems for coupled electro-thermal device simulations where the heat generation is not necessarily known ahead of time. In addition, Bagnall demonstrated that in cases where the temperature distribution along the convective boundary is less uniform this approach becomes less accurate. This is simply because a single temperature is chosen around which to linearize the Kirchhoff transformation, and temperatures far away from the chosen value become more inaccurate, again due to the nonlinearity.

### 2.3 Iterative Method.

To overcome the discussed restrictions related to the Kirchhoff transformation, a classical iterative approach to solving the EBE has been implemented. This takes advantage of the fact that the thermal conductivity is allowed to vary from one cell to another to incorporate the temperature-dependent values, and a sequence of constant-conductivity problems is solved, where the thermal conductivity in each cell is updated with each new solution according to the local temperature. In other words, we solve a system of the form
$∇2Tμ,Ci=−Pμ(r)κμ,C(Ti−1)$
(7)

where the new temperature $Ti$, for iteration i, on the left-hand side is obtained using the temperature-dependent thermal conductivity in each cell obtained from the previous temperature solution from iteration i − 1. This iterative approach uses the same power law fit previously discussed for the thermal conductivity; however, it relaxes the requirement that the exponent α in the temperature dependence be the same for different materials.

In addition, while the convective boundary condition in Eq. (6) is linear in and of itself, it becomes nonlinear under the Kirchhoff transformation as mentioned previously. However, convective boundary conditions are accounted for in a straightforward manner under the iterative approach, i.e., in the same manner Dirichlet and Neumann conditions are handled. Hence, this approach allows us to accurately simulate heat transport across interfaces possessing nonideal effects, including realistic convective heat sinks and thermally resistive material interfaces.

This is significant for devices such as a high electron mobility transistor, which is studied in Sec. 3.4, where material interfaces are used to create a 2DEG conduction channel, as well as in layouts attempting to use a high thermal conductivity material to transport heat away from the active device, where the interface can have an associated thermal resistance.

## 3 Simulation Results

In order to validate our iterative approach, we have reproduced simulation layouts from Refs. [17] and [19] in one-dimensional (1D) and two-dimensional (2D) cases. In addition, a validation of the simulation of interface resistances was performed on a ten-layer structure similar to that of Hickson et al. [20], demonstrating not only its accuracy but its robustness as well. Finally, we apply this approach to a GaN/AlGaN HEMT to obtain a validation in fully coupled electro-thermal simulations.

### 3.1 One-Dimensional and Two-Dimensional.

To verify the proposed iterative approach, we have duplicated the simulation tests performed in Refs. [17] and [19]. Since the results in Ref. [17] were used as a validation of the Kirchhoff transformation approach itself, here we have only shown the Kirchhoff transformation results when they differ significantly from the iterative solver due to the nonlinearity. Instead, we use as a reference the results obtained using the full nonlinear FEA solver in the commercial matlab solver [21] that is able to solve the nonlinear problem exactly.

In the simplest case of a 1D simulation, we use the material parameters of a 100-μm-thick piece of silicon. The thermal conductivity of silicon at 300 K is taken as 150 W/m K, and the exponent in the temperature power law of Eq. (4) is taken as $α=−1.3$ [18]. In other words, the temperature dependent thermal conductivity is taken as $κμ(T)=150(T/300)(−1.3)$.

The boundary conditions used, also indicated in the schematic inset in Fig. 1, are a constant inward heat flux all along the top boundary of 5 × 106 W/m2, while a convective boundary is set along the bottom of the material characterized by a heat transfer coefficient H =105 W/m2 K and an ambient temperature $T∞=350$ K. As seen in Fig. 1, our iterative method agrees with the full nonlinear solution using both a constant and temperature-dependent thermal conductivity.

Fig. 1
Fig. 1
Close modal

Next, we consider a more complicated 2D heat transfer problem. Here, rather than a constant heat flux across the entire top boundary of the material, a single 2-nm-wide heat source of 1 × 109 W/m2 is centered on the top surface at x =250 μm. All other simulation parameters are the same as used previously, i.e., the layout size and the material properties. In Fig. 2, we show the result of the iterative solver compared with that of the full nonlinear FEA. The inset shows the area directly around the central heat source along the top of the material where the temperature gradient is largest. Again, the iterative approach shows excellent agreement with the results of the nonlinear FEA. Note that this is also seen to be true when using the Kirchhoff transformation in this particular layout [17], and hence this serves as a benchmark with the advantage being that the iterative approach does not require any a priori knowledge of the temperature distribution near the convective boundary.

Fig. 2
Fig. 2
Close modal

However, as stated previously, Bagnall et al. [17] found that the approach of effectively linearizing the Kirchhoff transformation around a suitably chosen temperature becomes inaccurate at smaller thicknesses. For this reason, the iterative approach was tested on the same 500-μm-long silicon slab with its thickness reduced to 5 μm. In Fig. 3, we show a comparison between the full nonlinear FEA, the linearized Kirchhoff transformation, and the iterative approach solutions, respectively. At this scale, the iterative approach is seen to be more accurate than the linearized Kirchhoff transformation, while still requiring no a priori knowledge of the temperature distribution near the convective boundary.

Fig. 3
Fig. 3
Close modal

### 3.2 Interface Resistance.

To implement the material interface condition, we follow the work of Hickson et al. [20]. However, here we are interested in solving the steady-state problem rather the time-dependent problem studied in Ref. [20]. In this work, we use a first-order expression for the interface coefficients as the second order expressions are much more cumbersome [20], and, since we generally work with nonhomogeneous grids, we expect first-order accuracy anyway. Considering only the 1D case on the x-axis, as it is a straightforward generalization in the other directions, the second derivatives on each side of the interface are written
$∂Ti−1∂x2≈ΔxwTL−(ΔxL+Δxw)Ti−1ΔxLΔxw(ΔxL+Δxw)$
(8)
$∂Ti+1∂x2≈ΔxeTR−(ΔxR+Δxe)Ti+1ΔxLΔxe(ΔxL+Δxe)$
(9)
where i − 1 denotes the nearest grid point to the left of the interface, i +1 the nearest grid point to the right of the interface, ΔxL the distance between point i − 1 and the interface, ΔxR the distance between point i +1 and the interface, Δxw the distance between points i − 1 and i − 2 to the left, Δxe the distance between points i +1 and i +2 to the right. In Eqs. (8) and (9), we have four unknowns: $Ti−1, Ti+1, TL$, and $TR$. The temperatures at the grid points everywhere in the simulation, such as $Ti−1$ and $Ti+1$, are solved for using numerical methods such as successive over-relaxation, and so we need to find expressions to substitute for $TL$ and $TR$. To find these expressions, we use the fact that the flux must be continuous on each side of the interface along with the boundary condition itself
$κL∂Ti−1∂x=κR∂Ti+1∂x$
(10)
$κL∂Ti−1∂x=H(TR−TL)$
(11)
We can now obtain the required expressions for $TL$ and $TR$ as
$TL=(HΔxR+κR)Ti−1+HΔxLκRTi+1(HΔxR+κR)κL+HΔxLκR$
(12)
$TR=(HΔxL+κL)Ti+1+HΔxRκLTi−1(HΔxL+κL)κR+HΔxRκL$
(13)

To verify the robustness of the iterative approach in the presence of interfaces, a ten-layer structure similar to that in Ref. [20] has been simulated. The boundary conditions are modified from Ref. [20] as we are interested in the steady-state problem, and hence we have imposed a 400 K Dirichlet condition on the left and a 300 K Dirichlet condition on the right of Fig. 4. The alternating layers are labeled and denoted by the light dashed lines. The parameters used are $κA=1$ W/m K, $κB=0.1$ W/m K, and the heat transport coefficient at each interface is Hi = 0.5 W/m2 K.

Fig. 4
Fig. 4
Close modal

Interestingly, we see that the temperature-dependent solution here follows that of the constant conductivity quite closely, which is reasonable as each layer is only 1 nm and the total size of the simulation is 10 nm. The Kirchhoff transformation solution, however, becomes more inaccurate the further it moves away from the Dirichlet conditions due to the nonlinearity.

### 3.3 Multifinger GaN High Electron Mobility Transistor.

One of the most robust tests found in literature for thermal simulation at the package scale, including convective boundary conditions, is that of a GaN RF power amplifier (represented by a constant heat flux) on top of a SiC substrate [17,19]. In Refs. [17] and [19], the RF power amplifier package's thermal characteristics were simulated in a three-dimensional (3D) layout, whereas here we have simplified it to a 2D layout for practical reasons. The generalization to the 3D case is straightforward.

The simulated layout consists of 11 heat sources of 0.5 μm width located at the very top of the material (y =102 μm), each representing the heat generated by a GaN HEMT with an inward flux of 1010 W/m2. The first source is located at x =25 μm, and the sources are separated by 50 μm, with the final one located at x =525 μm. The total size of the entire layout is 1000 by 102 μm. In the y-direction, the lower 100 μm layer is composed of a SiC substrate, while the top 2 μm layer is GaN. Finally, the ambient temperature in the convective boundary condition used here is $T∞=293$ K.

As an initial benchmark, we first performed a constant-conductivity simulation using both the matlab FEA solver and the iterative approach proposed here. This was done to ensure the accuracy of the numerical grid used, and the full convergence of the solutions. The full nonlinear problem was then solved using the same thermal conductivities and temperature dependences found in Refs. [17] and [19] for each material, namely
$κGaN=150(T300)−1.3, κSiC=420(T300)−1.4$
(14)

It should be noted again that when using the Kirchhoff transformation, the exponent must be identical for each material in the simulation, while this is not true neither in the fully nonlinear FEA nor the iterative approach used here. Due to this, if using the Kirchhoff transformation it is necessary to adjust the exponent for GaN to $α=−1.4$ to match the value for SiC, as in Ref. [17]. The result of this simulation is shown in Fig. 5, where the maximum absolute temperature discrepancy between the two approaches is 0.7K.

Fig. 5
Fig. 5
Close modal

### 3.4 GaN/AlGaN High Electron Mobility Transistor Device Simulation.

As mentioned previously, an area of particular interest in electro-thermal simulation is that of device reliability. Self-heating effects are well known to cause a decrease in operating current, as seen in phenomena such as the well-known current collapse in HEMTs. The ability to design devices accounting for their thermal characteristics, rather than simply electrical ones, would help in addressing reliability concerns prior to fabrication.

In our simulations, the CMC framework for charge transport takes into account scattering due to both optical and acoustic phonons (both polar and deformation potential), ionized impurities, and impact ionization. The material properties used for GaN are the same as those in Yamakawa et al. [22]. The 2DEG channel is created by placing a polarization charge sheet at the heterointerface according to Ambacher's formalism [23].

To perform the fully coupled electro-thermal simulation, we first allow the electrical CMC to reach steady-state just as in an isothermal simulation. Then, we average the power density transferred from the charge carriers to the respective phonon modes over a sufficiently large window to reduce statistical noise. That is, we obtain the value of $∂Wμ/δt|e−p$ in Eq. (2) directly from the CMC simulation. The energy exchange between the acoustic and optical phonon modes is treated via the relaxation time approximation
$∂Wμδt|p−p=Ci(Ti−Tjτi,j)$
(15)

where Ci is the volumetric heat capacity of the phonon mode, $T(i,j)$ the temperature of the respective mode, and $τ(i,j)$ the decay relaxation time from, for example, optical phonons to acoustic phonons. Since the emphasis of this work is on the model itself, this relaxation time for optical phonons to decay into acoustic phonons is taken as a constant 5 ps [24].

At this point, a solution for Eq. (2) can be found by assuming that the lattice is initialized at 300 K and that the values of the thermal conductivity at 300 K are used. The boundary condition at the bottom edge of the device is taken either as a Dirichlet condition at 300 K, or as the convective condition of Eq. (6) with the ambient temperature assumed to be a constant 300 K. All other boundary conditions are assumed to be perfectly insulating, i.e., a zero flux exiting in the direction normal to the boundary. In other words, we solve
$∇2Tμ,Ci=−Pμ(r)κμ,C(Ti−1)$
(16)
subject to the boundary condition
$−κ∇T=H(TB,C−300 K)$
(17)

where $TB,C$ denotes the temperature in each respective boundary cell C along the bottom of the device.

Fig. 6
Fig. 6
Close modal
Fig. 7
Fig. 7
Close modal

The scattering rates are then updated in every cell of the simulation domain according to the local temperature, and the process is repeated until both the optical and acoustic temperatures have converged to a steady-state value everywhere. The described simulation procedure is depicted in the flowchart of Fig. 6. Note that the process depicted in the thermal solver is solved self-consistently in an inner loop, and then convergence is tested based on change between the current temperature solution found and the previous one.

Fig. 8
Fig. 8
Close modal

Here, we have chosen the GaN/AlGaN structure characterized experimentally by Altuntas et al. [25], and previously simulated using the Kirchhoff transformation by our group [26]. The layout of this device is shown in Fig. 7, consisting of a T-gate structure with a relatively large 1-μm buffer layer. A thin 1-nm AlN barrier is used to create a 2DEG conduction channel within the device due to a triangular well that is formed because of the conduction band mismatch of the heterojunction. In addition, thread dislocations are included in the simulations at a 2D density of $TDD=5×109$ cm−2.

Fig. 9
Fig. 9
Close modal

The final temperature maps obtained at the bias point VDS = 2 V, VGS = 10 V, where both high currents and high electric fields are present, are shown in Figs. 8 and 9, respectively. Here, the temperatures are shown as a color map on top of the 3D surface representing the conduction band to illustrate where the gate, source, drain, and channel regions are located in the device. These can be directly compared to the previous result obtained using the Kirchhoff transformation in Ref. [26], and we see good agreement between the two electro-thermal results.

Fig. 10
Fig. 10
Close modal
Fig. 11
Fig. 11
Close modal

That is, the optical mode temperature peaks at near 950 K and is highly localized to the drain edge of the gate as should be expected. This is because the optical phonons have very low group velocity, and hence a low thermal conductivity. Hence, the dominant decay path available to them is to decay into acoustic phonons, which generally happens at a rate slower than new optical phonons are generated. On the other hand, the acoustic temperature map is much more spread because of their high thermal conductivity, and the peak acoustic temperature is seen to be roughly 380 K, also to the drain side of the gate.

Next, we illustrate the effect of using both convective boundary conditions and convective interface conditions, as this is one of the main features of the proposed iterative method.

In Fig. 10, we show the IV curves for three separate cases of modeling the heat sink; (1) using a simple 300 K Dirichlet boundary condition, (2) using a convective boundary condition with H =108 W/m2 K, and (3) using a convective boundary condition with H =107 W/m2 K. Note that the perfect heat sink represented by the Dirichlet boundary condition would be represented by $H→∞$, and also that the most extreme condition of H =107 W/m2 K is an unrealistically poor heat sink, as evidenced by the extremely large temperatures seen in Fig. 11.

As the value of the heat transport coefficient, H, of the heat sink decreases, its ability to dissipate the heat generated in the device is reduced, and we see an associated large decrease in current due to increase in electron–phonon scattering rates with temperature. In addition, the current degradation is much more pronounced going from H =108 W/m2 K to H =107 W/m2 K than is seen between the Dirichlet and H =108 W/m2 K boundary conditions. This is due to the large temperature increase when using the H =107 W/m2 K heat sink, which in turn greatly increases the scattering rates seen in the device leading to reduced current. The ability to simulate nonideal thermal contacts is highly valuable as their effect on electrical characteristics are significant even for high values of H, as seen in Fig. 10.

## 4 Conclusion

We have demonstrated the validity of a proposed iterative approach for solving an energy balance equation to include temperature-dependent properties within a Monte Carlo simulation. The iterative approach was chosen for two reasons: (1) it overcomes the drawbacks of using the Kirchhoff transformation, and (2) it allows us to solve a Poisson-like equation. The importance of solving a Poisson-like equation is that any device simulator already must have an efficient solver for Poisson's equation to obtain the electrostatics in the device from the charge distribution. Hence, this allows for easier implementation within an existing finite difference device simulator, and for reuse of large amounts of already existing work. Within the iterative approach, deviations from the ideal case such as convective boundary conditions and material interface conditions are easily modeled while incorporating highly realistic nonlinear material properties such as the thermal conductivity and electron–phonon scattering rates. It is seen that this approach works extremely well for thermal larger-scale thermal spreading problems as compared to full nonlinear FEA results, as well as at the device scale when coupled to the existing CMC framework for full nonlinear electro-thermal simulation of a GaN/AlGaN HEMT.

## Funding Data

• Air Force Office of Scientific Research (AFOSR) (Grant No. FA9550-16-1-0406; Funder ID: 10.13039/100000181).

## References

1.
Artaki
,
M.
, and
Price
,
P. J.
,
1989
, “
Hot Phonon Effects in Silicon Field-Effect Transistors
,”
J. Appl. Phys.
,
65
(
3
), pp.
1317
1320
.10.1063/1.343027
2.
Semenov
,
O.
,
Vassighi
,
A.
, and
Sachdev
,
M.
,
2006
, “
Impact of Self-Heating Effect on Long-Term Reliability and Performance Degradation in CMOS Circuits
,”
IEEE Trans. Device Mater. Reliab.
,
6
(
1
), pp.
17
27
.10.1109/TDMR.2006.870340
3.
Fischetti
,
M. V.
, and
Laux
,
S. E.
,
1988
, “
Monte Carlo Analysis of Electron Transport in Small Semiconductor Devices Including Band-Structure and Space-Charge Effects
,”
Phys. Rev. B
,
38
(
14
), pp.
9721
9745
.10.1103/PhysRevB.38.9721
4.
Hess
,
K.
,
1991
,
Monte Carlo Device Simulation: Full Band and Beyond
,
Springer
,
New York
.
5.
Jacoboni
,
C.
, and
Reggiani
,
L.
,
1983
, “
The Monte Carlo Method for the Solution of Charge Transport in Semiconductors With Applications to Covalent Materials
,”
Rev. Mod. Phys.
,
55
(
3
), pp.
645
705
.10.1103/RevModPhys.55.645
6.
Jacoboni
,
C.
, and
Lugli
,
P.
,
1989
,
The Monte Carlo Method for Semiconductor Device Simulation
,
Springer
,
Vienna, Austria
.
7.
Saraniti
,
M.
, and
Goodnick
,
S.
,
2000
, “
Hybrid Fullband Cellular Automaton/Monte Carlo Approach for Fast Simulation of Charge Transport in Semiconductors
,”
IEEE Trans. Electron Devices
,
47
(
10
), pp.
1909
1916
.10.1109/16.870571
8.
Wachutka
,
G.
,
1990
, “
Rigorous Thermodynamic Treatment of Heat Generation and Conduction in Semiconductor Device Modeling
,”
IEEE Trans. Comput.-Aided Des. Integr. Circuits Syst.
,
9
(
11
), pp.
1141
1149
.10.1109/43.62751
9.
Fushinobu
,
K.
,
Majumdar
,
A.
, and
Hijikata
,
K.
,
1995
, “
Heat Generation and Transport in Submicron Semiconductor Devices
,”
ASME J. Heat Transfer
,
117
(
1
), p.
25
.10.1115/1.2822317
10.
Raleva
,
K.
,
Vasileska
,
D.
,
Goodnick
,
S. M.
, and
Nedjalkov
,
M.
,
2008
, “
Modeling Thermal Effects in Nanodevices
,”
IEEE Trans. Electron Devices
,
55
(
6
), pp.
1306
1316
.10.1109/TED.2008.921263
11.
,
B.
,
Vasileska
,
D.
, and
Goodnick
,
S. M.
,
2012
, “
Is Self-Heating Responsible for the Current Collapse in GaN HEMTs?
,”
J. Comput. Electron.
,
11
(
1
), pp.
129
136
.10.1007/s10825-012-0385-z
12.
,
T.
,
Kelsall
,
R. W.
, and
Pilgrim
,
N. J.
,
2006
, “
Investigation of Self-Heating Effects in Submicrometer GaN/AlGaN HEMTs Using an Electrothermal Monte Carlo Method
,”
IEEE Trans. Electron Devices
,
53
(
12
), pp.
2892
2900
.10.1109/TED.2006.885099
13.
Hao
,
Q.
,
Zhao
,
H.
, and
Xiao
,
Y.
,
2017
, “
A Hybrid Simulation Technique for Electrothermal Studies of Two-Dimensional GaN-on-SiC High Electron Mobility Transistors
,”
J. Appl. Phys.
,
121
(
20
), p.
204501
.10.1063/1.4983761
14.
Hao
,
Q.
,
Zhao
,
H.
,
Xiao
,
Y.
, and
Kronenfeld
,
M. B.
,
2018
, “
Electrothermal Studies of GaN-Based High Electron Mobility Transistors With Improved Thermal Designs
,”
Int. J. Heat Mass Transfer
,
116
, pp.
496
506
.10.1016/j.ijheatmasstransfer.2017.09.048
15.
Bonani
,
F.
, and
Ghione
,
G.
,
1995
, “
On the Application of the Kirchhoff Transformation to the Steady-State Thermal Analysis of Semiconductor Devices With Temperature-Dependent and Piece-Wise Inhomogenous Thermal Conductivity
,”
Solid-State Electron.
,
38
(
7
), pp.
1409
1412
.10.1016/0038-1101(94)00255-E
16.
Hahn
,
D. W.
, and
Ozisik
,
M. N.
,
2012
,
Heat Conduction
,
Wiley
,
Hoboken, NJ
.
17.
Bagnall
,
K. R.
,
Muzychka
,
Y. S.
, and
Wang
,
E. N.
,
2014
, “
Application of the Kirchhoff Transform to Thermal Spreading Problems With Convection Boundary Conditions
,”
IEEE Trans. Compon., Packag. Manuf. Technol.
,
4
(
3
), pp.
408
420
.10.1109/TCPMT.2013.2292584
18.
Palankovski
,
V.
, and
Selberherr
,
S.
,
1999
, “
Thermal Models for Semiconductor Device Simulation
,”
Third European Conference on High Temperature Electronics
(
HITTEN 99
), Berlin, July 7, pp.
25
28
.10.1109/HITEN.1999.827343
19.
Garven
,
M.
, and
Calame
,
J.
,
2009
, “
Simulation and Optimization of Gate Temperatures in GaN-on-SiC Monolithic Microwave Integrated Circuits
,”
IEEE Trans. Compon. Packag. Technol.
,
32
(
1
), pp.
63
72
.10.1109/TCAPT.2008.2004586
20.
Hickson
,
R.
,
Barry
,
S.
,
Mercer
,
G.
, and
Sidhu
,
H.
,
2011
, “
Finite Difference Schemes for Multilayer Diffusion
,”
Math. Comput. Modell.
,
54
(
1–2
), pp.
210
220
.10.1016/j.mcm.2011.02.003
21.
The Mathworks, Inc.
,
2018
, “
MATLAB Release 2018
,”
The Mathworks
, Natick, MA.
22.
Yamakawa
,
S.
,
Aboud
,
S.
,
Saraniti
,
M.
, and
Goodnick
,
S. M.
,
2003
, “
Fast Full-Band Device Simulator for Wurtzite and Zincblende GaN MESFET Using a Cellular Monte Carlo Method
,”
J. Comput. Electron.
,
2
(
2–4
), pp.
481
485
.10.1023/B:JCEL.0000011475.74817.6e
23.
Ambacher
,
O.
,
Majewski
,
J.
,
Miskys
,
C.
,
,
A.
,
Hermann
,
M.
,
Eickhoff
,
M.
,
Stutzmann
,
M.
,
Bernardini
,
F.
,
Fiorentini
,
V.
,
Tilak
,
V.
,
Schaff
,
B.
, and
Eastman
,
L. F.
,
2002
, “
Pyroelectric Properties of Al(in)GaN/GaN Hetero- and Quantum Well Structures
,”
J. Phys.: Condens. Matter
,
14
(
13
), pp.
3399
3434
.10.1088/0953-8984/14/13/302
24.
Tsen
,
K. T.
,
Ferry
,
D. K.
,
Botchkarev
,
A.
,
Sverdlov
,
B.
,
,
A.
, and
Morkoç
,
H.
,
1997
, “
Direct Measurements of Electron-Longitudinal Optical Phonon Scattering Rates in Wurtzite GaN
,”
Appl. Phys. Lett.
,
71
(
13
), p.
1852
.10.1063/1.119420
25.
Altuntas
,
P.
,
Lecourt
,
F.
,
Cutivet
,
A.
,
Defrance
,
N.
,
,
E.
,
Lesecq
,
M.
,
Rennesson
,
S.
,
Agboton
,
A.
,
Cordier
,
Y.
,
Hoel
,
V.
, and
De Jaeger
,
J.-C.
,
2015
, “
Power Performance at 40 GHz of AlGaN/GaN High-Electron Mobility Transistors Grown by Molecular Beam Epitaxy on Si(111) Substrate
,”
IEEE Electron Device Lett.
,
36
(
4
), pp.
303
305
.10.1109/LED.2015.2404358
26.
Latorre-Rey
,
A. D.
,
Merrill
,
K.
,
Albrecht
,
J. D.
, and
Saraniti
,
M.
,
2019
, “
Assessment of Self-Heating Effects Under Lateral Scaling of GaN HEMTs
,”
IEEE Trans. Electron Devices
,
66
(
2
), pp.
908
916
.10.1109/TED.2018.2888812