A Simulator for Investigation of Breakdown Characteristics of SiC MOSFETs

: Breakdown characteristics play an important role in silicon carbide (SiC) power devices; however, the wide bandgap of SiC poses a challenge for numerical simulation of breakdown characteristics. In this work, a self-developed simulator employing a novel numerical processing method to prevent convergence issues, based on semi-classical transport models and including several kinds of mobility, generation and recombination models, is used to investigate the performance and break-down characteristics of 4H-SiC MOSFETs in high-power applications. Good agreement between our simulator and an experiment and commercial TCAD was achieved. The simulator has good stability and convergence and can be used as a powerful tool to design and optimize semiconductor devices. Further, the breakdown characteristics are evaluated with different factors, including lattice temperature, device structure and doping profiles. Our results show that the doping profile plays the most important role in the breakdown voltage, followed by the device structure, while the impact of lattice temperature is found to be minimal.


Introduction
Silicon carbide (SiC) has emerged as an attractive candidate for high power-density switching applications due to its superior material properties.SiC has three times the bandgap width, ten times the critical electric field, three times the thermal conductivity and nearly twice the saturation velocity of silicon (Si).In addition, the three-times bandgap width of SiC results in an intrinsic carrier concentration 18 orders of magnitude lower than that of Si [1][2][3][4].These attributes demonstrate significant potential for developing power devices capable of operating at higher power levels, elevated temperatures and increased frequencies, while maintaining reduced leakage currents.Among the various polytypes, the 4H-SiC crystal structure stands out as a commercially accessible option due to its wider bandgap, higher carrier mobility compared to other polytypes, exhibiting superior electrical properties [5,6].
SiC-based power metal oxide-semiconductor field-effect transistors (MOSFETs) are favored for their low gate input impedance, fast switching speed and reduced on-resistance relative to other power devices [7,8].The high electric breakdown field of SiC enables these MOSFETs to switch power much faster than comparable silicon devices, resulting in lower energy losses.Consequently, breakdown voltage (BV) is considered a critical parameter for power devices [9,10].In order to accurately and comprehensively evaluate the breakdown voltage of 4H-SiC devices and further design and optimize device performance, accuracy and efficiency simulating the I-V characteristics and breakdown characteristics are essential.
However, the wide bandgap of 4H-SiC results in an extremely low intrinsic carrier density of about 5 × 10 −9 cm −3 at room temperature [11], compared to 1 × 10 10 cm −3 for Si.This low intrinsic carrier concentration causes underflow, convergence issues and sharp carrier density variations.These make the simulation of 4H-SiC devices, especially breakdown voltage, challenging for numerical simulations.While using higher arithmetic precision and adaptive mesh dissection strategies [12] could improve the convergence problem, it is still worth exploiting to achieve better numerical stability.Moreover, to achieve more accurate device simulation, modifications to existing models or the development of new ones are often required.Independent and flexible device simulators are needed to incorporate these models.Therefore, this work demonstrates a robust simulator for 4H-SiC with high numerical accuracy and convergence, to provide an accurate analysis of the breakdown characteristics and device performance.
The main aim of this paper is to employ a self-developed simulator to investigate the breakdown characteristics of 4H-SiC MOSFETs.The remaining contents of the paper are arranged in the following manner.Section 2 outlines the simulator's framework, explicates the numerical methods for scaling independent variables in fundamental equations to improve numerical stability, discusses the discrete methods used in the simulator, and presents the physical models included.For validation, comparisons are made with both experimental results and commercial TCAD.Section 3 discusses the breakdown characteristics' dependence on the lattice temperature, device geometry and doping profile.Finally, conclusions are drawn in Section 4.

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 ψ, n and p, while for hydrodynamic (HD), they are ψ, n, p, Tl, Tn and Tp, where ψ 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.
Appl.Sci.2024, 14, x FOR PEER REVIEW 3 of 13 which gives the safer expression for scaled electron density: where ψscale represents the scaled electrostatic potential and Ф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.In each step of the Newton iteration, manual differentiation is used to compute the elements of the Jacobian matrix due to its fast speed.Then, the solving process is performed using the Portable, Extensible Toolkit for Scientific Computation (PETSc) library [20].
Various models, including mobility models, such as bulk mobility, doping and temperature-dependent mobility (the Masetti mobility model), surface degradation mobility and high-field saturation mobility model [21][22][23][24][25][26][27], Auger recombination [28], Shockley-Read-Hall recombination [29], impact ionization [30], incomplete ionization [31] and bandgap narrowing [32], are incorporated into the simulator to provide a comprehensive understanding of the device's behavior and to ensure that the relevant factors have been taken into account.This has enabled us to simulate the electrical and breakdown characteristics more accuracy.The incorporation of these diverse models into the simulation framework is particularly instrumental in accurately characterizing the breakdown voltage.
If a high electric field is present alongside a long space charge region, the carrier acceleration length will be larger than the mean free path between two impact ionizations, Here, k B T/q, √ N c N v and 1/300 are used to scale the electrostatic potential, carrier densities and temperature, respectively, where k B represents the Boltzmann constant, q represents the elementary charge, N c represents the density of state (DOS) at conduction band and 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 ) and the scaled intrinsic carrier concentration with the ln value is given by which gives the safer expression for scaled electron density: where ψ scale represents the scaled electrostatic potential and Φ 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.In each step of the Newton iteration, manual differentiation is used to compute the elements of the Jacobian matrix due to its fast speed.Then, the solving process is performed using the Portable, Extensible Toolkit for Scientific Computation (PETSc) library [20].
Various models, including mobility models, such as bulk mobility, doping and temperaturedependent mobility (the Masetti mobility model), surface degradation mobility and highfield saturation mobility model [21][22][23][24][25][26][27], Auger recombination [28], Shockley-Read-Hall recombination [29], impact ionization [30], incomplete ionization [31] and bandgap narrowing [32], are incorporated into the simulator to provide a comprehensive understanding of the device's behavior and to ensure that the relevant factors have been taken into account.This has enabled us to simulate the electrical and breakdown characteristics more accuracy.
The incorporation of these diverse models into the simulation framework is particularly instrumental in accurately characterizing the breakdown voltage.
If a high electric field is present alongside a long space charge region, the carrier acceleration length will be larger than the mean free path between two impact ionizations, resulting in carrier multiplication.The ionization coefficient, α, is defined to represent the reciprocal of the mean free path.Therefore, the generation rate can be expressed as where J n and J p are the electron and hole current density as in the carrier transport model and based on the van Overstraeten-de Man model, the ionization coefficient α is with where T is the lattice temperature, T 0 the fixed temperature of 300 K, E the parallel electric field and the factor γ is expressed by the optical phonon energy hω op .

