Modeling and simulation of bulk gallium nitride power semiconductor devices

Bulk gallium nitride (GaN) power semiconductor devices are gaining significant interest in recent years, creating the need for technology computer aided design (TCAD) simulation to accurately model and optimize these devices. This paper comprehensively reviews and compares di ff erent GaN physical models and model parameters in the literature, and discusses the appropriate selection of these models and parameters for TCAD simulation. 2-D drift-di ff usion semi-classical simulation is carried out for 2.6 kV and 3.7 kV bulk GaN vertical PN diodes. The simulated forward current-voltage and reverse breakdown characteristics are in good agreement with the measurement data even over a wide temperature range. C 2016 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http: // creativecommons.org / licenses / by / 4.0 / ). [http: // dx.doi.org / 10.1063 / 1.4948794]


I. INTRODUCTION
Gallium nitride (GaN) is a highly promising wide bandgap semiconductor material to succeed silicon in high frequency power electronics applications. 1-3 While the lateral high electron mobility transistor (HEMT) remains the most popular GaN device type, vertical bulk GaN power devices have also considerable interest in recent years. 4-11 A GaN HEMT fabricated on a heteroexpitaxial layer on a substrate of sapphire, silicon, or silicon carbide offers potentially low fabrication cost, but currently suffers from high threading dislocation density (∼10 10 cm −2 ) and the associated reliability concerns, owing to the inherent lattice mismatch between the epitaxial film and foreign substrate. Vertical bulk GaN power devices, which can circumvent the limitations of the lateral GaN HEMTs such as the current collapse phenomenon and lack of avalanche energy capability, have recently become a possibility with high-quality free-standing bulk GaN materials being increasingly available with a low defect density (10 4 -10 6 cm −2 ). 4-6 For both lateral and vertical GaN devices, technology computer aided design (TCAD) simulation is a vital tool to accurately predict the behavior and optimize the design of these devices, following the well-proven practice of the silicon industry.
TCAD simulation solves the Poisson and current continuity equations for both electrons and holes in two or three dimensions to provide terminal current-voltage characteristics and physical insights of the semiconductor devices. Selection of proper physical models and model parameters is vital for TCAD simulation. Physical models and model parameters pertinent to energy bandgap, incomplete ionization, electron and hole mobility, impact ionization and carrier recombinationgeneration in GaN are reported in the literature. Most model parameters for GaN are determined from Monte Carlo simulations and a few are available from experimental data. However, due to different computational methods and non-uniform selection of physical mechanisms for Monte Carlo simulations, these parameters have a wide spread of values. Also the experimental values cannot be used perśe as different growth techniques, implantation methods and other fabrication steps yield widely fluctuating numbers. Different Monte Carlo calculation methods like the full potential linearized augmented plane wave (FLAPW), 12 orthogonalized linear combination of atomic orbitals (OLCAO) 13 or the empirical pseudo potential method (EPM) 14 are used to determine the band structure and hole and electron effective masses. 12,[14][15][16][17] Similar Monte Carlo calculations are performed to determine the bulk electron [18][19][20][21][22][23][24][25][26][27] and hole [28][29][30][31] mobility employing different band structure and scattering mechanisms. Although experimentally determined values of electron mobility [32][33][34][35] are reported there are few available experimental data for hole mobility. The Chynoweth law 36 is accepted as an accurate representation of the avalanche effect in GaN but impact ionization coefficients reported from Monte Carlo simulations 24,25,27 as well as experimental results 36,37 differ from one work to another. In case of high doping concentrations incomplete ionization of donor (silicon) 38,39 and acceptor (magnesium) [40][41][42][43] needs to be taken into account and their corresponding activation energies also show a wide spread. Radiative constant 44 and Auger coefficients 45,46 for GaN are studied in detail for the blue LED and their values are well accepted, but general SRH recombination parameters like electron and hole lifetime are still an unexplored area mainly due to the randomness of the distribution of traps from crystal to crystal. A range of GaN models and parameters have been used for two-dimensional (2D) drift-diffusion modeling for GaN HEMTs 47-52 over the past decade, However, they showed significant inconsistency from one work to another, creating the need to comprehensively compare and calibrate the models and parameters for better simulation accuracy and consistency. Furthermore, there have been very few reports on modeling of vertical bulk GaN devices, which operate very differently from lateral GaN HEMTs. Many of the physical models implemented in lateral HEMTs such as piezoelectric polarization and the 2D electron gas are of little relevance with bulk vertical GaN devices. On the other hand, bulk electron mobility, bulk hole mobility and impact ionization are of particular interest for vertical architectures that need to be carefully addressed. Baik et al. 53 reported on modeling of bulk GaN Schottky and PN diodes using drift-diffusion simulations, but the models were not validated against measurement data. In this work, we provide a comprehensive review and comparison of the physical models and associated model parameters for bulk GaN materials in the literature, and discuss the appropriate selection of models and parameters for TCAD simulation. We model 2.6 kV and 3.7 kV bulk GaN vertical PN diodes with Synopsys Sentaurus TCAD tools 54 using selected model and parameter sets. The simulated forward current-voltage and reverse breakdown characteristics are compared with measured data. Good agreement is found between the simulation and experimental results even over a wide temperature range. It is worth noting that trap-related effects are not considered in our simulation due to the lack of detailed knowledge of trap distribution in these devices. The models and model parameters proposed in this work can be used in other TCAD software.

