TCAD Device Modelling and Simulation of Wide Bandgap Power Semiconductors

Abstract

Technology computer-aided Design (TCAD) is essential for devices technology development, including wide bandgap power semiconductors. However, most TCAD tools were originally developed for silicon and their performance and accuracy for wide bandgap semiconductors is contentious. This chapter will deal with TCAD device modelling of wide bandgap power semiconductors. In particular, modelling and simulating 3C- and 4H-Silicon Carbide (SiC), Gallium Nitride (GaN) and Diamond devices are examined. The challenges associated with modelling the material and device physics are analyzed in detail. It also includes convergence issues and accuracy of predicted performance. Modelling and simulating defects, traps and the effect of these traps on the characteristics are also discussed.

1. Introduction

1.1. Materials and devices

Power devices made from wide bandgap materials have superior performance compared to those made from silicon. They can sustain higher voltages or endure smaller losses, they can operate at much higher junction temperatures and can switch faster. This is because, as seen in Figure 1, of their superior electrical properties. This in turn can revolutionize how applications work and possibly make new applications possible. The most mature materials in terms of process and device technology are the 4H-SiC and GaN-based devices. There is however, an increased interest and effort in advancing 3C-SiC and overcoming the technical challenges associated with the development of good-quality 3C-SiC wafers because of lower processing costs and the possibility for high channel mobility MOSFETs. Ultra-wide bandgap materials such as diamond are also being explored.

Figure 1.

Material properties of significant power semiconductors normalized against Si, data taken form [1, 2].

Silicon carbide is similar to silicon in terms of device design and optimization strategies. Most devices made in silicon can also be made in SiC. This includes Schottky Barrier Diodes, MOSFETs, IGBTs, and Thyristors. However, devices made in GaN are hetero devices and base their operation on the two-dimensional electron gas (2DEG) that is formed in the quantum well between the heterojunction interfaces. This quantum well provides electrons with a highly conductive channel, allowing high electron mobility. It is reported in literature that the electron mobility in GaN HEMT devices can be upwards of 1500 cm²/Vs [3].

1.2. Technology computer-aided design

TCAD is an engineering computer-aided tool which enables the physics-based modelling of semiconductor devices and their fabrication process. Due to the excellent predicting capability, semiconductor process and device engineers use TCAD for virtual prototyping and optimization of devices, to reduce the number of experimental cycles and, therefore, to reduce the production cost. TCAD can also be used to study the performance of devices when used in new applications or environments, to find performance limits and to analyze failures [4, 5]. Modern TCAD suites are compiled by several tools. A typical example of a TCAD suite is diagrammatically described in Figure 2 and it consists of the following elements:

• Device design tool

Figure 2.

Exemplary combination of tools making up a typical TCAD suite framework.

The device design tool allows the rapid creation of device structures via the use of scripting language or a graphical user interface without the need to know the process recipe. At this stage, the device geometry, the material and doping profile, and concentration of regions is defined. Commercial tools include Synopsys Sentaurus Structure Editor [6] and SilvacoDevEdit™ [7]. These tools allow device designers to parameterize device aspects and features to optimize their design or to assess the performance dependence on the parameters of question.

• Process simulation tool

The process simulation tools allow for the virtual fabrication of devices and the emulation of fabrication steps and conditions. They typically make use of a scripting language and require knowledge of a process recipe. These tools allow process engineers to fine tune their recipe and to analyze the effect of each process step and condition on the resulting device structure. Commercial tools include Synopsys Sentaurus Sentaurus Process [8] and Silvaco Athena [9].

• Device (and circuit) simulation tool

Device simulation tools have the capability to simulate the electrical, thermal, and optical properties and performance of devices. They can also account for the circuitry that surrounds a device when used in real applications. Therefore, they typically have SPICE capabilities too. They predict the device performance by executing finite element analysis and the solution of fundamental semiconductor physics equations. They make use of numerical devices created either through a device design tool or a process simulation tool. They take into consideration the materials incorporated in the device and they have a database with physics equations and the equivalent material parameters. Commercial tools include Synopsys Sentaurus Device [10] and Silvaco Atlas [11].

1.3. Summary

This chapter aims to highlight issues and present solutions and methods for achieving accurate models of WBG power devices. This includes proper modelling the material physics equations and their parameters, creating the finite elements and geometry and simulation of device characteristics including numerical issues and the effect of traps.

2. Material physics modelling

TCAD tools solve fundamental semiconductor physics equations to predict the performance of semiconductor devices. Accurately modelling the material properties through appropriate physics equations and parameters is essential for reliable simulations. These can depend on how advanced and mature the process, growth and technology is. For silicon, the material properties have been studied extensively and the technology is mature. Therefore, the material models and their parameters are typically considered accurate and able to predict the silicon-based devices performance with high level of accuracy. However, for WBG semiconductor-based devices, a lower confidence level exists. In the following subsections, we discuss the physics equations and the parameter sets required for accurate physical modelling of 3C-SiC, 4H-SiC, GaN and diamond power semiconductors. The equations used and parameters are compatible with most TCAD products.

2.1. Silicon carbide

There are hundreds of Silicon Carbide polytypes. From those, 4H-SiC is the most mature and studied polytype. The cubic polytype, 3C-SiC is also of particular interest because of its special properties, such as the ability to grow this compound semiconductor on large area silicon wafers, the lower temperatures required when processing the material and the ability to make MOSFETs with much higher channel mobility. The former two results in a much cheaper starting wafer whereas the latter one can enable the development of devices with higher efficiency. The corresponding values required in a physics parameter file for 4H-SiC were developed and improved alongside the improvement of the technology. Recently there have been efforts to compile an equivalent set of parameters for 3C and to validate them [12, 13]. This section presents the required models and the parameters of both materials. Wherever necessary, a range of values is used. Further, each physical mechanism presented is accompanied with the corresponding identified limitations. SiC polytypes can experience anisotropic properties, therefore when aiming multidimensional device simulation, these must be accounted for in the material parameter file and physics equations [14]. 4H-SiC experiences such behaviour, whilst 3C-SiC experiences isotropic behaviour. The parameters are then included within the ‘Device (and Circuit) Simulation’ tool of TCAD tools.