Validation
The vertical diffused SiC MOSFET (VDMOSFET), shown in Figure 2, has excellent performance and is widely used.To verify and calibrate the simulator, we compared the output characteristics of VDMOSFET with those obtained from the experimental measurement [33].The structure parameters are listed in Table 1, where * represents parameters for verifying the experimental result and # represents parameters for other simulations.It should be noticed that the doping profile within the base region is assumed to be a Gaussian distribution with a peak doping concentration of 5 × 10 16 cm 3 .
where Jn and Jp are the electron and hole current density as in the carrier transport model and based on the van Overstraeten-de Man model, the ionization coefficient α is with where T is the lattice temperature, T0 the fixed temperature of 300 K, E the parallel electric field and the factor γ is expressed by the optical phonon energy ℏωop.

Validation
The vertical diffused SiC MOSFET (VDMOSFET), shown in Figure 2, has excellent performance and is widely used.To verify and calibrate the simulator, we compared the output characteristics of VDMOSFET with those obtained from the experimental measurement [33].The structure parameters are listed in Table 1, where * represents parameters for verifying the experimental result and # represents parameters for other simulations.It should be noticed that the doping profile within the base region is assumed to be a Gaussian distribution with a peak doping concentration of 5 × 10 16 cm 3 .
While performing the simulations, the Masetti mobility model, surface degradation mobility and high-field saturation mobility model, Auger recombination, impact ionization incomplete ionization and bandgap narrowing are used.
Figure 3 presents comparisons for the VDMOSFET, with the output characteristics at T = 300 K for Vgs changing from 7 V to 10 V, and the transfer characteristics at a Vds of 10 V.As shown in Figure 3a, when Vgs is low, good agreement is achieved, while Vgs = 10 V, the ON-current are consistent with the experimental result.The difference in the saturation region is due to the limitations of existing physical models, especially high-field transport models.Figure 3b shows that the results achieved good agreement.The model parameters used in the simulator are listed in Table 2.While performing the simulations, the Masetti mobility model, surface degradation mobility and high-field saturation mobility model, Auger recombination, impact ionization incomplete ionization and bandgap narrowing are used.
Figure 3 presents comparisons for the VDMOSFET, with the output characteristics at T = 300 K for Vgs changing from 7 V to 10 V, and the transfer characteristics at a Vds of 10 V.As shown in Figure 3a, when Vgs is low, good agreement is achieved, while Vgs = 10 V, the ON-current are consistent with the experimental result.The difference in the saturation region is due to the limitations of existing physical models, especially high-field transport Appl.Sci.2024, 14, 983 5 of 12 models.Figure 3b shows that the results achieved good agreement.The model parameters used in the simulator are listed in Table 2.The simulation results are also verified by the commercial TCAD tool [35] for the detail distribution.The TCAD simulation utilized the same physical model as the simulator described above.For this simulation task, which involves nearly 6000 mesh points, our simulator, running on CentOS 7, requires roughly 200 MB of memory and completes the process in about three minutes on an Intel ® (Santa Clara, CA, USA) Xeon ® Gold 5218 Processor.
Figure 4 shows the simulation results of the 4H-SiC VDMOSFET using this simulator, where the potential distribution is at x = 20 µm in the vertical junction line between P Base The simulation results are also verified by the commercial TCAD tool [35] for the detail distribution.The TCAD simulation utilized the same physical model as the simulator described above.For this simulation task, which involves nearly 6000 mesh points, our simulator, running on CentOS 7, requires roughly 200 MB of memory and completes the process in about three minutes on an Intel ® (Santa Clara, CA, USA) Xeon ® Gold 5218 Processor.
Figure 4 shows the simulation results of the 4H-SiC VDMOSFET using this simulator, where the potential distribution is at x = 20 µm in the vertical junction line between P Base and N Drift, and the electron density distribution y = 12.515 µm in the horizontal midline of N Drift; the simulation using a commercial TCAD is also shown.
Appl.Sci.2024, 14, x FOR PEER REVIEW 6 of 13 and N Drift, and the electron density distribution y = 12.515 µm in the horizontal midline of N Drift; the simulation using a commercial TCAD is also shown.
From Figure 4, it can be seen that excellent agreements with potential distribution, electron density distribution, output characteristics and transfer characteristics are archived.

