Monte Carlo simulations of electron transport in 4H-SiC using the DFT-calculated density of states

We carry out Monte Carlo simulations of electron transport in 4H-silicon carbide (4H-SiC) based on the numerically calculated density of states (DOS) to obtain the electron mobility at low electric fields. From the results, it can be concluded that a correct calculation of the DOS requires a very dense wavevector k-mesh when low electron kinetic energies are considered. The crucial issue is the numerical efficiency of the DOS calculation. We investigate the scaling efficiency when different numbers of cores are used.


Introduction
Modern electronics requires simulation tools that can combine atomic-scale simulations and the continuous domain of semiconductor devices. We focus herein on 4H-SiC, a wide-bandgap semiconductor used to design power devices with increased efficiency [4]. The Monte Carlo method [13] takes into account scattering events and a band-structure model to investigate electron transport. The macroscopic mobility is then obtained by averaging the drift velocity with respect to an external electric field. As described in Ref. [13], the simplest approach considers analytical formulas for the energy band structure and the density of states. However, this approach is correct only for very low energies and ideal crystalline structures. To model the transport of electrons with even moderate energies above ca. 0.15 eV in 4H-SiC, one cannot use analytical equations for the DOS and the energy band-structure model. In Ref. [1], the authors suggested using simple equations for the scattering rates and to combine these formulas with the DOS calculated numerically, using the density functional theory (DFT) method. The DOS was calculated therein for surprisingly high energies (up to 10 eV), but no consideration was given to the low energy range. We propose to extend this approach and calculate the DOS very accurately for energies between 0 and 0.25 eV. We also suggest combining the Brooks-Herring (BH) and Conwell-Weisskopf (CW) models for scattering on ionized impurities to keep the Monte Carlo approach valid for the whole range of doping concentrations.

Monte Carlo method
In the Monte Carlo method, an electron is considered to be a point charge with an effective mass. The free flight of electrons is randomly interrupted by one of the scattering events considered during the simulation. The scattering phenomena are therefore an essential element of the Monte Carlo model. For the presented simulations of 4H-SiC, the following scattering models are included: -Acoustic phonon scattering: where T is the temperature, c L = 2 S with s is the speed of sound in 4H-SiC, is the mass density, d is the acoustic phonon deformation potential, and E is the electron energy -Ionized impurity scattering: -Brooks-Herring (BH) model [13] where , N I is the density of ionized impurities, S is the dielectric permittivity of 4H-SiC, Z = 1 , e is the elementary charge, kis the wavevector, and D is Debye length, combined with -Conwell-Weisskopf (CW) model where m d is the effective density of states and is the nonparabolicity parameter of the band shape.
-Polar optical phonon scattering [1] for absorption and emission: where ℏ op is the polar optical phonon energy, N ℏ op is the phonon number calculated from the Bose-Einstein distribution, and where ℏ H is the nonpolar optical phonon energy, Z H = 3 for 4H-SiC, N ℏ H is the phonon number calculated from the Bose-Einstein distribution, and D H is the polar optical deformation potential. The material parameters are summarized in Table 1.

Model selection
The DOS for 4H-SiC is calculated using MedeA VASP software [6,8]. For the atomistic simulation, it is assumed that the elementary crystal cell has the following dimensions: a = b = 3.073 Å and c = 10.053 Å. To obtain the DOS, two approaches are compared. The first is the DFT method with generalized gradient approximation (GGA) PBEsol functionals [10]. The second algorithm uses MJLDA meta-GGA functionals [14]. For both methods, the plane-wave cutoff energy is set to 400 eV. The basic difference between these two approaches is that the first one, DFT with GGA, underestimates the energy gap E g for semiconductors [12]. The DFT with GGA returned E g = 1.99 eV, while for MJLDA with meta-GGA, we obtained E g = 3.17 eV, which is very close to the value reported in literature E g = 3.26 eV [3]. 19 Polar optical phonon energy, ℏ op (eV) 0.12 Low-frequency dielectric constant, S (1) 9.7 High-frequency dielectric constant, 0 (1) 6.5 Nonpolar optical phonon energy, ℏ H (eV ) 0.0854 Polar optical deformation potential, ℏ H (eV/m ) 9 × 10 10 The most critical parameter for the DOS calculation is the k-mesh spacing of the Brillouin zone. In Ref. [1], the DOS was calculated for a wide range of energy. In our opinion, the range does not have to be so large; however, the more accurate calculation should be done for the low-kinetic-energy region. This is mainly important for the range of energies where a linear dependence between the drift velocity and electric field is expected. From our previous simulations, we can expect that the energy of electrons would not exceed 0.15 eV [15]. Several atomistic simulations are performed to find an appropriate k mesh for the DOS calculation. Exemplary results are shown in Fig. 1. The DOS plots are shifted so that the zero energy on the x-axis corresponds to the bottom of the valence band. The values obtained from the analytical formula in Eq. (7) [13] for the DOS are also plotted to compare the results. The mesh parameters are summarized in Table 2. where m d = m 1 m 2 m 3 1∕3 , = 0.323 [9], m 1 = 0.58 , m 2 = 0.28 , m 3 = 0.31 of the electron mass m 0 [7].
From the presented calculations, one can see that a very dense k-mesh with 89 × 89 × 25 points is required to obtain a correct value of the DOS near the edge of the conduction band. One should also note that the analytical formula for the DOS in Eq. (7), which takes into account a single valley, is not valid above 0.15 eV. As shown in Fig. 2, although the correct value of E g can be obtained from the MJLDA meta-GGA functional, it overestimates the value of the DOS.

