Monte Carlo Modeling of Wurtzite and 4 H Phase Semiconducting Materials

We present a discussion of the complexities encountered in particle simulation models 
for noncubic symmetry semiconductors, focusing on the wurtzite and 4H polytypes 
of GaN and SiC. We have identified three general issues, band structure, scattering 
mechanisms, and band intersections, which in our opinion, constitute the most important 
modifications to conventional Monte Carlo simulators for cubic symmetry 
semiconductors. It is found that the band intersection points present the greatest 
modeling challenge. We discuss the effect of band intersections on the transport 
dynamics and our initial attempts at treating transport near these points. Comparison to 
experimental data is also made.


INTRODUCTION
Owing to their intrinsically larger band gap than conventional silicon and GaAs semiconductors, wide band gap semiconductors are better suited to many emerging device applications such as high radiation environments typically encountered in space flight, high power, high frequency amplifiers used in base stations for wireless telecommunica- tions networks, and high temperature conditions found in automotive and jet engines.The most important wide band gap semiconductors are the III-nitrides, GaN, InN, A1N [1] and their related ternary compounds and SiC [2].In each of these materials, the principal energy gap is about two to three times larger than that of Si or GaAs.As a result, their breakdown electric field strengths, immunity to radiation, and power handling capa- bility are expected to be significantly greater than Si or GaAs [3][4][5][6].
Though these materials are of great interest, their technological maturation is frustrated to a large extent by the lack of reproducible high quality material.As a result, experimental char- acterization of many of the intrinsic transport properties of these materials has been limited.The Corresponding author.Tel.: 404-894-6767, Fax: 404-894-0222, e-mail: kbrennan@ece.gatech.edulack of experimental data generally negatively im- pacts device simulation capability.In most macro- scopic simulation tools, the critical modeling parameters, such as carrier mobilities, lifetimes, velocity-field characteristics and ionization coeffi- cients, are obtained empirically.In order to pro- ceed in evaluating the ultimate performance levels of the wide band gap semiconductors prior to the technical maturation of these materials, it is clearly necessary to utilize a theoretical technique that can proceed relatively independently of experiment.To this end, we have developed a general methodo- logy, materials theory based modeling, that is based on fundamental materials theory results rather than on empirically determined materials characteristics.
It is the purpose of this paper to provide an overview of some of the important modifications that need be made to standard ensemble Monte Carlo simulations in order to properly examine transport in noncubic symmetry semiconductors.
In the next section, we discuss these modifications.In section three we present calculations of the impact ionization coefficients in 4H-SiC that illustrate the importance of these modifications in gaining agreement with experimental measure- ments.Finally, conclusions and suggestions for future work are presented in the conclusions.

MODEL DESCRIPTION
The details of our approach have been exhaus- tively reported elsewhere [7,8].Nevertheless, the salient features of the technique can be summar- ized as follows.The transport model utilizes a full band ensemble Monte Carlo simulation.Owing to its inherent flexibility and rigor, the Monte Carlo technique provides a detailed solution of the Boltzmann equation limited only by the descrip- tion of the underlying physics included within the model.Properly tailored Monte Carlo models can thus provide a highly reliable and relatively complete description of the transport dynamics of carriers within a material or device.For these reasons, we utilize the Monte Carlo model in the calculations presented herein.
The two major ingredients within the Monte Carlo code are the band structure and carrier scattering rates.Though many different approxi- mations are often employed to simplify the band structure and scattering rate descriptions used within Monte Carlo models, the most complete description that insures an accurate accounting of high field transport requires detailed numerical formulations of these features [9,10].Band

Structure
There exist numerous techniques for calculating the electronic band structure of semiconductors.These techniques can be categorized into two general strategies, ab initio or semi-empirical.The ab-initio techniques, as the name implies, are first principle approaches in which the band structure is determined directly from the numerical solution of the Schroedinger equation for the system.Unfortunately, the general problem is a many body problem that cannot be solved in closed form.As a result, several simplifications must be introduced.The most common approach is to simplify the many body problem using density functional theory [11].Using the linear augmented plane wave, LAPW, method [12] within the density functional theory, the details of the band structure can be calculated.However, these approaches are highly computationally expensive.As a result, usage of ab-initio techniques is currently prohibi- tive for scattering rate and impact ionization transition rate calculations wherein numerical expressions for the wavefunctions, sampled over millions of k-points, are required.For this reason, we utilize the empirical pseudopotential method, which is a more computationally efficient algorithm for generating the band structure data.
We have developed a genetic algorithm for optimizing the pseudopotential band structure calculation [13,14].The effective potentials used in the pseudopotential calculation are optimized through an iterative scheme in which the band structures are calculated recursively and selected features are compared to experimental and/or ab initio results.The resulting band structures are in excellent overall agreement with experimental/ab initio values for the energy gaps at high symmetry points, band ordering, and effective masses.In addition, the pseudopotential band structure cal- culation is highly computationally efficient.
Independent of the procedure adopted for calculating the band structure, the determination of the band structure for wurtzite and 4H polytypes is more challenging than that for cubic symmetry polytypes.Additional complexity arises in the wurtzite and 4H polytypes since the Brillouin zone is smaller and the number of atoms per unit cell is larger.As a result, more bands appear and the overall complexity of the band structure increases.This can be vividly seen from inspection of Figure 1, wherein the band structures of zincblende (ZB) and wurtzite (WZ) phase GaN and 4H-SiC are shown.As can be seen from the figures, as the materials progress from ZB GaN to 4H-SiC, the overall complexity of the band structure markedly increases.