Analysis of the Lattice Temperature
To investigate the BV of the VDMOSFET, Ids-Vds is simulated with different lattice temperatures and various device structures.During simulation, the maximum Vds is 2500 V, Vgs is fixed to 10 V, and the break criteria of the current density is fixed to 1 × 10 −4 A/µm.Due to our reasonable choice of scale factor and calculation method for the scaled electron density, there are no convergence issues encountered caused by a low intrinsic carrier concentration.
Figure 5 shows the ON-state breakdown characteristic curves and the transfer characteristics curves.With the increasing in lattice temperature from 300 K to 573 K, the BV decreases from 2450 V to 2435 V, which indicates that the BV is relatively stable with the lattice temperature increases, as shown in Figure 5a.
When Vds is between 0 and 100 V, increasing Vgs forms a channel in the base region, causing electrons in the N+ region to flow to the N+ substrate.At Vgs = 10 V, with Vds ranging from 0 to 100 V, the on-resistance values at three temperatures are 1.313 mΩ•cm 2 , 1.421 mΩ•cm 2 and 1.475 mΩ•cm 2 , respectively.The formula of on-resistance is given as [5]  = (7) where εS is the semiconductor electrical permittivity, µn is the electron mobility and EC is the critical electric field.From Figure 4, it can be seen that excellent agreements with potential distribution, electron density distribution, output characteristics and transfer characteristics are archived.