2.1.1. Band parameters

The desired properties of the WBG SiC originate from its bandgap characteristics. The value of $E_g$ shall not be considered constant, but variable with dependencies on both the concentration of impurities and temperature as shown in Figures 3 and 4 correspondingly. Increasing the temperature of the material essentially leads to lower gap values as described by Eq. (1). In addition, the phenomenon of bandgap narrowing causes band displacements in the energy scale directly related to doping according to Eq. (2), (3). These displacements represent potential barriers that may influence carrier transportation phenomena in the device. The bandgap dependence on doping in SiC can be described as in Eq. (4). The perception of the modeled effective bandgap value for the SiC can be expressed as $E_{g,\mathit{eff}}\left(T\right)=E_g\left(T\right)-E_{\mathit{bgn}}$ utilizing the parameters shown in Tables 1–3.

$$E_g\left(T\right)=E_g\left(0\right)-\alpha T^2/\left(T+\beta \right)$$E1

$$E_{\mathit{bgn}}^{\mathit{\text{cond}}}=\left\{\begin{array}{c}{A}_n\bullet N_{\mathit{tot}}^{1/3}+C_n\bullet N_{\mathit{tot}}^{1/2},N_{D,0}>N_{A,0}\\ B_p\bullet N_{\mathit{tot}}^{1/4}+D_p\bullet N_{\mathit{tot}}^{1/2},\mathit{\text{otherwise}}\end{array}\right.$$E2

$$E_{\mathit{bgn}}^{\mathit{val}}=\left\{\begin{array}{c}{B}_n\bullet N_{\mathit{tot}}^{1/4}+D_n\bullet N_{\mathit{tot}}^{1/2},N_{D,0}>N_{A,0}\\ A_p\bullet N_{\mathit{tot}}^{1/3}+C_p\bullet N_{\mathit{tot}}^{1/2},\mathit{\text{otherwise}}\end{array}\right.$$E3

$$∆E_g^0=-E_{\mathit{bgn}}^{\mathit{\text{cond}}}+E_{\mathit{bgn}}^{\mathit{val}}$$E4

Figure 3.

The 3C-SiC (top) and 4H-SiC (bottom) bandgap dependence on temperature as described with the models given in this chapter are in excellent agreement to literature data [15] for both β-SiC and α-SiC, respectively.

Figure 4.

The bandgap narrowing phenomenon as modeled for TCAD simulations assuming n-type SiC material. The measured energy displacements of the bands are sourced from [16].

Parameters description	3C-SiC	4H-SiC
Eg(0) [eV]	2.39	3.265
α [eV/K]	0.66 × 10^−3	3.3 × 10^−2
β [Κ]	1335	1.0 × 10^5

Table 1.

SiC parameter set related to bandgap.

Lindefelt model coefficients	3C-SiC semiconductor material
	n-type	p-type
$A_{n,p}$ [eV·cm]	−1.48 × 10^−8	1.3 × 10^−8
$B_{n,p}$ [eV·cm^3/4]	9.0 × 10^−7	−4.8 × 10^−7
$C_{n,p}$ [eV·cm^3/2]	−3.06 × 10^−12	1.43 × 10^−12
$D_{n,p}$ [eV·cm^3/2]	6.85 × 10^−12	−6.41 × 10^−13

Table 2.

3C-SiC band gap narrowing doping dependence.

Lindefelt model coefficients	4H-SiC semiconductor material
	n-type	p-type
$A_{n,p}$ [eV·cm]	−1.791 × 10^−8	3.507 × 10^−8
$B_{n,p}$ [eV·cm^3/4]	8.927 × 10^−7	−2.312 × 10^−6
$C_{n,p}$ [eV·cm^3/2]	−2.2 × 10^−12	6.74 × 10^−12
$D_{n,p}$ [eV·cm^3/2]	6.24 × 10^−12	−1.07 × 10^−12

Table 3.

4H-SiC bandgap narrowing doping dependence.

To model the intrinsic characteristics of the semiconductor, the temperature dependence of the density of states (DoS) for SiC in Eq. (5) is used. The parameters of Table 4 determine the intrinsic carrier concentration of the semiconductor as well as quantities like the effective carrier masses. In Figure 5, plotting the intrinsic carrier concentration against the temperature a discrepancy with measurements can be noticed for the case of 3C-SiC. This suggests that the TCAD simulations of 3C-SiC power devices utilizing this model should yield reasonably good results for limited temperature range of 200–500 K.

Parameters description	3C-SiC	4H-SiC
NC (300 K) [cm^−3]	1.5536 × 10^19	1.719 × 10^19
NV (300 K) [cm^−3]	1.1639 × 10^19	2.509 × 10^19

Table 4.

SiC temperature dependent density of states.

Figure 5.

The intrinsic carrier concentration as resulting from the model of DoS for both SiC cases in question. Comparison with literature data for 3C-SiC [18] and 4H-SiC [16] is performed. Assuming low doping levels (5 × 10¹⁵ cm⁻³) the bandgap narrowing is considered negligible.

Under low field conditions, the mobility of both types of carriers in SiC also depends on the doping concentration and on temperature. The first dependence is described by the Caughey-Thomas (C-T) model ${\mu }_0={\mu }_{\mathit{min}}+\left({\mu }_{\mathit{max}}-{\mu }_{\mathit{min}}\right)/\left(1+{\left(N/N_{\mathit{ref}}\right)}^{\alpha }\right)$ as illustrated in Figure 6. Each coefficient in the C-T equation changes with temperature as in Eq. (6) and the values in Table 6. Currently, there is an uncertainty for the holes’ mobility actual value in 3C-SiC. However, the values adopted in Table 5 are suggested from recent measurements [16, 17]. The mobility of carriers for the temperature range of 250–700 K was described in [12]. In high field conditions, with magnitude of values as shown in Figure 7, the mobility and velocity of carriers become inseparable directly affecting each other. The Canali model, Eq. (7) is utilized for these purposes and its parameters for SiC are presented in Table 7. In Eq. (8), the holes’ saturation velocity determines the range of electric field values at 200 kV/cm < E < 2000 kV/cm [12].

$$N_{C,V}\left(T\right)=N_{C,V}300\times {\left(T/300\right)}^{\raisebox{1ex}{$\text{3}$}\!\left/ \!\raisebox{-1ex}{$\text{2}$}\right.}$$E5

$$\mathit{Par}={\mathit{Par}}_0\times {\left(T/300K\right)}^{\gamma }$$E6

$${\mu }_0\left(E\right)=\frac{\left(\alpha +1\right){\mu }_{\mathit{low}}}{\alpha +{\left[1+{\left(\frac{\left(\alpha +1\right){\mu }_{\mathit{low}}E}{v_{\mathit{sat}}}\right)}^{\beta }\right]}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\beta $}\right.}}$$E7