Scattering Rates
At high electric field strengths, phonon scattering and impact ionization dominate the transport physics in semiconductors.In this paper, we restrict ourselves to discussion of these mechanisms.From its spherical symmetry properties, it can be shown [15] that the only allowed electron- optical phonon interaction within the gamma valley in the ZB phase is that involving the LO mode.In noncubic semiconductors, such as wurtzite and 4H, due to the different symmetry properties of the crystal, both the TO and LO optical phonon modes can contribute to electron- optical phonon scattering.Lee et al. [16] have recently examined the polar optical phonon scattering rate in wurtzite crystals.They showed that in general the matrix element governing the interaction between electrons and TO-like optical phonon modes is not zero since the TO-like mode is not purely transverse in wurtzite.As such the TO branch does interact with electrons in wurtzite polytypes unlike the situation in cubic semicon- ductors.Though the TO-like mode does contri- bute to the optical phonon scattering rate, its contribution in wurtzite phase GaN, as shown by Lee et al. [16], is very weak.Therefore, for bulk wurtzite phase GaN, the cubic interaction Hamiltonian well approximates the situation.We have used this approximation in our formulation of the polar optical phonon scattering rate in our Monte Carlo calculations.
For acoustic phonon scattering the different symmetries of the WZ and ZB phases result in different forms of the deformation potential tensor.In ZB phase material, the deformation potential in the central valley of the first conduc- tion band is a second rank tensor with diagonal form and equal diagonal elements [17].The deformation potential retains its diagonal form in the WZ phase, but due to its lower symmetry, the diagonal elements are not necessarily all equal.
The value of Ezz is generally expected to be different from E= Eyy.Though some informa- tion exists about the valence band deformation potentials [18] for wurtzite GaN, to our knowl- edge, there is no information available at room temperature for the conduction band.For this reason, and because the effective mass along the hexagonal axis is within 1% of that in the basal plane, in our calculations it is assumed that all of the diagonal elements of the deformation potential tensor are the same in WZ GaN.The parameters used for calculating the scattering rates for both the conduction and valence bands in WZ GaN have been collected elsewhere [8,19].In all of the calculations presented, the scattering rates are determined numerically by integrating the matrix elements for the transition rates over the full Brillouin zone.As a result, any anisotropy in the rates [20] is naturally included in our calculations.
For high field transport analysis, it is important to include the details of the impact ionization transition rates.In our simulation, the impact ionization transition rates are calculated numeri- cally [21,22].The calculation is performed by integrating Fermi's golden rule for a two-body, statically screened Coulomb interaction over the possible final states using a numerically generated dielectric function [23] and pseudowavefunctions.
Band Intersections In our opinion and experience, the most pressing problem introduced in simulating non-cubic sym- metry semiconductors is that of band intersections.
As discussed above, owing to the smaller Brillouin zone and larger number of atoms per unit cell, the band structure of these materials is highly complex leading to numerous band intersection points.
These points complicate the Monte Carlo simula- tion since a decision must be made as to which band the carrier will reemerge upon after passing through a band intersection region.As we will show, consideration of the transport dynamics at these points is critical to accurately predict the transport parameters in these materials.Standard Monte Carlo algorithms provide no means of treating these occurrences since a carrier is assumed to lie only within a single band that does not change while drifting within an applied field.
We have developed a new technique for deter- mining the trajectory of a carrier in the presence of band intersections.In our bulk Monte Carlo simulator, the trajectory of the carrier following a band intersection is determined stochastically following the usual Monte Carlo procedure, by evaluating the magnitude of the overlap integral between the cell periodic part of the initial state and cell periodic part of each of the possible final states [24,25].As a result, when drifting through a band intersection region, a carrier can potentially transfer from one band to another.
For lack of a better term, we refer to this process as interband tunneling.
Though we have utilized the overlap integral test for treating band intersections in bulk material, we have found that it is presently too computationally expensive for use in device simulations.As an alternative, we use a velocity continuation test at band intersections.As in the case for the overlap integral test, the final state is selected stochasti- cally, but with the difference that the choice is based on the requirement of the carrier velocity continuation during the drift.The final state is selected from the possible alternative paths weighted by the closeness of the corresponding final group velocities to that of the initial state.
Finally, it should be noted that neither approach can be considered to definitively describe the physics of band intersections.Instead, these two approaches should be viewed as first attempts to include the complexities of multiband transport in ensemble Monte Carlo simulations.We are presently examining more complete formulations of this problem.