Analysis of the Lattice Temperature
To investigate the BV of the VDMOSFET, Ids-Vds is simulated with different lattice temperatures and various device structures.During simulation, the maximum Vds is 2500 V, Vgs is fixed to 10 V, and the break criteria of the current density is fixed to 1 × 10 −4 A/µm.Due to our reasonable choice of scale factor and calculation method for the scaled electron density, there are no convergence issues encountered caused by a low intrinsic carrier concentration.
Figure 5 shows the ON-state breakdown characteristic curves and the transfer characteristics curves.With the increasing in lattice temperature from 300 K to 573 K, the BV decreases from 2450 V to 2435 V, which indicates that the BV is relatively stable with the lattice temperature increases, as shown in Figure 5a.
When Vds is between 0 and 100 V, increasing Vgs forms a channel in the base region, causing electrons in the N+ region to flow to the N+ substrate.At Vgs = 10 V, with Vds ranging from 0 to 100 V, the on-resistance values at three temperatures are 1.313 mΩ•cm 2 , 1.421 mΩ•cm 2 and 1.475 mΩ•cm 2 , respectively.The formula of on-resistance is given as [5] where ε S is the semiconductor electrical permittivity, µ n is the electron mobility and E C is the critical electric field.
The ionization coefficient, α, as in Equation ( 5), increases the electric field strength at a constant temperature.When Vds exceeds 2400 V, the device undergoes an avalanche breakdown due to the strong electric field strength.As the temperature rises, α exhibits a marginal decrease.Therefore, the generation rate in (4) undergoes little change, resulting in the breakdown voltage remaining relatively stable despite the increasing lattice temperatures within this temperature range.As depicted in Figure 5b, the transfer characteristics for the base region with a radius of 0 µm exhibit variations across different lattice temperatures, with the corresponding threshold voltage (Vth) and subthreshold swing (SS) values detailed in Table 3.The results in Table 3 show that with the increase in temperature, the change in Vth is less than 0.5 V, while the increase in SS mainly affects the leakage current, and they have little effect on the BV.The ON-state current density distributions of Vds = 2.4 kV and 2.5 kV and Vgs = 10 V when T = 300 K are plotted in Figure 6.In Figure 6a, it is evident that electrons flow from the channel into the JFET region, subsequently reaching the drain before breakdown.When the drain voltage is sufficiently high, i.e., a sufficiently high reverse voltage is applied at the collector junction in the parasitic NPN bipolar transistor at the bottom of the PBase region, as depicted in Figure 6b, a secondary breakdown occurs.This event is characterized by a rapid increase in the current density, resulting in the destruction of a MOSFET.In the range of 300 K to 573 K, the saturation current density decreases with increasing temperature as the device enters the saturation region.This reduction is attributed to enhanced lattice vibrations at elevated temperatures, which significantly diminish mobility and consequently lower the current density.
The ionization coefficient, α, as in Equation ( 5), increases the electric field strength at a constant temperature.When Vds exceeds 2400 V, the device undergoes an avalanche breakdown due to the strong electric field strength.As the temperature rises, α exhibits a marginal decrease.Therefore, the generation rate in (4) undergoes little change, resulting in the breakdown voltage remaining relatively stable despite the increasing lattice temperatures within this temperature range.
As depicted in Figure 5b, the transfer characteristics for the base region with a radius of 0 µm exhibit variations across different lattice temperatures, with the corresponding threshold voltage (Vth) and subthreshold swing (SS) values detailed in Table 3.The results in Table 3 show that with the increase in temperature, the change in Vth is less than 0.5 V, while the increase in SS mainly affects the leakage current, and they have little effect on the BV.The ON-state current density distributions of Vds = 2.4 kV and 2.5 kV and Vgs = 10 V when T = 300 K are plotted in Figure 6.In Figure 6a, it is evident that electrons flow from the channel into the JFET region, subsequently reaching the drain before breakdown.When the drain voltage is sufficiently high, i.e., a sufficiently high reverse voltage is applied at the collector junction in the parasitic NPN bipolar transistor at the bottom of the P Base region, as depicted in Figure 6b, a secondary breakdown occurs.This event is characterized by a rapid increase in the current density, resulting in the destruction of a MOSFET.

