## 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 [3–6]. 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 *I*–*V* 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.

*μ*is written

where $W\mu (r,t)=(1/\Omega )\Sigma kE\mu (k)f\mu (r,k,t)$ is the ensemble energy density within the volume Ω of the reciprocal space, $F\mu (r,t)=(1/\Omega )\Sigma kv(k)E\mu (k)f\mu (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 (*e*–*p*) and phonon–phonon (*p*–*p*) interactions, respectively.

where Fourier's law has been used to approximate the energy flux, $F\mu (r)=\u2212\kappa \mu (r,T)\u2207T$. 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 $\kappa \mu ,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.

*T*

_{ref}is a chosen reference temperature, often 300 K, at which the thermal conductivity value, $\kappa \mu ,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]

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.

*R*

_{th}, or a finite heat transfer coefficient,

*H*, where

*H*=

*1/*

*R*

_{th}. Clearly, as $Rth\u21920,\u2009H\u2192\u221e$ 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

which relates the heat flux approaching the interface, $\u2212\kappa \u2207T$, with the heat transport coefficient and the temperature jump across the interface, $(TL\u2212TR)$. Since the Kirchhoff transformation is nonlinear for any nonlinear functional form of the thermal conductivity, it is generally not true that $H(\Theta L\u2212\Theta R)=H(TL\u2212TR)$. 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(\Theta L\u2212\Theta R)\u2248H(TL\u2212TR)$. 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.

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 $\alpha =\u22121.3$ [18]. In other words, the temperature dependent thermal conductivity is taken as $\kappa \mu (T)=150(T/300)(\u22121.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 × 10^{6} W/m^{2}, while a convective boundary is set along the bottom of the material characterized by a heat transfer coefficient *H *=* *10^{5} W/m^{2} K and an ambient temperature $T\u221e=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.

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 × 10^{9} W/m^{2} 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.

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.

### 3.2 Interface Resistance.

*x*-axis, as it is a straightforward generalization in the other directions, the second derivatives on each side of the interface are written

*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, Δ*

*the distance between point*

_{xL}*i*− 1 and the interface, Δ

*the distance between point*

_{xR}*i*+

*1 and the interface, Δ*

*the distance between points*

_{xw}*i*− 1 and

*i*− 2 to the left, Δ

*the distance between points*

_{xe}*i*+

*1 and*

*i*+

*2 to the right. In Eqs. (8) and (9), we have four unknowns: $Ti\u22121,\u2009Ti+1,\u2009TL$, and $TR$. The temperatures at the grid points everywhere in the simulation, such as $Ti\u22121$ 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*

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 $\kappa A=1$ W/m K, $\kappa B=0.1$ W/m K, and the heat transport coefficient at each interface is *H _{i}* = 0.5 W/m

^{2}K.

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 10^{10} W/m^{2}. 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\u221e=293$ K.

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 $\alpha =\u22121.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.

### 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].

where *C _{i}* is the volumetric heat capacity of the phonon mode, $T(i,j)$ the temperature of the respective mode, and $\tau (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].

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

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.

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\xd7109$ cm^{−2}.

The final temperature maps obtained at the bias point *V*_{DS} = 2 V, *V*_{GS} = 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.

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 *I*–*V* 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 *=* *10^{8} W/m^{2} K, and (3) using a convective boundary condition with *H *=* *10^{7} W/m^{2} K. Note that the perfect heat sink represented by the Dirichlet boundary condition would be represented by $H\u2192\u221e$, and also that the most extreme condition of *H *=* *10^{7} W/m^{2} 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 *=* *10^{8} W/m^{2} K to *H *=* *10^{7} W/m^{2} K than is seen between the Dirichlet and *H *=* *10^{8} W/m^{2} K boundary conditions. This is due to the large temperature increase when using the *H *=* *10^{7} W/m^{2} 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).