II. MODELS AND PARAMTERS
Material parameters for GaN vary considerably with the crystal variants as well as the orientation of the epitaxial layer and the growth technique employed. Bulk GaN devices are fabricated in the Wurtzite crystal structure with the epitaxial growth along the c-axis of the crystal by metal organic chemical vapor deposition (MOCVD). 4,7 Fundamental GaN material parameters relevant to a MOCVD grown epitaxial layer along the c-axis on the Wurtzite variant are used in the simulations and listed in Table I. TCAD simulations solve the Poisson equation and the electron and hole current continuity equations to predict the behavior of semiconductor devices under electrical stress as shown in Eq. (1)-(3) respectively.
dp(x, y, z,t) dt  where ε s is the semiconductor permittivity, ϕ the electrostatic potential, q the electron charge, n and p the free electron and hole density, N D − the ionized donor density, N A + the ionized acceptor density, J n/p the electron and hole current density, G n/p the rate of regeneration for electrons and holes and R n/p the rate of recombination for electrons and holes. The pertinent models are selected from a vast array of available models in the Sentaurus TCAD suite.

A. Bandgap and electron and hole effective mass
GaN being a direct bandgap material, the minimum of the conduction band lie directly above the valence band maxima in the E-k plane at the Γ-region where k = 0. The direct energy bandgap for GaN at 0 K is reported anywhere between 3.38 eV to 3.56 eV. [58][59][60][61][62] A median value of 3.53eV at 0 K is used in this work for the bandgap model which is equivalent to 3.46 eV at 300 K by using the Varshni expression 63 in Eq. (4).
where E g (0) is the forbidden energy gap at 0 K at the Γ-region in the E-K plane and α and β are the Varshni parameters. The Varshni parameters are compared from Refs. 58, 60, 62, and 64 and the values of α = 9.09x10 −4 eV.K −1 and β = 830 K are used from Ref. 64 as they give the best fit for the change of knee voltage in the forward conduction mode of the GaN P-N diode over a temperature range of 300-425 K. For TCAD simulations a single averaged hole and electron effective mass is needed. The threefold degeneracy in the valence band at K = 0, owing to the heavy hole valence band, the light hole valence band and the split-off valence band lead to a complicated effective mass calculation. Santic et al. 56 performed a detailed comparison of the hole effective masses from the literature and proposed a consolidated value of 1.25m 0 at room temperature. At k = 0, there is only a single minimum in the conduction band and the low field electron effective mass at room temperature has an accepted value of around (0.20 ± 2)m 0 . 12,14,16,17,62 However, effective mass of the electron changes at high electrical fields due to the migration from the low lying Γ-valley to the higher M-L valleys. The change of electron effective mass due to migration into higher M-L valleys at higher fields is modeled later in the mobility section.