Evaluation of the Device Geometry
Since the lattice temperature has less impact on the BV, following simulation, we fixed the lattice temperature constant at 300 K. Figure 7 plots the BV curves for devices with different radii of the curvature at the base corner (R Base corner), as shown in Figure 2. The results show that a significant increase in the ON-state current occurring as the bottom of the P base region changes from a right angle to an arc.As the curvature of the arc increases, there is no significant change in the ON-state current.Alongside that, BV decreases significantly as the radius of curvature at the base corner increases when Vgs = 10 V and T = 300 K. Specifically, the BV with a radius of 5 µm undergoes a 21% decrease compared to that with a radius of 0 µm.The increase in curvature causes a higher current density, which increases the impact ionization rate, as shown in Equation ( 5), which in turn reduces the breakdown voltage.

Evaluation of the Device Geometry
Since the lattice temperature has less impact on the BV, following simulation, we fixed the lattice temperature constant at 300 K. Figure 7 plots the BV curves for devices with different radii of the curvature at the base corner (RBase corner), as shown in Figure 2. The results show that a significant increase in the ON-state current occurring as the bottom of the P base region changes from a right angle to an arc.As the curvature of the arc increases, there is no significant change in the ON-state current.Alongside that, BV decreases significantly as the radius of curvature at the base corner increases when Vgs = 10 V and T = 300 K. Specifically, the BV with a radius of 5 µm undergoes a 21% decrease compared to that with a radius of 0 µm.The increase in curvature causes a higher current density, which increases the impact ionization rate, as shown in Equation (5), which in turn reduces the breakdown voltage.Figure 8 illustrates the impact ionization rate distribution for different corner curvature radius at the base region.Figure 8a,b presents the impact ionization rate before and after breakdown for a radius of 0 µm, respectively.The peak impact ionization rate occurs near the top interface between the P base and N drift region, with peaks of 3.6 × 10 20 cm −3 /s and 7.7 × 10 23 cm −3 /s, respectively.The variation in the distribution of the impact ionization rate clearly elucidates the breakdown process as Vds increases from 2.4 kV to 2.5 kV at Vgs = 10 V.The contour where the impact ionization rate reaches 10 23 cm −3 /s, as demonstrated in Figure 8b, corresponds to the current flowing from the N+ region to the N+ substrate, as shown in Figure 6b. Figure 8c presents the impact ionization rate distribution