TRANSPORT CALCULATIONS
The modification of the ensemble Monte Carlo code for noncubic symmetry materials that has the greatest influence on the calculated results is the band intersection treatment.The effect of band intersections on the transport dynamics is most vividly illustrated by examining high field trans- port, particularly impact ionization.As is well known, the impact ionization coefficients are very sensitive to the details of the high energy tail of the energy distribution function.Comparison of the calculated ionization coefficients both with and without the band intersection treatment described above clearly demonstrates the importance of this effect as will be shown below.
As discussed above, little experimental infor- mation is presently available for the transport parameters for most of the wide band gap semiconductors.This is especially true for the im- pact ionization coefficients.However, there are some direct experimental measurements for the hole initiated impact ionization coefficient in bulk 4H-SiC [26,27].The Monte Carlo calculated results with and without tunneling at the band intersection points are plotted along with both sets of experi- mental measurements in Figure 2. As can be seen from Figure 2, for fields applied in directions within the basal plane (1000), the ionization coefficients are calculated to be nearly the same when tunneling is included or not in the model.The fact that the two models give about the same result stems from the fact that the carriers have many ways of drifting to high energy without passing through band inter- section points.However, when the field is applied along the c-axis, (0001), there is a significant difference between the calculated coefficients.Notice that the calculated coefficient when tunneling is not present (marked by circles in Fig. 2) is very much less than that when tunneling is present (marked by diamonds in Fig. 2).Inverse Electric Field (10 .7 cmN) FIGURE 2 Plot of the calculated and experimental measurements of the hole initiated impact ionization coefficients in bulk 4H-SiC for an applied electric field along the c-axis.Notice that good agreement between the calculated and experimental measurements of Ragunathan and Baliga occurs when tunneling effects at mixing and crossing points are included in the calculations.If tunneling effects at band intersection points are neglected, the otherwise identical Monte Carlo calculations deviate by about two orders of magnitude from the measured experimental results.Notice that there is relatively no dependence on the calculated ionization coefficients with tunneling effects at band intersection points for fields applied in the 1000 direction in the basal plane.
notice that the calculations including interband tunneling give much closer agreement to either set of experimental measurements.The calculations without tunneling are orders of magnitude lower than either set of experimental data.The much lower ionization coefficients calcu- lated in the absence of tunneling arise from the fact that the hole distribution is much colder than when tunneling is included.This can be under- stood as follows.The bands are labeled in ascending order at each discrete k point, independ- ent of the band's apparent continuity.Therefore, in the absence of interband tunneling, when a hole encounters a band intersection region, it continues along the lowest energy band by default and cannot drift into the higher energy band, even if the continuity of the derivative implies that it should.Consequently, the hole can only reach higher energy through scattering events and not through drift.Since scatterings are relatively more infrequent, this results in a colder hole distribution and a lower impact ionization coefficient when tunneling is not included.

CONCLUSIONS
New complexities arise in Monte Carlo based modeling of noncubic symmetry materials that require some modifications of the standard codes.
We have identified three general areas where important revisions must be made.These are the band structure, scattering rates and band intersec- tions.It is our opinion that the band intersections present the greatest challenge and have the greatest effect on the calculations.
We have shown that the inclusion of the possibility of interband tunneling in the vicinity of band intersection points is necessary in order to enable significant carrier heating in materials which have many band intersection points.Calcu- lations of the hole initiated impact ionization coefficient, which reflects the high energy dynamics of holes, show excellent agreement with experi- ment only when tunneling is included in the model.Though treatment of the band intersection points is crucial to explaining the transport physics of these materials, it should be noted that our model is at best preliminary.A more complete formula- tion of the band intersection physicslis presently underway.
FIGURESketch of the calculated energy band structure of the a) zincblende phase of GaN b) the wurtzite phase of GaN and c) the 4H phase of SiC.All of the band structures have been calculated using the empirical pseudopotential method. FIGURE(Continued).