B. Incomplete Ionization and Bandgap narrowing
GaN is doped with silicon (Si) as donor and magnesium (Mg) as acceptor species respectively, both of which have high activation energies and are not completely ionized. Especially for Mg-doped p-type GaN, proper modeling of incomplete ionization becomes extremely important as the doping efficiency of Mg drops down to 5-10%. The concentration of the ionized impurity atom for n-type Si-doped GaN and Mg-doped p-type GaN is incorporated in the incomplete ionization model in TCAD according to the Fermi-Dirac distribution shown in Eq. (5a) 38 where N A , N D are the original acceptor and donor concentrations, to the donor and acceptor activation energies, E c and E V being the conduction and valence band energy respectively. The degeneracy of the electronic state, g B has a value of 2 for GaN. N C and N V are the effective density of states in the conduction and valence band respectively calculated by equation (6).
where, N V /C are the valence and conduction band effective density of states, m h/e * the hole and electron effective mass and T is the temperature.
The donor and acceptor activation energies are effectively reduced by the high carrier concentration in semiconductors which is also incorporated in the incomplete ionization model. The dependence of ∆E D/ A on the carrier concentration is modeled according to Eq. (7) 39,43 where ∆E D,0 and ∆E A,0 are the donor and acceptor ionization energy at very low doping levels.
A value of 17 meV 38,66 for ∆E D,0 has been used in our simulations. ∆E A,0 however show a significant spread and ranges between 140 meV -245 meV. 41,43,67 Simulations performed using ∆E A,0 = 240 meV gives an effective Mg doping density of around 5.2x10 17 cm −3 for an initial doping of 1x10 19 cm −3 which concurs with the experimental data 4 and is used in our work. The value of α n and α p is 3.4x10 −6 meV-cm for Si-doped 39 and 3.14 x10 −5 meV. cm for Mg-doped 43 GaN respectively. Due density, a bandgap narrowing effect takes place and this is manifested by a cube root of the carrier density for Si-doped GaN which takes into account the exchange contribution of the electron-electron interaction. 39 For Mg-doped GaN an accurate bandgap narrowing (BGN) model is yet to be established, hence only a model for the Si-doped GaN is incorporated in the simulations. In n-type GaN, the BGN effect is modeled according to Eq. (8) 66 and is implemented by the Jain-Roulston model.
where ∆E g is the change in the bandgap energy, E g (0) the band gap energy of pure sample, E g (n) the doping dependent bandgap energy, K is a fitting parameter and n is the electron density. The fitting parameter K was estimated to be 1.15x10 −8 eV-cm. 39

C. Polarization
(0001)-oriented wurtzite GaN is a polar material having positive gallium (Ga) and negative nitrogen (N) atoms on two faces of the crystal. The two types of polarization manifested in Wurtzite GaN are piezoelectric polarization owing to mechanical perturbations in the crystal structure and spontaneous polarization owing to an intrinsic asymmetry of the gallium and nitrogen covalent bonding in the equilibrium crystal structure. The polarization effect in GaN is modeled by the set of Eq. (9a)-(9c). 68 where P is the total polarization, P s p the spontaneous polarization, P pz the piezoelectric polarization, e 31 , e 33 , c 13 , and c 33 are the elastic and piezoelectric coefficients, relax is the parameter which controls the extent of piezoelectric polarization, a o is the strained lattice constant (Å), a is the unstrained lattice constant (Å). The fabricated bulk GaN devices show very little piezoelectric strain and hence a fully relaxed crystal is considered and only the spontaneous polarization is taken into account. Spontaneous polarization has an accepted value of around 0.29 C/m 268,69 which is incorporated in our model. The relax parameter is set to 1 in order to remove the effect of piezoelectric strain. Although the elastic constants and the piezoelectric coefficients are not used for the simulated devices their values [68][69][70] are listed in Table II in case they might be relevant for other specific device designs in the future.