Evaluation of the Device Geometry
Since the lattice temperature has less impact on the BV, following simulation, we fixed the lattice temperature constant at 300 K. Figure 7 plots the BV curves for devices with different radii of the curvature at the base corner (RBase corner), as shown in Figure 2. The results show that a significant increase in the ON-state current occurring as the bottom of the P base region changes from a right angle to an arc.As the curvature of the arc increases, there is no significant change in the ON-state current.Alongside that, BV decreases significantly as the radius of curvature at the base corner increases when Vgs = 10 V and T = 300 K. Specifically, the BV with a radius of 5 µm undergoes a 21% decrease compared to that with a radius of 0 µm.The increase in curvature causes a higher current density, which increases the impact ionization rate, as shown in Equation ( 5), which in turn reduces the breakdown voltage.Figure 8 illustrates the impact ionization rate distribution for different corner curvature radius at the base region.Figure 8a,b presents the impact ionization rate before and after breakdown for a radius of 0 µm, respectively.The peak impact ionization rate occurs near the top interface between the P base and N drift region, with peaks of 3.6 × 10 20 cm −3 /s and 7.7 × 10 23 cm −3 /s, respectively.The variation in the distribution of the impact ionization rate clearly elucidates the breakdown process as Vds increases from 2.4 kV to 2.5 kV at Vgs = 10 V.The contour where the impact ionization rate reaches 10 23 cm −3 /s, as demonstrated in Figure 8b, corresponds to the current flowing from the N+ region to the N+ substrate, as shown in Figure 6b. Figure 8c presents the impact ionization rate distribution Figure 8 illustrates the impact ionization rate distribution for different corner curvature radius at the base region.Figure 8a,b presents the impact ionization rate before and after breakdown for a radius of 0 µm, respectively.The peak impact ionization rate occurs near the top interface between the P base and N drift region, with peaks of 3.6 × 10 20 cm −3 /s and 7.7 × 10 23 cm −3 /s, respectively.The variation in the distribution of the impact ionization rate clearly elucidates the breakdown process as Vds increases from 2.4 kV to 2.5 kV at Vgs = 10 V.The contour where the impact ionization rate reaches 10 23 cm −3 /s, as demonstrated in Figure 8b, corresponds to the current flowing from the N+ region to the N+ substrate, as shown in Figure 6b. Figure 8c presents the impact ionization rate distribution with a corner curvature of 5 µm.Notably, the contour range with a value of 10 23 cm −3 /s (dashed line in Figure 8b,c) is significantly broader than that in Figure 8b, indicating that increasing the curvature of the base region corner leads to the electrons being more concentrated on the lower side of the base region, thereby reducing the BV.
Figure 9 presents both the current density distribution associated with corner curvature radii of 5 µm and 1.2 µm at the base region, as well as the distribution of the electric field at y = 20 µm.An evident observation is the substantial widening of the current flow path in Figure 9a,b, in contrast to the case with a radius of 0 µm, as depicted in Figure 6b.Comparing the dashed lines in Figure 9a,b, a wider current density in the region of 15-17.5 µm can be seen when radius = 5 µm.This observation further explains the difference in BV when the radius exceeds 0 µm, as illustrated in Figure 7.The electric field distribution in Figure 9c indicates that the electric field strength does not decrease significantly as the radius increases.While there is a minor reduction in the peak electric field intensity at a radius of 5 µm, the overall electric field distribution demonstrates a leftward shift in the region of 10-25 µm.On the other hand, Figure 6b reveals that upon breakdown, the current density beneath the base region is predominantly distributed in the region of 12.5-17.5µm.Moreover, Figure 9c shows that within this region, the electric field intensity at a radius of 5 µm is higher than that at a radius of 1.2 µm.This accounts for the broadening of the current density beneath the Base Region when the radius is set to 5 µm.The alteration in radius modifies the electric field distribution, subsequently influencing the current density.with a corner curvature of 5 µm.Notably, the contour range with a value of 10 23 cm −3 /s (dashed line in Figure 8b,c) is significantly broader than that in Figure 8b, indicating that increasing the curvature of the base region corner leads to the electrons being more concentrated on the lower side of the base region, thereby reducing the BV. Figure 9 presents both the current density distribution associated with corner curvature radii of 5 µm and 1.2 µm at the base region, as well as the distribution of the electric field at y = 20 µm.An evident observation is the substantial widening of the current flow path in Figure 9a,b, in contrast to the case with a radius of 0 µm, as depicted in Figure 6b.Comparing the dashed lines in Figure 9a,b, a wider current density in the region of 15-17.5 µm can be seen when radius = 5 µm.This observation further explains the difference in BV when the radius exceeds 0 µm, as illustrated in Figure 7.The electric field distribution in Figure 9c indicates that the electric field strength does not decrease significantly as the radius increases.While there is a minor reduction in the peak electric field intensity at a radius of 5 µm, the overall electric field distribution demonstrates a leftward shift in the region of 10-25 µm.On the other hand, Figure 6b reveals that upon breakdown, the current density beneath the base region is predominantly distributed in the region of 12.5-17.5µm.Moreover, Figure 9c shows that within this region, the electric field intensity at a radius of 5 µm is higher than that at a radius of 1.2 µm.This accounts for the broadening of the current density beneath the Base Region when the radius is set to 5 µm.The alteration in radius modifies the electric field distribution, subsequently influencing the current density.
This finding signifies the regulatory influence of the device structure on the current density distribution, leading to a notable reduction in BV.The broader current flow path indicates a more intricate modulation of the carrier distribution, which contributes to the observed decrease in BV.These results provide valuable insights into the interplay between device geometry and BV, thereby facilitating the optimization and design of devices.