$$v_{\mathit{sat}}=v_{\mathit{sat},0}{\left(\frac{300K}{T}\right)}^{v_{\mathit{sat},\mathit{exp}}}$$E8

Figure 6.

Based on experimental data [19], the SiC models utilized result in a very good accuracy. Increased the doping concentration, more scattering occurs and the mobility drops. The maximum carrier mobility values in SiC range and depend on crystal thickness [20] and impurities level.

Parameter	3C-SiC		4H-SiC [perpendicular to c-axis]		4H-SiC [parallel to c-axis]
	Electrons	Holes	Electrons	Holes	Electrons	Holes
μmax [cm^2/Vs]	650	70	910	114	1100	114
μmin [cm^2/Vs]	40	15	40	0	40	0
Nref [cm^−3]	1.5 × 10^17	5 × 10^19	2.0 × 10^17	2.4 × 10^18	2.0 × 10^17	2.4 × 10^18
α	0.8	0.3	0.76	0.69	0.76	0.69

Table 5.

SiC parameters for low field mobility and coefficients used to express doping dependence.

Parameter	3C-SiC		4H-SiC; [same for both directions]		Corresponding parameter from Table 5
	Electrons	Holes	Electrons	Holes
γmax	−0.3	−2.5	−2.4	−2.6	μmax [cm^2/Vs]
γmin	−0.5	−0.5	−1.536	−0.57	μmin [cm^2/Vs]
γNref	2	0	0.75	2.9	Nref [cm^−3]
γα	0	0	0.722	−0.2	α

Table 6.

SiC parameters for low field mobility and coefficients used to express temperature dependence.

Figure 7.

The models for SiC electrons’ velocity allow accurate TCAD simulations for temperatures up to 600 K. Comparison with literature data for 3C-SiC [21] and 4H-SiC [19] is performed. A doping level of 5 × 10¹⁵ cm⁻³ is assumed.