D. Electron and Hole mobility
A number of Monte Carlo, 18-21,23-27,71 experimental 32-35 and modeling 71-75 studies have been performed to understand electron mobility in GaN. Unlike Si, GaN exhibits a negative differential mobility for electrons at very high fields due to the movement of the lower conduction band electrons to the higher M and L valleys, which have a higher electron effective mass, effectively reducing the mobility. So a dual mobility model approach is taken for GaN, one for the low field positive differential region and the other for the high field negative differential region. Proper selection of the peak electron drift velocity becomes important as it determines the transition between the two regions. Most Monte Carlo simulation overestimates the peak velocity against experimental data. According to Schwierz et al., 74 instead of the exact curve fitting of any set, it is better to find a reasonable compromise between the mathematical fit, the theory based expectations, and the tendencies observed for the experimental and theoretical data. The same methodology is followed in this work. Hole mobility, on the other hand has been much less explored and most work has been theoretical. [29][30][31]67,[76][77][78][79] Unlike electrons, the hole mobility does not show a negative differential mobility region and tends to saturate at high fields. At low field both electron and hole mobility is given according to the doping and temperature dependent Arora model according to Eq. (10): where µ min , µ ma x are the minimum and maximum mobility determined from Monte Carlo simulations or experimental data, α, β 1 , β 2 , β 3 , β 4 , are all temperature dependent fitting parameters, N r e f is the reference doping level, N the donor concentration and T is the temperature. The high field electron mobility is modeled according to Eq. (11) and implemented in Sentaurus by the Transferred Electron Effect 2 model. 71,73,74 This second model is triggered when the peak drift velocity is reached at the specified electrical field E c which takes into account the migration of the electrons from the lower Γ valley to the higher M-L valleys. The high field hole mobility is modeled according to the Caughey-Thomas relation 79 shown in Eq. (12) and implemented as the Canali model in Sentaurus TCAD.
where µ n0 , µ h0 are the previously determined value of the low field mobility µ 0(n/h) from equation (10) for electrons and holes respectively. E C , v n, sat , n 1 , n 2 , a are used in equation (11) and v h, sat and β in equation (12) are fitting parameters determined from Monte Carlo Simulations and E is the electrical field in both the cases.

E. Impact Ionization parameters
Defect induced high leakage current rather than avalanche breakdown is typically observed in GaN lateral HEMT power devices. However, avalanche breakdown is observed in vertical GaN PN diodes in recent experimental studies, 4-7,9 exhibiting a positive temperature coefficient of the breakdown voltage and hence proper modeling is vital for vertical power devices. The GaN impact ionization coefficient, which gives an estimate of the avalanche voltage has been studied both experimentally 36,37 and by Monte Carlo simulations. 24,25,27,[80][81][82] The avalanche multiplication is calculated by using the impact ionization coefficients according to equation (13). 83 where G ii is the generation rate of carriers due to impact ionization, α n/p are the impact ionization coefficients for electrons and holes and J n/p the electron and hole current. The impact ionization coefficient in GaN is calculated by the VanOverstraetendeMan model based on Chynoweth's equation 36 as shown in Eq. (14): where a n/p and b n/p the impact ionization rate parameters for electrons and holes and E is the electrical field. The selection of values for these parameters is explained in the next section.