Evaluation of the Doping Profile
In addition to investigating the impact of the corner curvature structure, we also assess the influence of the doping profile on BV.We further evaluate the variation in Baliga's Figure-of-Merit (BFoM) [36] with the doping profile, which is a valuable factor in comparing the performance of power MOSFETs; it can be expressed as follows [36]: Figure 10 plots the BV and BFoM dependence on the doping concentration of the drift region (NDrift).There is a pronounced decrease in BV as the doping concentration of the drift region increases, as shown in Figure 10a.In Figure 10b, the BFoM initially de- This finding signifies the regulatory influence of the device structure on the current density distribution, leading to a notable reduction in BV.The broader current flow path indicates a more intricate modulation of the carrier distribution, which contributes to the observed decrease in BV.These results provide valuable insights into the interplay between device geometry and BV, thereby facilitating the optimization and design of devices.

Evaluation of the Doping Profile
In addition to investigating the impact of the corner curvature structure, we also assess the influence of the doping profile on BV.We further evaluate the variation in Baliga's Figure-of-Merit (BFoM) [36] with the doping profile, which is a valuable factor in comparing the performance of power MOSFETs; it can be expressed as follows [36]: Figure 10 plots the BV and BFoM dependence on the doping concentration of the drift region (N Drift ).There is a pronounced decrease in BV as the doping concentration of the drift region increases, as shown in Figure 10a.In Figure 10b, the BFoM initially decreases with the increasing dopant concentration, but changes very little as the dopant concentration exceeds 10 16 cm −3 .This is because the depletion region area reduces with the increasing doping concentration, which provides a wider current flow path.Therefore, the doping concentration in the N drift region should be as low as possible in order to achieve better device performance.

Evaluation of the Doping Profile
In addition to investigating the impact of the corner curvature structure, we also assess the influence of the doping profile on BV.We further evaluate the variation in Baliga's Figure-of-Merit (BFoM) [36] with the doping profile, which is a valuable factor in comparing the performance of power MOSFETs; it can be expressed as follows [36]: Figure 10 plots the BV and BFoM dependence on the doping concentration of the drift region (NDrift).There is a pronounced decrease in BV as the doping concentration of the drift region increases, as shown in Figure 10a.In Figure 10b, the BFoM initially decreases with the increasing dopant concentration, but changes very little as the dopant concentration exceeds 10 16 cm −3 .This is because the depletion region area reduces with the increasing doping concentration, which provides a wider current flow path.Therefore, the doping concentration in the N drift region should be as low as possible in order to achieve better device performance.Moreover, Figure 11 shows a significant reduction in BV when a Gaussian doping profile is incorporated between the drift region and the substrate.This observation highlights the substantial influence of increasing the dopant concentration in the N+ drift Moreover, Figure 11 shows a significant reduction in BV when a Gaussian doping profile is incorporated between the drift region and the substrate.This observation highlights the substantial influence of increasing the dopant concentration in the N+ drift region, which leads to a reduction in the depletion region and thereby affects the breakdown characteristics.

Conclusions
In order to improve the convergence in the numerical simulation of 4H-SiC devices, in this paper, we use a novel numerical processing method to develop a device simulator for 4H-SiC with high efficiency, high numerical accuracy and convergence.The simulator incorporates various physical models to investigate breakdown characteristics and device

Figure 1 .
Figure 1.Framework of the simulator.

Figure 1 .
Figure 1.Framework of the simulator.

Figure 2 .
Figure 2. Schematic cross-section of the 4H-SiC VDMOSFET used for the simulation.Green represents N-type doping, yellow represents P-type doping, and shades of color represent heavy and light doping.Brown for oxide and gray for contact.