Parameter	3C-SiC		4H-SiC
	Electrons	Holes	Electrons [┴ to c-axis]	Electrons [// to c-axis]	Holes
beta (${\beta }_0$)	0.75	2.5	1.2	1.2	1.2
γbeta	−0.9	0	1.0	1.0	1.0
alpha ($\alpha$)	0	0	0	0	0
vsat,0 [cm/sec]	2.5 × 10^7	1.63 × 10^7	2.2 × 10^7	1.9 × 10^7	2.2 × 10^7
vsat,exp	−0.65	1.55	0.44	0.44	0.44

Table 7.

SiC parameters for high field mobility and saturation velocity along with coefficients used to express temperature dependence.

2.1.2. Doping and incomplete ionization

Compared to silicon, dopants in wide bandgap semiconductors have larger ionization energies, making activation of the doping species an issue. The dopant impurities are better modeled as traps to account for the phenomenon of incomplete ionization. As shown in Table 8, the formed energy levels depend on the polytype and the impurity. To model this behaviour, Eqs. (9) and (10) are utilized, where N_A,D is the doping concentration, N_A,D′ is the effective doping concentration, E_A,D and E_f are the activation energy, and the Fermi level, respectively, G_A is the degeneracy factor, and N_crit is the value that determines where the metallic-type conduction mechanism starts [22]. The degeneracy factor temperature dependence can be expressed by Eq. (11), whilst the typical values for g_A is 4 and for g_D is 2 for the impurity levels of acceptors and donors in SiC, respectively [23]:

$${\mathrm{N}}_{A,D^{\prime }}={\mathrm{N}}_{A,D}/\left(1+{\mathrm{G}}_A\ast \left({\mathrm{E}}_{A,D}-{\mathrm{E}}_f\right)/\mathrm{kT}\right)N_{A,D}<N_{\mathit{\text{crit}}}$$E9

$${\mathrm{N}}_{A,D^{\prime }}={\mathrm{N}}_{A,D}{N}_{A,D}\ge N_{\mathit{\text{crit}}}$$E10

$${\mathrm{G}}_{A,D}\left(T\right)={\mathrm{g}}_{A,D}·\mathit{exp}\left(\frac{\Delta E_{A,D}}{\mathit{kT}}\right)$$E11

Impurity	Species type	Energy levels [eV]
		3C-SiC	4H-SiC
Nitrogen (N)	Donora	0.057	0.071
Vanadium (V)	Donora	0.660	0.800
Aluminium (Al)	Acceptorb	0.260	0.265
Gallium (Ga)	Acceptorb	0.343	0.300
Boron (B)	Acceptorb	0.735	0.293

Table 8.

SiC impurities/shallow traps due to doping.

^a

The formed energy level is considered from the conduction band (E_C).

^b

The formed energy level is considered from the valence band (E_V).

2.1.3. Impact ionization

The impact ionization coefficients of electrons and holes need to be determined to specify the breakdown voltage [24]. For 3C-SiC, the Chynoweth law $\alpha \left(E_{\mathit{ava}}\right)=\gamma \bullet \alpha \bullet \mathit{exp}\left(-\mathit{\gamma b}/E_{\mathit{ava}}\right)$ can be adopted. The parameter values are determined in Table 9, following the work of van Overstraeten-de Man [25]. The temperature dependence of these parameters is expressed by determining the optical phonon energy as indicated by $\gamma =\text{tanh}\left(ħ{\omega }_{\mathit{op}}/2k{T}_0\right)/\text{tanh}\left(ħ{\omega }_{\mathit{op}}/2\mathit{kT}\right)$. Notably, it has been found in [26] that these values are relatively insensitive to temperature variations in the range of 300 K < T < 500 K. For 4H-SiC, a slightly different model is utilized after Hatakeyama’s work [27] to effectively describe the anisotropic behaviour of the avalanche coefficients. The avalanche force is considered to have two components to account for the anisotropic structure of 4H-SiC [28], satisfying $F^2=F_{//}^2+F_{┴}^2$. Utilizing the projections of these electric field components, the avalanche coefficients can be computed as indicated in Eqs. (12)–(15), with the default value of $\theta =1$:

$$a=a_{┴}^{{\left(\frac{B·F_{┴}}{b_{┴}·F}\right)}^2}·a_{//}^{{\left(\frac{B·F_{//}}{b_{//}·F}\right)}^2}$$E12

$$b=B\sqrt{1-\theta {{\rm A}}^2{\left(\frac{{\rm B}{F}_{┴}{F}_{//}}{F{b}_{┴}{b}_{//}}\right)}^2}$$E13

$$A=\mathit{log}\frac{a_{//}}{a_{┴}}$$E14

$$B=\frac{F}{\sqrt{{\left(\frac{F_{┴}}{b_{┴}}\right)}^2+{\left(\frac{F_{//}}{b_{//}}\right)}^2}}$$E15

Parameters description	Parameter name	3C-SiC		4H-SiC [┴ to c-axis]		4H-SiC [// to c-axis]
		Electrons	Holes	Electrons	Holes	Electrons	Holes
Ionization coefficients for electrons and holes	an,p [cm^−1]	1.28 × 10^6	1.07 × 10^7	1.76 × 10^8	3.41 × 10^8	2.1 × 10^7	2.96 × 10^7
	bn,p [V/cm]	5.54 × 10^6	1.12 × 10^7	3.3 × 10^7	2.5 × 10^7	1.7 × 10^7	1.6 × 10^7
Low field range up to this value	E0 [V/cm]	4.0 × 10^5	4.0 × 10^5	—	—	—	—
Optical phonon energy	$ħ{\omega }_{\mathit{op}}$ [eV]	0.120	0.120	0.190	0.190	0.190	0.190

Table 9.

SiC impact ionization coefficients for $T_0=300K$

2.1.4. Generation-recombination

The doping dependence of the SRH lifetime is modeled with Scharfetter in Eq. (16). The band-to-band non-radiative recombination is expressed with the Auger recombination rates, $R_{\mathit{net}}^A=\left(C_{n}n+C_{p}p\right)\left(\mathit{np}-n_{i,\mathit{eff}}^2\right)$. As shown in Eqs. (17) and (18), the effect of temperature and doping is accounted for. Typical values for the Scharfetter and Auger models are shown in Tables 10 and 11. However, since these are heavily process dependent, they need to be adjusted every time the process conditions change:

$${\tau }_{\mathit{dop}}={\tau }_{\mathit{min}}+\frac{{\tau }_{\mathit{max}}-{\tau }_{\mathit{min}}}{1+{\left(\frac{N_{A,0}+N_{D,0}}{N_{\mathit{ref}}}\right)}^{\gamma }}$$E16

$$C_n\left(T,n\right)=\left(A_n+B_n\left(\frac{T}{T_0}\right)+C_n{\left(\frac{T}{T_0}\right)}^2\right)\left[1+H_n\mathit{exp}\left(-\frac{n}{N0,n}\right)\right]$$E17

$$C_p\left(T,p\right)=\left(A_p+B_p\left(\frac{T}{T_0}\right)+C_p{\left(\frac{T}{T_0}\right)}^2\right)\left[1+H_p\mathit{exp}\left(-\frac{p}{N0,p}\right)\right]$$E18

Parameter	SiC
	Electrons	Holes
τmin [sec]	0	0
τmax [sec]	2.5 × 10^−6	0.5 × 10^−6
Nref [cm^−3]	1 × 10^17	1 × 10^17
γ	0.3	0.3
Tα	1.72	1.72
Τcoef	2.55	2.55
Etrap [eV]	0	0

Table 10.

SiC SRH lifetime parameter set.

Parameter	3C-SiC		4H-SiC
	Electrons	Holes	Electrons	Holes
A [cm^6/sec]	0.3 × 10^−31	0.2 × 10^−31	6.7 × 10^−32	7.2 × 10^−32
B [cm^6/sec]	0	0	2.45 × 10^−31	4.5 × 10^−33
C [cm^6/sec]	0	0	−2.2 × 10^−32	2.63 × 10^−32
H	1.9	1.9	3.47	8.26
N0 [cm^−3]	1.0 × 10^19		1.0 × 10^18
T0 [K]	300	300	300	300

Table 11.

SiC Auger recombination rates.

2.2. Gallium nitride

GaN-based devices utilize GaN, AlGaN and AlN materials. AlGaN is a molar fraction of AlN and GaN. GaN and AlGaN naturally form Wurtzite crystal structures with the ability of forming different Ga and N faces. For this chapter, only the Ga-face is considered. In TCAD, it is necessary to define the material properties of AlN and GaN separately. AlGaN material properties are thereafter approximated through a linear interpolation, depending on the molar fraction of AlN and GaN. A more accurate result can be yielded wherever the molar compounds are known and experimental evaluation of their properties has been performed. In those cases, new material parameters can be made which would reflect the exact molar compound of interest, which in turn would give more accurate simulations.

Similar to modelling SiC, the important parameters include modelling the bandgap, doping and incomplete ionization, impact ionization, mobility and generation-recombination.

2.2.1. Band parameters

Like all III-Nitride semiconductors, GaN and AlN are direct bandgap semiconductors, that is, the maximum valley in their valance band is directly below the minimum conduction valley. Similar to SiC examined earlier, E_g cannot be assumed constant in value and will vary heavily with temperature and doping concentration. Again, like SiC, thermally exciting GaN leads to lower band gap energies. According to Eq. (1) the variation of bandgap value can be described utilizing the parameters in Table 12 [29].

Parameters description	GaN	AlN
Eg(0) [eV]	3.507	6.23
α [eV/K]	9.1×10^−4	1.79×10^−3
β [Κ]	836	1.462×10^3

Table 12.

GaN/AlN parameter set related to bandgap.

Further, utilizing Eq. (19), the intrinsic carrier concentration of the GaN material can be calculated. This premises that the DoS and/or the effective carrier masses in this WBG material are known. Equations (20) and (21) adequately describe this temperature dependent procedure. The intrinsic characteristics of GaN can then be given by [30]:

$${\mathrm{n}}_{\mathrm{i}}={\left({\mathrm{N}}_{\mathrm{c}}·{\mathrm{N}}_{\mathrm{v}}\right)}^{1/2}\text{exp}\left(-{\mathrm{E}}_{\mathrm{g}}/\left({\text{2k}}_{\mathrm{B}}\mathrm{T}\right)\right)$$E19

$${\mathrm{N}}_{\mathrm{c}}˜=4.82{10}^{15}·{\left({\mathrm{m}}_{\mathrm{\Gamma }}/{\mathrm{m}}_0\right)}^{3/2}{\mathrm{T}}^{3/2}\left({\mathrm{cm}}^{-3}\right)˜=4.3\times {10}^{14}\times {\mathrm{T}}^{3/2}\left({\mathrm{cm}}^{-3}\right)$$E20

$${\mathrm{N}}_{\mathrm{v}}=8.9\times {10}^{15}\times {\mathrm{T}}^{3/2}\left({\mathrm{cm}}^{-3}\right)$$E21

2.2.2. Mobility

The low field carriers’ mobility for the bulk WBG compound depends on the carriers’ density following the C-T formula, as illustrated in SiC physical model. However, variations of this model exist to fit better some compound semiconductors behaviour. In SiC above, the so-called Arora model is utilized originating from the C-T model. It provides additional control on the temperature dependence of SiC low field mobility parameters. For GaN and AlN, the Masetti model, in Eq. (22), is preferred to describe the doping dependence, with the fitting parameters shown in Table 13. In this table, ${\gamma }_{{\mu }_{\mathit{max}}}$ indicates the coefficient for the temperature dependency of the maximum mobility parameter, as described by Eq. (6). This is the only temperature dependent parameter that the Masetti model accounts for:

$${\mu }_{\mathit{dop}}={\mu }_{\mathit{min}1}\text{exp}\left(\frac{P_c}{N_{A,0}+N_{D,0}}\right)+\frac{{\mu }_{\mathit{\text{const}}}-{\mu }_{\mathit{min}2}}{1+{\left(\frac{N_{A,0}+N_{D,0}}{C_r}\right)}^{\alpha }}-\frac{{\mu }_1}{1+{\left(\frac{C_s}{N_{A,0}+N_{D,0}}\right)}^{\beta }}$$E22

Parameter	GaN		AlN
	Electrons	Holes	Electrons	Holes
μconst [cm^2/Vs]	1500 [31] – 1800 [29]	20	300	14
${\gamma }_{{\mu }_{\mathit{max}}}$	1	2.1	1	2.1
${\mu }_{\mathit{min}1}$ [cm^2/Vs]	85	33	20	11
${\mu }_{\mathit{min}2}$ [cm^2/Vs]	75	0	65	0
${\mu }_1$ [cm^2/Vs]	50	20	20	10
$P_c$ [cm^−3]	6.5×10^15	5×10^15	8×10^17	5×10^18
$C_r$ [cm^−3]	9.5×10^16	8×10^16	7×10^16	8×10^17
$C_s$ [cm^−3]	7.2×10^19	8×10^20	5.2×10^17	8×10^18
α	0.55	0.55	0.88	1.05
β	0.75	0.7	0.75	0.75

Table 13.

Coefficients used for temperature and doping (Masetti model) dependency of low field mobility in GaN and AlN.

The calculated mobility in low field conditions is utilized in the Canali model, as discussed in SiC section with Eq. (7). The parameter values presented in Table 14 enable modelling the mobility in high field conditions while taking into consideration the temperature dependence with Eq. (6).

Parameter	GaN		AlN
	Electrons	Holes	Electrons	Holes
beta (${\beta }_0$)	7.2044	4	17.3681	4
γbeta	6.1973	0	8.7253	0
alpha ($\alpha$)	0.7857	0	0.8554	0
vsat,0 [cm/sec]	1.3×10^7	1.7×10^7	1.5×10^7	1.25×10^7
vsat,exp	0.7	0.725	2	2

Table 14.

GaN/AlN parameters for high field mobility and saturation velocity along with coefficients used to express temperature dependence.

2.2.3. Doping and incomplete ionization

Dopants in this WBG semiconductor would not fully ionize even at high temperatures [32] because the impurities form deep levels, as further observed in Table 15 [33, 34, 35, 36, 37, 38]. The models as presented in Eqs. (9)–(11) can be utilized here following the same degeneracy factor values for the conduction and valence bands [39]. That is, g_A equals 4 for shallow acceptors and g_D equals 2 for shallow donors.

Impurity	Species type	Energy levels [eV]
		GaN		AlN
		Ga	N	Al	N
Silicon (Si)	Donora	0.017	—	—	—
Nitrogen (N)	Donor^a	—	—	1.4–1.85	—
Vacancy (VN)	Donor^a	—	0.03–0.1	—	0.17
Carbon (C)	Donor^a	0.11–0.14	—	0.2	—
Magnesium (Mg)	Donor^a	—	0.26, 0.6	—	—
Vacancy (VGA)	Acceptorb	0.14	—	—	—
Silicon (Si)	Acceptor^b	0.19	—	—	—
Magnesium (Mg)	Acceptor^b	—	0.14–0.21	0.1	—
Zinc (Zn)	Acceptor^b	0.21–0.34	—	0.2	—
Mercury (Hg)	Acceptor^b	0.41	—	—	—
Cadmium (Cd)	Acceptor^b	0.55	—	—	—
Beryllium (Be)	Acceptor^b	0.7	—	—	—
Lithium (Li)	Acceptor^b	0.75	—	—	—
Carbon (C)	Acceptor^b	—	0.89	—	0.4
Gallium (Ga)	Acceptor^b	—	0.59–1.09	—	—
Aluminium (Al)	Donor^a	—	—	—	3.4–4.5
Vacancy (VAl)	Acceptor^b	—	—	0.5	—

Table 15.

Impurities and shallow traps due to doping in GaN/AlN

^a

The formed energy level is considered from the Conduction band (E_C)

^b

The formed energy level is considered from the Valence band (E_V)

2.2.4. Impact ionization

The Van Overstraetan-de Man expression, described in SiC section, with the parameters in Table 16 can be used to model the impact ionization phenomenon [29]. It is worth noticing that for the case of high electric fields (i.e. larger than the E₀ value) the impact ionization coefficients are predicted to be similar in GaN [40].

Parameters description	Parameter name	GaN		AlN
		Electrons	Holes	Electrons	Holes
Ionization coefficients for electrons and holes	an,p [cm^−1]	1.5×10^5	6.4×10^5	2.9×10^8	1.34×10^7
	bn,p [V/cm]	1.41×10^7	1.46×10^7	3.4×10^8	2.03×10^8
Low field range up to this value	E0 [V/cm]	4.0×10^5	4.0×10^5	4.0×10^5	4.0×10^5
Optical phonon energy	$ħ{\omega }_{\mathit{op}}$ [eV]	0.035	0.035	0.035	0.035

Table 16.

GaN/AlN impact ionization coefficients for $T_0=300K$

2.2.5. Generation-recombination

The non-radiative recombination [41], SRH, is described by the Scharfetter model in Eq. (15) is considered a dominant process in the bulk. Table 17 includes the relevant parameter values. Auger recombination process is realized utilizing Eqs. (17)–(18) along with the coefficient values shown in Table 18 [42]. Auger recombination in nitrides is responsible for the loss of quantum efficiency in InGaN-based light emitters [43]. This non-radiative loss mechanism is due to the large values of the Auger coefficients (2 × 10⁻³⁰ cm⁶/s) for specific parts of the emission spectrum, like blue to green region.

Parameter	GaN		AlN
	Electrons	Holes	Electrons	Holes
τmin [sec]	0	0	0	0
τmax [sec]	0.7×10^−11	2.0×10^−11	1.0×10^−9
Nref [cm^−3]	1×10^16
γ	1	1	1	1
Tα	−1.5	−1.5	−1.5	−1.5
Τcoef	2.5	2.5	2.5	2.5

Table 17.

GaN/AlN SRH lifetime parameter set.

Parameter	GaN
	Electrons	Holes
A [cm^6/sec]	3.0×10^−31	3.0×10^−31
B [cm^6/sec]	0	0
C [cm^6/sec]	0	0
H	0	0
N0 [cm^−3]	1.0×10^18	1.0×10^18
T0 [K]	300	300

Table 18.

GaN Auger recombination rates.

2.3. Diamond

Diamond (C) is considered to be a strong contester for high power due to its outstanding electrical and thermal material properties [44]. However, the realization of power semiconductor devices based on this material is extremely difficult and such devices are currently at the research level. The main reason is that there are not enough activated carriers at room temperature, leading to poor device performance. Furthermore, diamond is also extremely expensive to fabricate due to large scale material growth constrains.

It is worth mentioning that since diamond is at very early stages of development, it is also at very early stage at the material characterization and modelling. Currently, diamond does not exist in the material libraries of the commercial TCAD simulation packages. As a result, all material properties and fitting parameters of various materials interface with diamond need to be added manually.

2.3.1. Band parameters

The experimentally extracted value of diamond bandgap (Eg) is about 5.47 eV at room temperature [45]. This value directly translated into high material critical electric strength. High breakdown field strength means that the material can withstand high potential drops across very thin layers thus minimizing the on-state resistance of the device allowing for the fabrication of highly energy efficient high voltage high current devices.

The values for the density of states for the valence and the conduction band (Nv,Nc) are given by expressions Eq. (5) where in this case N_v300 = 1.8e19 and N_c300 = 5e18 [cm⁻³], as reported in [22, 44, 45, 46, 47]. Therefore, the appropriate levels of density of state could be calculated using Eq. (19). For CVD diamond this value is around 1.2 × 10⁻²⁷ cm⁻³ [22] at 300 K, Eg is the bandgap, T is the absolute temperature in Kelvin and k the Boltzmann constant (1.38 × 10⁻²³ J/K). This value is extremely small for any meaningful numerical analysis using TCAD simulations. Therefore, possible strategies to facilitate and match experimental results include adding a fitting coefficient/value for intrinsic concentration, activating the constant carrier generation models or including trap levels and recombination centres in the materials bandwidth [22].

2.3.2. Doping and incomplete ionization

Diamond doping either of n-type or p-type to a less extent is challenging; in order to develop accurate and reliable simulation models one has to include the incomplete ionization models with the appropriate coefficients. It is therefore necessary to use the ‘incomplete ionization model’ [48] (Figure 8). The equations implemented follow the model described in Eqs. (9) and (10). Note that the degeneracy factor in this particular case is a function of temperature, given by the model in Eq. (23):

$$g_{A,D}=4+2\mathit{exp}\left(-\Delta g_{A,D}/\left(\mathit{kT}\right)\right)$$E23

Figure 8.

(a) Incomplete ionization model implemented for boron in single crystal diamond and (b) variation of activation energy with doping concentration of boron [22].

2.3.3. Mobility

For an intrinsic diamond, the hole mobility is up to 3800 cm²/(Vs) whereas for boron-doped diamond the mobility is dominated by the concentration dependent scattering. At elevated temperatures, above 325 K it has been shown that the mobility is decreased due to acoustic phonon scattering [48, 49]. The combination of p+ doped diamond leads to lattice scattering. It therefore is advisable that these complex relationships between mobility, doping and temperature are taken into consideration when diamond simulations are implemented. A fairly accurate model, the unified mobility model (UniBo) [50] accounts for the dependence of the mobility on the doping and on the temperature using with a single model.

2.3.4. Impact ionization

Impact ionization integral coefficients for electron and holes are very important in the design of diamond devices with the appropriate values been included in the impact ionization models parameters. Under reverse bias conditions, the increase of the ionization integral towards unity indicates the increase in the generation of avalanche carriers due to impact ionization, leading eventually to the device breakdown. The avalanche coefficients are given by Eq. (24) where in which a_n,p, b_n,p and c_n,p are fitting parameters and E is the lateral electric field through the semiconductor layer. The numerical values for both diamond and silicon are given in Table 19:

$$\mathrm{{\rm A}}={\mathrm{a}}_{\mathrm{n},\mathrm{p}}\text{exp}\left(-{\mathrm{b}}_{\mathrm{n},\mathrm{p}}/\mathrm{E}\right){\mathrm{c}}_{\mathrm{n},\mathrm{p}}$$E24

Ionization coefficients for electrons and holes	an (cm^−1)	ap (cm^−1)	bn (Vcm^−1)	bp (Vcm^−1)
Diamond [22]	1.89×10^5	5.48×10^6	1.7×10^7	1.42×10^7

Table 19.

Avalanche coefficients for diamond.

3. Device modelling and simulation

Due to the very low intrinsic carrier concentration of WBG semiconductors, the expected leakage current is very low ($\sim {10}^{-20}A{\mathit{cm}}^{-2}$), which causes converge issues. To alleviate from this, the numerical precision should be significantly high. This is achieved with the inclusion of certain keywords in the ‘Device (and Circuit) Simulation’ tool command file (e.g. ‘ExtendedPrecision’ for SDevice tool of the Synopsys Sentaurus TCAD). Simulations that use extended arithmetic precision are computational more intensive, therefore, the arithmetic precision should be increased in a trade-off manner up to a level that is able to provide a solution.

A further method utilized to improve convergence issues, especially when simulating the blocking characteristics, is to add extra carriers’ generation through an equivalent keyword in the command file. The latter can increase the leakage current to levels that are high enough (e.g. $>{10}^{-10}A{\mathit{cm}}^{-2}$) for the simulations to converge at lower precision level. This method not only helps the simulations convergence but it also helps with accounting the realistic leakage currents in WBG devices, which are orders of magnitudes higher than the predicted ideal ones. The reason for the much higher real leakage currents is attributed to the immaturity of the technology, the high defects and traps density and due to background irradiation. Their effect on the leakage current is therefore, considered with this constant carrier generation statement.

The choice of solvers is also important when simulating WBG devices. Some linear solvers will be more suited for small to medium 2D simulations and others for medium to large 3D simulations due to their superior parallel performance and significantly smaller memory usage. In addition, relaxing the numerical setting for the linear solver constitutes another trade-off which may improve the operation of the solver, however, convergence complications may come as a cost.

3.1. Silicon carbide

For accurate simulation results, modelling the defects and traps is imperative. They directly influence the performance and strongly affect its reliability [51]. To highlight the effect of traps on the device performance and electrical characteristics, the P-i-N rectifier structure of Figure 9 was prepared for modelling and simulation. A Gaussian energetic distribution of deep levels is considered distributed spatially uniformly.

Figure 9.

Simplified device schematic used for simulations of a SiC P-i-N power rectifier following the fabricated diode in [52]. The Gaussian energetic distribution of the deep levels considered is also illustrated.

For most applications, the linear region of the forward I-V characteristics is important [14], the sub-threshold characteristics, however, are indicators of the material quality [53]. Both regions of the device characteristics are affected by the presence of traps. Figure 10 depicts the effect of the traps capture cross section, whereas Figure 11 depicts the effect of the concentration of deep level traps on the sub-threshold current-voltage (IV) forward characteristics. Because generation-recombination is the main carrier transport mechanism in the sub-threshold region, an increased concentration of deep levels or a larger capture-cross-section results in a higher magnitude for leakage current [14].

Figure 10.

Sub-threshold region of a SiC P-i-N structure with 10 um drift region at 1.4e16 cm⁻³. The defects concentration is fixed with a variable cross section value. The simulated 4H-SiC material contains donor traps (left), whereas the 3C-SiC contains acceptor traps (right).

Figure 11.

Sub-threshold region of a SiC P-i-N structure with 10 um drift region at 1.4e16 cm⁻³. The defects’ cross section is fixed with a variable concentration value. The simulated 4H-SiC material contains donor traps (left), whereas the 3C-SiC contains acceptor traps (right).

The forward the linear region is governed by the recombination-generation and the drift-diffusion. In this case, the defects have an opposite effect on the IV characteristics. An increased concentration of the traps intensifies carriers scattering, thus effectively reducing the mobility of carriers, which in turn leads to a decreased conductivity. This behaviour is less significant at elevated temperatures as the carriers gain enough kinetic energy to remain unaffected from the presence of nearby defects in the bulk. Consequently, the on-resistance of the material can decrease, as illustrated in Figure 12.

Figure 12.

The impact of carriers’ scattering lowers with increased temperature. This behaviour is independent of the type and the concentration of the deep level as can be seen for the 4H-SiC (left) and the 3C-SiC (right) case.

Traps also need to be included when modelling the interfaces between SiC and other materials. The traps’ energetic profile at or near the interface in those cases need to be identified and modeled appropriately. For SiC Schottky interfaces, the combined effect of tunnelling and traps is modeled with the quantum tunnelling and trap-assisted tunnelling models [54]. The effect of traps also needs to be modeled at the SiC/SiO2 interface, in particular when modelling SiC MOSFETs [55].

3.2. Gallium nitride

The GaN High Electron Mobility Transistor (HEMT) is regarded as the most successful attempt at harnessing the superior material properties of Gallium Nitride. The AlGaN HEMT is a heterostructure formed through the union of Al_XGa_1-XN and GaN. The inherent spontaneous and mechanically induced piezoelectric polarization charges a dictate the formation of a 2D Electron Gas across the heterojunction interface [56]. It is worth noting that the 2DEG channel is inherently present across the device and therefore, means that the device is naturally on. This has caused engineers to develop normally off GaN HEMT one being the p-type Gate GaN HEMT device. This device contains a p-typed GaN region beneath the gate which depletes the 2DEG of carriers and therefore, effectively stops the channel underneath the gate. In this section, device modelling will focus upon this device. The complexities associated with modelling GaN HEMT devices in TCAD are summarized below [57]:

• Reduced crystal symmetry compared to silicon due to the Wurtzite crystal structure.

• The polarization charges and subsequent effects on the device performance.

• The lower intrinsic carrier concentration associated with wide bandgaps semiconductors.

• The exchange of the inter-valley at high field that modulates the current through the Gunn effect.

• The abrupt nature of the heterojunction between the semiconductors and partially floating regions.

• The significant and extensive quantity of traps and their subsequent characteristics.

The forward operation of the GaN HEMT device is dependent on the characteristics of the 2DEG channel. There are multiple methods that can be employed to model the 2DEG channel. The first is the placement of an interface charge across the AlGaN/GaN interface. The second is a simplified strain model in which the TCAD calculates the polarization charges based on the material which calculates the spontaneous and piezoelectric charges through the molar fraction of the AlGaN layer and the corresponding strain based on the GaN layer. Another technique is applying a full strain model which calculates the polarization charges and are expressed through the local strain tensor. Finally, a stress model can be applied which calculates the polarization charges and are expressed in the local stress tensor.

TCAD simulations with the simplified strain model utilized, describe how the threshold voltage changes with the thickness of AlGaN. This behaviour is in accordance to Eq. (25) [58] where; V_th is threshold voltage, ${\Phi }_B$ is the height of the Schottky barrier,$E_C$is the conduction band offset, d is the thickness of the AlGaN barrier, N_D is the 2DEG concentration and ε is the relative dielectric constant of AlGaN. To increase the threshold voltage of the GaN HEMT, three parameters can be changed, the work function of the Schottky contact, the aluminium mole concentration in the AlGaN and the thickness of the AlGaN region. In this model, we have chosen to decrease the thickness of the AlGaN barrier whilst maintaining the same parameters for the rest of the device. Figure 13 depicts the simplified equivalent schematic representation of the simulated device and the transfer characteristics for three different AlGaN thicknesses of 0.20, 0.23 and 25 μm. As shown, the threshold of the device varies even with a very small change in the AlGaN thickness. This also demonstrates how sensitive the 2DEG is to the process. For that reason, a fixed interface charge across the AlGaN/GaN interface, is many times the preferred modelling approach, the concentration of which can be used as a fitting parameter:

$$V_{\mathit{th}}={\Phi }_B-\frac{∆E_C}{q}-{\iint }_d\frac{q{N}_D\left(x\right)}{\epsilon }{\mathit{dx}}^2$$E25

Figure 13.

Simplified schematic representation of the simulated HEMT device (left) and its simulated transfer characteristics for 0.20, 0.23 and 0.25 μm AlGaN barrier layers (right).

4. Conclusions

Adequate modelling and simulation of WBG power devices and their performance with TCAD presents challenges and complexities. It includes modelling the material physical properties, improving the numerical accuracy of simulations, taking special care for the device structure design and incorporating the effect of defects, often in the form of traps. It also includes understanding the complexities and trade-offs between convergence and simulation speed, and how these are affected by the choice of solvers and numerical accuracy. This chapter gave an overview of those for 3C-SiC, 4H-SiC, GaN and diamond-based devices.

Acknowledgments

The authors would like to thank the Research Institute for Future Cities and Transport, Coventry University and the Royal Society for funding this research.