F. Recombination-Regeneration
The three main processes for recombination/generation in bulk GaN are radiative, Shockley-Read-Hall (SRH), and Auger. The ABC model 84 used to study the carrier recombination processes in GaN in a shown in Eq. (15).
where J is the injected carrier density, d is the width of the quantum well, An is the SRH recombination rate, Bn 2 is the rate of radiative recombination and Cn 3 is the Auger recombination rate. The radiative recombination constant B (C rad for simulation purposes) is 1.1x10 −10 cm 3 s −1 for GaN. 44 SRH recombination is the dominant process of recombination for bulk GaN power devices. The SRH recombination-generation rate is given by Eq. (16). 53,85 where n and p are the electron and hole density, n i the intrinsic doping density, τ n and τ p the electron and hole lifetimes. Although SRH recombination is a trap assisted process in GaN, due to a lack of in-depth information on the trap level and density of the experimental data, in our simulations a constant electron lifetime of 0.7 ns and hole lifetimes of 2 ns are used. Auger recombination becomes the dominant method of recombination for high carrier densities which leads to non-radiative carrier loss and reduced high power efficiency. 53 Auger recombination in GaN occurs through phonon assistance unlike the direct Auger recombination observed in Silicon. 84 The Auger recombination rate is modeled according to Eq. (17). 53,85 where C p and C n are the Auger coefficients for holes and electrons, n and p are the electron and hole density and n i the intrinsic doping density. The value of the Auger coefficients for both holes and electrons used in this work is 3x10 −31 cm 6 -s −1 . 45,46 The models and parameters discussed in the above section are collated and the final values used in our simulations are presented in Table II. Tables I and II. The forward conduction and reverse breakdown simulations are performed for PN diodes with active area of 0.11 mm 2 (3.7 kV) and 0.72 mm 2 (2.6 kV) respectively. The 3.7 kV device has a 40 µm n-drift layer with a doping concentration of 6x10 15 cm −3 , a thin 0.5 µm n+ cathode region with a doping concentration of 1x10 19 cm −3 and a 1.5µm thick p+ layer with an effective doping of 5.2x10 17 cm −3 . The 2.6 kV device has a 15 µm thick n-drift layer with a doping concentration of 2x10 16 cm −3 , a 0.5 µm n+ cathode region with a doping concentration of 1x10 19 cm −3 and a 1.5 µm thick p+ layer with an effective doping of 5.2x10 17 cm −3 . The choice of a thicker device for forward simulations helps to minimize the contribution of parasitic and contact resistances when compared to the bulk resistance contribution from the thick n-drift layer. This leads to a more accurate modeling of the effect of mobility, lifetime and the bandgap on the I-V characteristics for the diode. However, detailed reverse bias breakdown voltage data for the 40 µm device was not reported which was available for the 15 µm device. Unfortunately, forward bias data for the 15 µm device wasn't reported. As a result, two different devices had to be simulated to compare both the reverse and the forward bias characteristics. The schematics of the two GaN P-N diodes designed are shown in Fig. 1.

A. Forward conduction
The forward bias simulations for the thicker device are compared with experimental data obtained from GaN P-N diodes 4 and is shown in Fig. 2. The p-ohmic contribution of 0.35 Ω as reported for this device is included in the simulations as a temperature dependent series resistance. The simulated data fits well with the experimental data for the intended current range. The forward knee voltage is around 3.1 V at 300 K. Very low current (below 1 µA.mm −2 ) is not fitted, as the physics pertaining to very low current needs further theoretical work in order to be properly implemented as physical models. The low current density (>10µA/mm 2 ) is shown in the inset of Fig. 2. Hole lifetime in the n-drift region has a major impact on the low current knee region and the electron lifetime does not affect the current significantly in this region. However, once the current is in full swing electron mobility dominates. Hole lifetime and hole mobility is assumed to have little impact on the I-V characteristics at high current densities owing to minimal conductivity modulation. This matches with the observations reported in Ref. 4.
The simulated forward I-V characteristics of the 3.7 kV GaN PN diode are compared with the experimental data 4 in Fig 3 for a junction temperature of 300-425 K with 25 K increments. The temperature dependent bandgap and incomplete ionization models discussed in the previous section are triggered for the different temperatures. The knee voltage decreases with increase of temperature and shows a reasonable fit with the experimental data over the entire temperature range. The slope of the current at higher densities decreases with increase in temperature pointing to a decrease in the mobility. We assume a T −3/2 temperature dependence of electron mobility in our simulations which results in good agreement with the measurement data. The effect of temperature on minority carrier lifetime is still not well understood, so we assume a temperature independent constant minority carrier lifetime in our simulation. The temperature dependence of electron mobility and bandgap are shown in Fig 4. The forward I-V characteristics of the GaN PN diodes are also simulated using different electron mobility model parameters reported by Farahmand et al., 71 Mnatsakanov et al., 73 and Schwierz et al. 74 as shown in Fig. 5. A similar model to one which we have used for electron mobility was used by all the three reports, the only change being the parametric values. It is observed in Fig. 4 that our mobility model parameters provide a much better fit with the measurement data than the other reported parameters. The electron mobility was underestimated in the literature, probably because they were based on the characterization of GaN materials with a higher defect density. Note that the defect density of the fabricated GaN PN diode from which experimental data is extracted shown in Fig. 4 is in a range of 10 4 -10 6 cm −2 , which is the lowest ever reported. simulations [19][20][21]23,27 and actual measurement. 33,34 The curves show two distinct regions across a peak drift velocity signifying the migration of electron into the higher M-L valleys for high electric fields. The low field drift velocity is similar for modeling and MC data except for the two experimental results from Wraback et al. 33 and Barker et al. 34 At higher fields however, the results seem to deviate from one another in assessing both the peak drift velocity as well as the value of electric field at which the migration starts to occur. Since there is a wide spread among the different works, instead of matching any particular curve we use a methodology which involved finding a median range. Although the simulated curve shows similar trend to the one reported by Kolnik et al. 19 at high fields, it is purely coincidental. A peak electron drift velocity value of 2.48x 10 7 cm.s −1 at the electric field of 1.76x10 5 V.cm −1 is assumed. The electron velocity decreases with increasing electric field after the peak velocity and finally saturates at 1.7x 10 7 cm.s −1 for very high fields. Due to the inadequacy of research for hole mobility, we can only compare the dependence of GaN hole drift velocity on electric field, using our selection of model parameter to theoretical Monte Carlo simulations 24,27,29,31 in Fig 7. Unlike electrons, there is no negative differential mobility region for holes and the hole drift velocity saturates with increasing electric field to a value of 5.6x10 6 cm.s −1 .