Figure 2 .
Figure 2. Schematic cross-section of the 4H-SiC VDMOSFET used for the simulation.Green represents N-type doping, yellow represents P-type doping, and shades of color represent heavy and light doping.Brown for oxide and gray for contact.

Figure 5 .
Figure 5. (a) Breakdown characteristic curves (inset: output characteristics of Vds between 0 and 100 V at Vgs = 10 V) and (b) transfer characteristic curves with different lattice temperatures.

Figure 5 .
Figure 5. (a) Breakdown characteristic curves (inset: output characteristics of Vds between 0 and 100 V at Vgs = 10 V) and (b) transfer characteristic curves with different lattice temperatures.

Figure 7 .
Figure 7. Breakdown characteristic curves with different radii of the corner curvature at the base region at T = 300 K and Vgs = 10 V. (insets: schematic cross-sections with R of 0 µm and 5 µm).

Figure 7 .
Figure 7. Breakdown characteristic curves with different radii of the corner curvature at the base region at T = 300 K and Vgs = 10 V. (insets: schematic cross-sections with R of 0 µm and 5 µm).

Figure 7 .
Figure 7. Breakdown characteristic curves with different radii of the corner curvature at the base region at T = 300 K and Vgs = 10 V. (insets: schematic cross-sections with R of 0 µm and 5 µm).

Figure 8 .
Figure 8. Impact ionization rate distributions with radius of (a,b) 0 µm and (c) 5 µm of the corner curvature at the base region.The breakdown process is shown from (a,b).The dashed line in (b) and (c) represents the contour where the impact ionization rate equals to 1 × 10 23 cm −3 /s.

Figure 8 .
Figure 8. Impact ionization rate distributions with radius of (a,b) 0 µm and (c) 5 µm of the corner curvature at the base region.The breakdown process is shown from (a,b).The dashed line in (b) and (c) represents the contour where the impact ionization rate equals to 1 × 10 23 cm −3 /s.Appl.Sci.2024, 14, x FOR PEER REVIEW 10 of 13

Figure 9 .
Figure 9. Distribution of the current density with (a) radius of 5 µm and (b) radius of 1.2 µm of the corner curvature at the base region.(c) The simulated electric field at y = 20 µm (inset: schematic cross section showing y at 20 µm).The dashed lines in (a,b) represent the contour in the N drift region where the current density equals to 10 4 A/cm 2 .

Figure 9 .
Figure 9. Distribution of the current density with (a) radius of 5 µm and (b) radius of 1.2 µm of the corner curvature at the base region.(c) The simulated electric field at y = 20 µm (inset: schematic cross section showing y at 20 µm).The dashed lines in (a,b) represent the contour in the N drift region where the current density equals to 10 4 A/cm 2 .

Figure 9 .
Figure 9. Distribution of the current density with (a) radius of 5 µm and (b) radius of 1.2 µm of the corner curvature at the base region.(c) The simulated electric field at y = 20 µm (inset: schematic cross section showing y at 20 µm).The dashed lines in (a,b) represent the contour in the N drift region where the current density equals to 10 4 A/cm 2 .

Figure 10 .
Figure 10.The relationship of (a) BV and (b) BFoM with the doping profile of the drift region.

Figure 10 .
Figure 10.The relationship of (a) BV and (b) BFoM with the doping profile of the drift region.
Appl.Sci.2024, 14, x FOR PEER REVIEW 11 of 13 region, which leads to a reduction in the depletion region and thereby affects the breakdown characteristics.

Figure 11 .
Figure 11.BV versus the depth of the Gaussian doping profile.

Figure 11 .
Figure 11.BV versus the depth of the Gaussian doping profile.

Table 1 .
Device parameters for the simulations.: Parameters for validating the experiment; # : Parameters for other simulations.
*: Parameters for validating the experiment; # : Parameters for other simulations.*

Table 2 .
Model parameters for the simulations.

Table 3 .
Vth and SS values under different lattice temperatures.

Table 3 .
Vth and SS values under different lattice temperatures.