Proceedings of ASME 2012 Summer Heat Transfer Conference
July 8-12, 2012, Puerto Rico, USA
Angie Fan, Calin Tarau & Richard Bonner – Advanced Cooling and Technologies Inc. Lancaster PA.
Tomas Palacios – Massachusetts Institute of Technology Cambridge, MA, USA
Massoud Kaviany – University of Michigan Ann arbor, MI, USA
In this study, a thermal and electrical coupled device solver is developed to simulate the energy transfer mechanism within a GaN FET with a gate length of 0.2µm. The simulation simultaneously solves a set of hydrodynamic equations (derived from the Boltzmann Transport Equation) and the Poisson equation for electron, optical phonon and acoustic phonon energies, electron number density, electric field and electric potential. This approach has been previously established for gallium arsenide (GaAs) devices [36,37], but has not been extended to GaN due to the lack of readily available property values for GaN devices that are required. Via extensive literature study, high-fidelity properties for GaN were collected in analytical forms with respect to many dependencies, e.g. lattice temperature, electrical field, electron number density, doping rate, defects rate. These properties are then implemented into the developed code to provide a high accuracy sub-micron GaN device simulation.
Simulations show that non-equilibrium heat generation is exhibited in a typical device while the drain current is reduced due to the decrease in electron mobility. Future analysis is needed to quantify the hot-electron effect on reducing the drain current and to discover more effective ways of heat removal.
Gallium nitride (GaN) structures have out standing electronic properties. GaN devices operate at higher voltages, higher frequencies, and higher power densities than other devices due to its inherently high critical electric field, electron mobility, and electron saturation velocity. However, the
relatively poor thermal conductivity of GaN and higher heat flux operation make thermal management an issue. Moreover, in sub-micron sized GaN devices, interactions between the electrical and thermal properties become significant due to the strong electron-phonon coupling near the device junction. Therefore, in order to accurately simulate the electrical or thermal performance of a GaN device, the thermal and electrical governing equations need to be solved concurrently.
The thermal energy transfer mechanism near the junction in GaN devices is via scattering of hot electrons with longitudinal optical (LO) phonons, followed by LO-phonon decay to longitudinal acoustic (LA) phonons, which ultimately transport the heat through the substrate. Dissipation of the heat is dominated by the lifetime of LO-phonon. The LO-phonon lifetime is related to electron density, hot electron temperature, and supplied electric power. The strong coupling of electrical and thermal characteristics suggests that GaN properties need to be expressed with many dependencies.
On the device simulation side, most commercially available simulators have been developed over two decades based on Si and GaAs. In the past, the power densities of the devices were not high enough to cause the thermal characteristics to interfere strongly with the electrical behavior. Therefore, the drift-diffusion model was sufficient, which assumes constant temperature while calculating the electrons and holes transport rates. Material properties were also simplified since they were independent of the thermal characteristics. Temperature change was calculated independently using the Fourier’s law of heat conduction. These assumptions are no longer valid with the advent of the high-power devices. At high power densities the thermal and electrical characteristics depend on each other and must be treated as a coupled problem.
The aim of this study is to consider the nonequlibrium nature of the electrons and the phonons and demonstrate the appropriate thermal and electrical characteristics in submicron FETs. We choose MESFET made of GaN for our study due to its simple structure and representation in the signal processing and communication device families.
In section 3, the governing equations of the non-equilibrium semiconductor devices are discussed. Section 4 shows the material properties study. Section 5 discusses the computational method of this work, and Section 6 discusses the results. Finally, Section 7 states the conclusions.
1. GOVERNING EQUATIONS
The Boltzmann’s transport equation (BTE) is valid for general inhomogeneous materials with arbitrary band structure . To solve this equation numerically by discretization of the differential and integral operators in seven-dimensional space is computationally very expensive. A widely used numerical method for solving the BTE is the Monte Carlo (MC) method. This method has been proven to give accurate results, but is still computationally expensive. A good engineering-oriented approach is to transform the BTE into a macroscopic transport model as long as meaningful expressions of electron transport in terms of average electron velocity and effective electron temperature are valid. Blotekjaer (1970) derived conservation equations by taking the moments of the BTE using the weight functions without imposing any assumptions on the form of the distribution function. The resulting equations represent the electron charge, momentum, and energy conservation, respectively . They are
where n is the electron number density, v, is electron is electron drift velocity, E is electrical field, and Te is electron temperature. The electron momentum density p and energy density We can be written as
The last terms in the equations are collision terms. A macroscopic relaxation time approximation is employed to model the collisions as 
which introduces the momentum and electron energy relaxation times Tm and Te-LO. A discussion on this approximation is given in . For example, Eq. (6) represents that the electron generation, recombination, and impact ionization in this article are neglected for simplicity. The electron momentum conservation Eq. (2) is the analog of the momentum equation of the Navier-Stokes equations for fluid flow. The driving forces in Eq. (2) are the electric field, the electron number density, and the electron pressure (a.k.a temperature gradients), while the drag force is electron-phonon collision.
For the heat flow vector Q is Eq. (5), the Fourier law is applied,
where Ke is the electron thermal conductivity. The Fourier law may not be the best way to get a close form. A more accurate method is to consider the third moment of the BTE. But more unknown microscale coefficients such as the relaxation times will be introduced. Under a high electric field E > 20 kV/cm, the electrons become energetic enough to produce optical phonons. They first lose energy to optical phonons and optical phonons decay to acoustic phonons. The process can be described using the first law of thermodynamics. The resulting energy conservation equations for optical and acoustic phonons are
where WLO and Wa are optical and acoustic phonon energy densities, which are related to the heat capacity C as dW=CdT, and Ka is the lattice thermal conductivity. The optical phonon collision in Eq. (10) and Eq. (11) is expressed as
where Top-a is optical phonon relaxation time, TLO and Ta are the optical and acoustic phonon temperatures. If E < 20 kV/cm, the electrons only interact with acoustic phonons, thus the energy conservation equation become
By applying Eq. (4,5,6,7,8,9,12) into Eq. (1,2,3,10,11) and rearranging the equations, also including the Poisson equation for the electric field, the closed governing equations of the submicron semiconductor devices when E > 20 kV/cm become
To solve the governing equations numerically, the expressions for the electron mobility, the electron drift velocity, the thermal conductivities for electron and phonons, the heat capacity for the phonons, the momentum and energy relaxation times, and the electron effective mass in GaN are required. An extensive literature study is performed to obtain all the parameters. The results are presented in the following section.
Optical and acoustic phonon heat capacity
The optical phonon (Einstein) heat capacity and the acoustic phonon (Debye) heat capacity have the following definitions
where NA is the Avogadro number, θE is the Einstein temperature, θA is the Debye temperature, and T is the lattice temperature.
The bulk lattice heat capacity is the sum of both the heat capacities:
The bulk heat capacity is experimentally measurable. The two components, CpA and CpLO, can be found by guessing the characteristic temperatures to fit the experimental bulk value.
In Figure 1, measurements of the bulk lattice heat capacity of the hexagonal GaN by several authors [4,5,6,7,8,9] are displayed. Danilchenko’s measurement was selected for fitting due to the broadest temperature range of 5-300K. The best fit with the accuracy of 3 % was obtained for Debye’s temperature KA365=θand Einstein’s temperature KE880=θ. The results are shown in Figure 1. The agreement between the experimental data obtained by these authors for the GaN heat capacity and the fit curve based on the sum of the two components (shown as Danilcenko in Figure 1) is reasonable especially in the 100-700 K temperature range. In conclusion, our simulations use Eq. 6 and Eq. 7 for the optical and acoustic phonon heat capacities where KD365=θandKE880=θ.
Electron drift velocity
The electron velocity in GaN increases for low electric fields and then becomes quasi-saturated close to 20 kV/cm10. After this quasi-saturation, the electron velocity increases again, reaching a maximum for an electric field of 100 kV/cm and then it decreases to its saturated value of 1.2 x 107 cm/s. A very important difference between the high field transport model of GaN and that of other semiconductors is that in the latter, there is no quasi-saturation of the electron velocity. The quasi-saturation behavior is generally associated with the onset of optical phonon emission in nitride based semiconductors [11,12]. In other semiconductors like Si or GaAs, the energy for optical phonon emission is much lower and the quasi-saturation is not very different from the standard saturation. Some authors have also associated the early saturation of the electron velocity in AlGaN/GaN structures with hot phonon scattering [13,14]. It is obvious that the need for a valid GaN electron transport model under high electric field is significant and literature offers both experimental and simulation data [15,16,17,18,19,20,21,22,23,24,25,26,27].
Schwierz (2005) proposed a complete analytical model for the electron mobility in wurtzite GaN.
where µo is the low-field mobility, Vsat, Ec, n1 n2 n3, µmin, µmax, nref and α are fitting parameters. µmax represents the mobility of undoped or unintentionally doped samples where lattice scattering is the main scattering mechanism. µmin is the mobility in highly doped material, where impurity scattering is dominant. The parameter α is a measure of how quickly the mobility changes from µmin to µmax – nref is the carrier concentration at which the mobility is half way between µmin to µmax. The temperature dependence of the low-field mobility can be modeled by making the four fitting parameters in Eq. 23 temperature dependent according to
where Par is the parameter of interest (i.e. µmin, µmax, nref and α), Par0 is the value of the parameter Par at room temperature. For the rest parameters Vsat, Ec, n1 n2 n3 in Eq. 22, the temperature dependent relationship can be derived by
where a, b, and c are constants that have to be determined by fitting.
Through a series development of the parameters, the final expression for the drift velocity and the electron mobility become
Thermal conductivity–phonon and electron
Thermal conductivity of GaN mainly consists of two components: KA– the acoustic phonon thermal conductivity (or lattice thermal conductivity) and Ke– the electron thermal conductivity.
In a theoretical investigation  of the thermoelectric effects in wurtzite GaN crystals the two main contributions to the GaN thermal conductivity were evaluated. The dependence of the phonon contribution (KA) to the thermal conductivity on the dislocation and point defect densities has been studied in detail by Zou  and co-workers . It has been established that the room temperature thermal conductivity in wurtzite GaN is sensitive to the dislocation line density Ndis only for very large values of Ndis>5×1010 cm-2. At this high dislocation range the variation of the dislocation line density by two orders of magnitude can bring about a large (factor of two) variation in the thermal conductivity value. For GaN materials with lower dislocation density, the room temperature thermal conductivity is more sensitive to the concentration of point defects impurities, dopants, etc. For example, the change in the doping density can lead to a variation of A k from 1.77 W/cm-K to 0.86 W/cm-K at 300 K. In their work they selected Ndis= 5 x 1010 cm -2 and determine kA from the virtual-crystal model . It is well known that the kA contribution is much larger than Ke for most practical situations even at very large doping densities.
Eq. 29~31 are used in our device simulator. It is observed that electron number density is not among the dependencies for the thermal conductivity of the electron. However, it is shown (Liu-2005) that electron thermal conductivity becomes a significant contribution for the overall GaN thermal conductivity only for doping rates larger than 1019 cm-3. Since in our simulations the doping rates are smaller, the lack of electron number density as dependency is not a concern.
Electron effective mass
Various methods were used to find the effective mass of the electron in GaN and the resulted electron effective mass considered in our device simulator for GaN is:
Electron momentum relaxation time
The electron momentum relaxation time can be calculated from the relation 
Since the electron mobility is a function of electric field, electron number density and temperature (as seen in Eq. 27), the electron momentum relaxation time τm will have the same dependencies and the analytical representation is:
where µ(E, n, T) is the electron mobility described by Eq. 27.
Electron- optical phonon energy relaxation time
Electron-longitudinal optical phonon scattering rates in wurtzite GaN have been directly measured  by sub- picosecond time-resolved Raman spectroscopy. It was found that the total electron-LO phonon scattering rate is given by Γ=(4 ±0.8)x1013 s-1. Since the electron-LO phonon scattering rate is about Γ=5×1012 s-1 in GaAs  the observed total electron-LO phonon scattering rate in GaN is almost one order of magnitude larger than that in GaAs. One of the possible explanations for this enormous increase of electron-LO phonon scattering rate in GaN is its much larger ionicity as reflected from the much larger LO-TO phonon splitting. In general, the strength of the electron-LO phonon coupling is set by the lattice-dipole interactions. This splitting is directly proportional to the lattice polarization, and is a measure of the effective charge. In GaAs, the LO-TO splitting is about 3 meV, while it is about 25 meV in wurtzite GaN. However, in the two-mode system of hexagonal GaN, this must be split between the two modes to compute the contribution of each of the dielectric function. This together with the much smaller dielectric constant of GaN, lead to an expected increase of a factor of 9 in the scattering strength, which is quite close to that observed experimentally. It is estimated  that the effective charge of GaN is about twice that of GaAs, while the dielectric constant is almost half. These two factors above would suggest a factor of 8 increase in the electron-LO phonon scattering strength. This provides a consistent interpretation of the scattering increase as due to the ionicity of GaN relative to GaAs. In conclusion, the total electron-LO phonon scattering time is:
3. NUMERICAL STUDY OF A SUB-MICRON GAN MESFET
Metal-semiconductor field-effect transistors (MESFET) Figure 2 are usually made of III-V materials and commonly used in high-speed communication devices. Heat generation and transport in gallium arsenide (GaAs) MESFETs using the hydrodynamic method has been studied [36,37]. In this work, the heat generation and transport in 2D GaN MESFET with the updated properties is presented.
The computational domain of the 2D transistor is shown in Figure 2. Finite difference method is used to discretize the set of governing equations. Second order central difference scheme is applied to all terms except for the second term on the left- hand side of Eq. 16, where a first-upwind scheme is applied. Square mesh is applied to the entire domain with a uniform grid size of 5 nm. A constant voltage is applied between the drain and source contacts. Zero voltage and a Schottky barrier height VB8.0−=φare applied on the gate contact. At the rest of the boundary, zero gradient normal to the boundary is used. For electron number density DNn=is used at each of the source and the drain contacts, and zero particle flux normal to the boundary is applied at the rest of the boundary. For LO- phonon and acoustic phonon temperatures, constant heat transfer coefficient h=103W/m2 – K is given at the top and bottom boundaries. At the rest boundaries, an adiabatic condition is applied to simulate the shoulder to shoulder chip packaging arrangement. For electron temperature, an adiabatic condition is applied to the entire boundary. Implicit time-marching scheme is applied for all five equations. The LU Decomposition Conjugate Gradient Squared (LUCGS) method is used to solve the matrix. Two time steps of 0.01 ps and 4.5 ps are used for iteration. The first time step is used at the beginning of the iteration when the electron temperature responds rapidly to the applied field. Once the electron temperature reaches a quasi-steady state with very little change, a larger time step is used.
4. RESULTS AND DISCUSSIONS
Figure 3 gives the computational results for the voltage, electric field, electron mobility and reduced electron number density in the GaN device where Vsd=10V, Vg=0V, t=17ps.
Figure 3(a) shows the voltage distribution. Although there is no voltage applied to the gate contact, underneath the gate there is a narrow zone of negative potential. This is due to the Schottky barrier. It is also shown that the major voltage drop occurs at the drain side, which means higher electric field will occur in this region. This is confirmed by the plot of the electric field distribution, Figure 3(b). The highest electric field locates under the gate at the drain side.
Figure 3(c) shows the computed profile of the electron mobility. There is a large zone of low mobility under the gate and drain. This zone coincides with the major voltage drop zone in Figure 3(a). The physics behind this is that the strong hot-electrons and LO-phonons scattering occurring in this region due to the extremely high electric field makes the electron velocity saturate. The electron mobility model adopted here depends not only on the average energy but also on the electric field and the electron number density; therefore it prevents the velocity overshoot  in the non-equilibrium device.
Figure 3(d) shows the reduced electron number density n/ND. A wide depletion zone is clearly shown under the gate. Therefore, the current channel is narrow which leads to large current density.
Figure 4 shows the concurrently computed heat generation rate, electron temperature, LO-phonon temperature and acoustic phonon temperature.
Figure 4(a) shows the heat generation rate. The heat generation rate is the product of the current density and the electric field. The current density and the electric field in the gate channel at the drain side are higher; therefore the highest heat generation rate appears around the region.
Figure 4(b) shows the distribution of the electron temperature. The hot electrons appear under the gate at the drain side corresponding to the high electric field. The peak electron temperature is 1400 K. The electron temperatures in the rest of the devices remain around 300 K. Experimental measurement is needed for verification.
Figure 4(c) (d) show the LO-phonon and acoustic phonon temperatures. The highest temperature region appears under the drain. Comparing to the electron temperature, the two lattice temperatures rise much slower. For the given computing time (t=17ps), the lattice temperatures are slightly increased. It is understandable that the electrons pick up the energy much more efficient from the applied electric field than transfer the energy to the lattice. But experimental data need to be compared in the future work.
In this work, a 2-D hydrodynamic equations based GaN device simulator is established to capture the coupling effects of the electrical and thermal characteristics in high power GaN devices. High fidelity material properties are established from an intensive literature study. It reveals the unique electron transfer feature in GaN. Future work includes validate the simulator with broad experimental data, improve the computational efficiency, and discover more effective ways of heat removal.
This work was supported by the National Reconnaissance Office’s (NRO) Applied Technology Department (ATD) through the Jalapeno program under contract NRO000-09-C-0348.
 Tibor Grasser, Ting-wei Tang, Hans Kosina, and Siegfried Selberherr, “A Review of Hydrodynamic and Energy-Transport Models for Semiconductor Device Simulation”, Proceedings of the IEEE, vol. 91, No. 2, Feb. 2003
 Jie Lai, and Arun Majumda, “Concurrent Thermal and Electrical Modeling of Submicrometer Silicon Devices”, J. Appl. Phys. 79 (9), 1May 1996
 M. Lundstrom, Fundamentals of Carrier Transport. Reading, MA:Addison-Wesley, vol. 10, Modular Series on Solid State Devices, 1990
 Danilchenko, B. A. , T. Paszkiewicz, S. Wolski, A. Jeżowski, and T. Plackowski, “Heat capacity and phonon mean free path of wurtzite GaN”, Appl. Phys. Lett. 89, 2006
 Chen, X., Lan, Y., Liang, J., Cheng X., Xu Y., Xu T., Jiang P., and Lu K., “Structure and Heat Capacity of Wurtzite GaN from 113 to 1073 K”, Chinese Phys. Lett. 16, 1999
 Itagaki, K. et.al., Thermochim. Acta, 163, 1990
 Leitner, J., Strejc, A., Sedmidubský, D., Růžička, K., “High temperature enthalpy and heat capacity of GaN”, Thermochimica Acta, Volume 401, Issue 2, Pages 169–173, 2003
 Ziborak-T I., et.al., Journal of Thermal Analysis and Calorimetry, 91, 2008
 Sanati, M., Albers, R. C., and Pinski, F., J., “ω-phase formation in NiAl and Ni2Al alloys”, J. Phys.: Condens. Matter 13, 2001
 Palacios T. et.al. Nitride Semiconductor Devices, Piprek, 2007
 Singh M., J. Appl. Phys. 94, 2003
 Gelmont B., J. Appl. Phys. 74, 1993
 Ardaravicius L. et.al. Appl.Phys.Lett.83, 2003
 Ridley B.K. et.al., J. Appl. Phys. 96, 2004
 Caughey D.M.,et.al., Proc. IEEE 55, 1967
 Benbakhti B. et. al., J. Phys. 193, 2009
 Kabra S. et.al. Microelectronics Journal 37 2006
 Schwierz F., Solid-State Electronics 49, 2005
 Mnatsakanov T.T., et. al. Solid-State Electron, 47, 2003
 Farahmand M., et.al. Solid-State Electronics, 49, 2005
 Arora N.D. et.al. IEEE Trans Electron Dev, 29, 1982
 Farahmand M., et.al. Solid-State Electronics, 49, 2005
 Schwierz F.: A compilation of more than 200 references. 2004
 Look D.C., et.al. Appl Phys Lett, 79, 2001
 Schwierz F., Solid-State Electronics, 49, 2005
 Mnatsakanov T.T., et. al. Solid-State Electron, 47, 2003
 Farahmand M., et.al. Solid-State Electronics, 49, 2005
 Liu W., et al., J. Appl. Phys. 97, 2005
 Zou J., et al., J. Appl. Phys. 92, 2002
 Kotchetkov D., et al., Appl. Phys. Lett. 79, 2001
 Liu W., et.al., Appl. Phys. Lett. 85, 2004
 Kittel, 1986
 Tsen K.T., et. al., Appl. Phys. Lett., 71, 1997
 Kash, J. A. , et al. Phys. Rev. Lett. 54, 1985
 Phillips, J. C., et.al. Bands in Semiconductors ~Academic, New York, 1973
 Majumdar A., Fushinobu K., and Hijikata, K., “Effect of gate voltage on hot-electron and hot-phonon interaction and transport in a submicrometer transistor”, J. Appl. Phys., 21, 1995
 Fushinobu K., Majumdar, A., and Hijikata, K., “Heat generation and transport in submicron semiconductor devices”, J. of Heat Transfer, 31, 1995