B. Reverse bias breakdown simulation
A comparison of the impact ionization coefficients from both Monte Carlo simulations 24,25 and experimental study 36 using the Chynoweth's equation for holes and electrons is shown in 7(a) and 7(b) respectively. It can be seen from Fig 8, at high electrical fields both the hole and the electron impact ionization parameters reported by Baliga et al. 36 is at least 5x-10x lower than the Monte Carlo simulated data. At low values of electric field, the hole impact ionization dominates while the electron impact ionization dominates at higher fields (>3MVcm −1 ) reported by Chen et al. 24 Baliga et al. 36 report that the dominating factor is hole impact ionization over the entire range of electrical field while Bellotti et al. 25 point out electron impact ionization as the major contributor to avalanche breakdown in a similar range considering both band to band tunneling and its absence at high electric fields. Hence, due to the reported discrepancy in the numbers as well as the trend, breakdown simulations are performed with all the reported impact ionization coefficients from Refs. 24, 25, and 36.  Reverse breakdown characterstics for the 2.6 kV device was simulated with the reported parameters and the resultant I-V characteristics is shown in Fig 9. The experimentally measured breakdown voltage for the 15µm diode with a n-drift doping of 2x10 16 cm −3 is 2600V. The breakdown voltage using parameters from Baliga et al. 36 yielded a value of 2605V while parameters from Chen et al. 24 give us a breakdown voltage of 2270V. Both tunneling and non-tunneling parameters from Bellotti et al. 25 were implemented and the resultant breakdown voltages were found to be around 1760V and 2190V respectively. Baliga et al. 36 reports that their parameters give a good fit for the breakdown voltage over a doping range of 1x10 14 cm −3 to 1x10 17 cm −3 and since it gives near exact fit with the experimental results, it is used in our simulation by itself.

IV. CONCLUSION
In this paper, physical models for energy bandgap, incomplete ionization, electron and hole mobility, polarization, impact ionization, radiative and non-radiative recombination and the associated model parameters in the literature are reviewed and compared. These mostly empirical model parameters have a wide spread and lead to TCAD simulation results with considerable disagreement among themselves and with the experimental results. We report a methodology of choosing physical models and model parameters based on the measurement data of the latest bulk GaN devices reported in the literature. Our TCAD simulation accurately models both forward and reverse I-V characteristics of high-voltage bulk GaN PN diodes, which are important power devices in their own right but also serve as a basic building block of other more complex GaN devices. It is worth noting that many second order effects such as the temperature dependence of electron and hole lifetimes, distribution of traps, etc. are not included in this work mainly because of the limited study on these subjects. Nevertheless, the good agreement between our modeling and experimental data provide us certain confidence that the physical models and parameters used in our work should provide a basis for TCAD simulation of other power devices fabricated in low dislocation density bulk GaN materials in the near future.