# A Simulator for Investigation of Breakdown Characteristics of SiC MOSFETs
## 2 Simulation Method and Validation  
### Simulation Method  
The framework of the simulator we developed is shown in Figure 1. It contains four main modules, including the input, material database, solver, and outputs. The device structure and doping profile are provided, as well as the solving command which is parsed by a self-developed input parser through input module, combined with the relevant material and physical model parameters from the material database loaded into the solver, which is the key module. In the solver module, the basic steps are to discretize the equations to obtain a system of nonlinear equations and then select an appropriate numerical method to solve them. It is worth noting that the independent variables of the basic equations vary greatly in numerical magnitude and show different orders of magnitude and have strongly different behaviors in regions with varying space charge densities. For drift-diffusion (DD), the independent variables are \(\psi\), n and p, while for hydrodynamic (HD), they are \(\psi\), n, p, Tl, Tn and Tp, where \(\psi\) is the electrostatic potential, n is the electron density, p is the hole density, Tl is the lattice temperature, Tn is the electron temperature and Tp is the hole temperature. The variables in these equations must be appropriately scaled before discretizing the basic equations to prevent the fluctuations of the potential from being masked by the larger carrier densities. Figure 1: Framework of the simulator. Here, \(\rm{k_{B}T/q}\), \(\sqrt{N_{c}N_{v}}\) and \(1/300\) are used to scale the electrostatic potential, carrier densities and temperature, respectively, where \(\rm{k_{B}}\) represents the Boltzmann constant, \(\rm{q}\) represents the elementary charge, \(\rm{N_{c}}\) represents the density of state (DOS) at conduction band and \(\rm{N_{v}}\) represents the DOS at the valance band. In order to avoid using the extremely low intrinsic carrier concentration caused by the wide bandgap directly, we choose the ln value of the scaled intrinsic carrier concentration to calculate the initial guess of the carrier densities. The intrinsic carrier concentration can be expressed as  
\[n_{i}=\sqrt{N_{c}N_{v}}\cdot\exp(-\frac{E_{g}}{2k_{B}T}) \tag{1}\]  
and the scaled intrinsic carrier concentration with the ln value is given by  
\[\ln(n_{i})_{scale}=\ln[(N_{c}N_{v})^{0.5}]-\frac{E_{g}}{2k_{B}T} \tag{2}\]  
which gives the safer expression for scaled electron density:  
\[n_{scale}=\exp[\psi_{scale}-\psi_{n,\ scale}+\ln(n_{i})_{scale}] \tag{3}\]  
where \(\psi_{scale}\) represents the scaled electrostatic potential and \(\Phi_{n,scale}\) is the quasi-Fermi level for electrons. This numerical treatment method can effectively prevent the occurrence of underflow, and thus avoid convergence issues and improve the stability of the numerical simulation. Then, the box discretization method, where an integral form is used to describe the basic laws of physics, and this description is later applied to each subdomain throughout the mesh in order to discretize the Poisson equation, carrier continuity equations and the energy conservation equations [13; 14; 15; 16]. Additionally, the center difference method is not utilized in discretizing the transport and the energy flux equations. This is because the simple discretization format of the center difference method is highly restrictive and can only maintain numerical stability if the drift current induced by the electric field is significantly smaller than the diffusion current. Therefore, to ensure the numerical stability in the presence of strong drift currents, the Scharfetter-Gummel method [17; 18; 19] is employed. After discretization, the Poisson equation coupled with the current continuity equations based on the DD or the HD transport model are solved self-consistently. The Newton method is utilized to solve the nonlinear equations due to its excellent stability and rapid convergence speed.