The effect of variation in the DOS
The scattering rate is an increasing function of the DOS, as shown by Eqs. (1)(2)(3)(4)(5)(6). Underestimating the DOS leads to lower scattering rates and higher values of the calculated mobility. The total scattering rate ( Tot ), i.e., the sum of scattering by all mechanisms calculated for different sets of DOS, is shown in Fig. 3. The MJLDA meta-GGA model overestimates the scattering rates for energies above ca. 0.15 eV, but underestimates the scattering for low energies. The accuracy of the electron mobility thus depends also on the external electric field applied. A histogram of the electron energies for an external field E = 10 6 V/m is shown in Fig. 3 (right y-axis). There are almost no particles with energies above 0.15 eV. For this value of the electrical field, the velocity and external field are proportional. For this range of energies, both the MJLDA meta-GGA model and the DFT with GGA-PBEsol with a coarse k-mesh underestimate the DOS. This leads to erroneous, too low values of the scattering rates.

The scalability of the calculations
All the calculations are performed on a Linux system with 128 GB RAM, 2× Intel Xeon CPU E5-2680 v3 (12 cores, corresponding to 24 logical cores each) using MedeA VASP 5.3 software. VASP 5.3, as described in the documentation [8], is compiled using the Intel Fortran compiler (IFORT v.14) and uses the Intel Math Kernel Library. The most crucial parameter is the time of computing (ToC). The VASP calculations can benefit from parallel computing. The ToC is shown as a function of the number of logical cores used during the calculations in Fig. 4. The ToC is reported for a single core and for even numbers of cores from 2 to 36. Missing data mean that the calculation process failed and no DOS was returned.
The ToC to calculate the DOS by means of MJLDA meta-GGA is more than three times longer than when using the DFT GGA-PBEsol method. However, the optimal number of parallel cores is ca. 4-8 for the calculation for the eight-atom elementary cell of 4H-SiC is similar when using both approaches. Using more than eight cores does not affect the ToC. Moreover, for the MJLDA meta-GGA model, the ToC increases when more than 24 cores are used. The synchronization between so many threads delays the execution of the routines.
The DFT GGA-PBEsol calculation failed when more than 16 threads were requested. Since only logs can be accessed, one can conclude that the problem occurred inside one of the Intel MKL routines due to incorrect communication between different threads.

The scattering, Brooks-Herring, and Conwell-Weisskopf models applied
The essential quantities for the MC simulations are the scattering rates. The calculated scattering rates are shown in Fig. 5. Two models (BH and CW) are used for the ionized doping scattering rates. The scattering rates for the BH and CW models are shown in Fig. 6.
The Monte Carlo simulations switch from BH to CW when the following condition is fulfilled:

Incomplete ionization
Note that the results shown above assume full ionization of dopants. However, this assumption is not valid for widebandgap semiconductors. Consider nitrogen, the most popular n-type dopant for 4H-SiC. A nitrogen atom may be inserted at two different locations, which results in two different dopant ionization energies: E Dh = 0.042 eV and E Dc = 0.083 eV [2]. The total ionization rate I for a donor atom density of N D is with [16] where N C is the effective density of states for electrons and g D = 2 is the degeneracy factor [16]. Assuming that N Dh = N Dc = 0.5N D and an electron density of n = N + D , the total ionized density can be expressed as Equation (12) is a third-order polynomial in N + D . Figure 7 shows the ionization rate of nitrogen dopants at 300 K, 500 K, and 700 K. At 300 K, the ionization level is lower than 50% for a doping concentration of 5 × 10 18 cm −3 .

Electron mobility
The mobility values calculated using the full and incomplete ionization models are shown in Fig. 8. In both cases, the combination of the CW and BH ionized dopant scattering models is used, as described in Sect. 4.1. Using only the BH approach, which is appropriate for low and intermediate doping concentrations, leads to an incorrect increase of the mobility for high doping concentration levels. We show this effect in Ref. [5]. The mobility for the incomplete ionization Fig. 5 The scattering rates for a doping concentration of N D = 10 17 cm −3

Fig. 6
The BH and CW scattering rates as a function of the donor concentration 1 3 model is slightly higher. This is caused by the fact that the ionization scattering rate is less significant due to the smaller number of ionized atoms, and nonionized impurity scattering is not considered in these simulations.

Conclusions
It is shown that calculations of the DOS in the range up to 0.25 eV must be carried out using a dense mesh in the k-domain. We recommend using more than 50 k-points in directions parallel to the c-axis and 25 in the c-axis. The DFT GGA PBEsol approach is more efficient than the MJLDA meta-GGA model and also calculates the DOS more accurately. However, to properly report the E g value, the MJLDA meta-GGA method must be used first to correct the DFT GGA PBEsol calculations of E g . The optimal number of cores is from four to eight when the elementary cell of 4H-SiC is considered. One should not use more than 16 cores since this may lead to failure of the calculations.
With the appropriate use of BH and CW, it is possible to obtain a correct value for the mobility over the entire range of doping concentrations. However, it is also vital to include the incomplete ionization model. The routine proposed in this article can be used for any semiconductor, including also crystals with defects. The current study demonstrates that utilizing DFT simulations and the Monte Carlo method enables the creation of a link between the world of atoms and the continuous domain of solid-state electronics.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflicts of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.

Fig. 7
The incomplete ionization rate for nitrogen dopants in 4H-SiC Fig. 8 A comparison of the mobility values calculated using the full and incomplete ionization model. The measurement data are from Ref. [11]