A New Anisotropic Driving Force Model for SiC Device Simulations

We present a new anisotropic driving force model to compute the anisotropy of the mobility and impact ionization rate in TCAD simulations of silicon carbide (SiC) devices. We derive the model from the relationship between mobility, driving force, and carrier heating. The obtained driving force along the <inline-formula> <tex-math notation="LaTeX">${c}$ </tex-math></inline-formula>-axis and that perpendicular to the <inline-formula> <tex-math notation="LaTeX">${c}$ </tex-math></inline-formula>-axis depend on the gradient of the quasi-Fermi potential and the mobility anisotropy factor. Since mobility also depends on the driving force, we solve for the mobility anisotropy factor to determine the anisotropic mobility and driving force together. We also propose a new interpolation scheme to consider the anisotropy of the impact ionization rate using the derived anisotropic driving force. We perform the simulation of a 4H-SiC vertical power MOSFET device and compare the present model with the existing models in the commercial simulator.

Hexagonal SiC polytypes like 4H-SiC and 6H-SiC have anisotropic material properties [13].In particular, the anisotropy of the carrier mobility [4], [14], [15] and the impact ionization rate [9], [16] significantly impact the terminal current and the breakdown characteristics.Therefore, it is essential to accurately model the anisotropic material properties in the TCAD device simulations.
In the literature, the anisotropy of the carrier mobility is typically modeled by a constant factor [4], [15].However, this approach cannot be applied when the mobility anisotropy factor (the mobility ratio along the perpendicular direction to that along the c-axis) depends on the bias conditions.For example, the low-field mobility of 4H-SiC along the c-axis is about 1.2 times larger than that perpendicular to the c-axis [8], [12].While the experimental value is not available, Monte Carlo simulations and anisotropy of the impact ionization rate suggest that the saturation velocity along the c-axis is about 0.6-0.8times smaller than that perpendicular to the c-axis [8], [12], [17].As a result, the mobility anisotropy factor would be significantly affected by the magnitude of applied electric field.
In the homogenous bulk condition, we can express the carrier mobility as a function of the electric field [18].In the drift-diffusion (DD) model, a physical quantity called driving force is introduced as a generalization of the applied electric field to compute the carrier mobility in nonhomogeneous conditions.For isotropic materials, the gradient of the quasi-Fermi potential is commonly used as a driving force [19].For anisotropic materials, we need to include the dependence of the mobility anisotropy factor on the driving force as the mobility, the mobility anisotropy factor, and the driving force are interrelated.As a result, we cannot directly use the gradient of the quasi-Fermi potential as a driving force for the mobility model.In the commercial device simulator [20], anisotropic mobility models are available to account for the dependence of the anisotropy factor on the driving force.However, the anisotropy of the driving force itself has never been studied in the literature.
This article presents a new anisotropic driving force model for the mobility and impact ionization rate calculations.Section II describes the anisotropic driving force model derived from the relationship between mobility, driving force, and carrier heating.We also propose a new interpolation scheme to consider the anisotropy of the impact ionization rate using the derived anisotropic driving force.We compare our model with the existing models in the commercial simulator in detail.We simulate a 4H-SiC vertical power MOSFET as an application example in Section III.

A. Derivation of Anisotropic Driving Forces
Our derivation of the anisotropic driving force model is inspired by [9].For simplicity, we consider a 2-D domain, where the c-axis is along the x-direction, but the model applies to 3-D simulations as well.In the derivation, we use the subscripts "A" and "I" to represent the "anisotropic" (along the c-axis) and "isotropic" (perpendicular to the c-axis) directions.We consider a homogenous bulk system in the derivation.Assume that a uniform driving field is applied along the x-direction, F = (F A , 0), such that the electron velocity is given by v = (−µ A (F A )F A , 0) where µ A is the electron mobility along the c-axis.Then, the relation between the driving field and the carrier temperature is given by [21], [22] where k B is the Boltzmann's constant, T n the electron temperature, T L the lattice temperature, τ ω the energy relaxation time, and q the elementary charge.Now, consider that a uniform driving field is applied to the y-direction, F = (0, F I ), where the magnitude of F I is adjusted to give the same electron temperature to (1).Then, the following relation holds: where µ I is the electron mobility perpendicular to the c-axis.
Finally, when we apply a field vector in arbitrary direction F = F x , F y such that the carrier temperature is the same as ( 1) and ( 2), the electron velocity becomes where the mobility along the x-(y-)direction is equal to that in (1) [( 2)] as the carrier temperature is unchanged.Therefore, we obtain [9] |F Hatakeyama et al. [9] further relate (4) to the anisotropic impact ionization rate to derive the anisotropic impact ionization model, assuming that the impact ionization rate should be independent of the field direction when the carrier temperature is constant.
Instead, here we rearrange (4) to obtain the desired anisotropic driving force expressions for the given field vector F = F x , F y as follows: where r n is the mobility anisotropy factor defined as follows: Therefore, for the given driving field F, we can compute the anisotropic driving forces from ( 5), (6), and (7).Since the mobility anisotropy factor is a function of the driving forces, we need to solve for the mobility anisotropy factor to satisfy ( 5), (6), and ( 7) simultaneously.For more details, see Appendix.As for the driving field F, we employ the gradient of the quasi-Fermi potential ( n ).Therefore, F x = ∂ x n and F y = ∂ y n in ( 5) and (6).
While the derivation relies on the homogenous bulk condition, we can apply the derived model to a nonuniform simulation domain by assuming local one-to-one correspondence between the driving field and the carrier temperature.Note that this is related to the fundamental approximation of the DD model, which ignores the non-local energy flux terms in the carrier temperature equation [21], [22].This oneto-one relation further allows us to express carrier mobility as a function of the driving force rather than the carrier temperature.

B. Details on Mobility Model
Here, we summarize the details of the mobility models employed in this article for 4H-SiC materials.We adopt the Caughey-Thomas expression [23] for bulk mobility with the model parameters reported in [12].We use the Lombardi model [24] for the inversion layer mobility with the surface roughness mobility parameters from [10].We calibrate the acoustic phonon mobility parameters to match the experimental data [25].More specifically, the low-field electron mobility is given by where N tot is the total doping concentration, F ⊥ the normal electric field divided by 1 V/cm, and the mobility parameters µ max I , µ min , B, C, and δ are 1010, 5, 10 6 , 1.4 × 10 5 , and 5.82 × 10 14 cm 2 /Vs, respectively.
To take into account the anisotropy of the low-field mobility [4], [8], [12], we assume that µ max A = 1.2µ max I while the other parameters are assumed to be isotropic.
With the employed parameter set, the simulated effective electron mobility values agree well with the experimental data [25], as shown in Fig. 1.To further match the roll-off behavior in the weak inversion conditions, we need to include a mobility degradation term related to the Coulomb scattering due to interface charges [10], [11], which is not considered in this work.
For the high field saturation, we apply the Canali model [18] with the parameters from [12] as follows: where the saturation velocities perpendicular and parallel to c-axis, v sat I and v sat A are 2.2 × 10 7 and 1.32 × 10 7 cm/s, respectively (v sat A = 0.6v sat I ) [12].

C. Comparison to Other Approaches
In the commercial device simulator [20], the anisotropy of the driving force is not considered.Instead, the following Fig. 1.Comparison of the effective electron mobility versus effective field curves at room temperature for different substrate doping concentrations obtained from simulations (lines) with the experimental data (symbols) [25].
expressions are used for the driving force when the anisotropic mobility model is activated [20]: or depending on the user option.
Here, we describe some qualitative differences between our model, ( 13) and ( 14) in the high field saturation condition where we can assume that µ A|I ≈ v sat A|I /F A|I .When a sufficiently large field F = (F x , 0) is applied along the c-axis, (13) would give the electron velocity in the saturation limit as follows: which is inconsistent with the fact that the saturation velocity along the c-axis should be v sat A .On the other hand, ( 6) and ( 14) would give the consistent saturation behavior as follows: All three options give the correct saturation behavior when the field is applied perpendicular to the c-axis ( v y ≈ v sat I ).Therefore, we can expect that (13) would give incorrect results when the current flows along the c-axis, which is common as SiC power devices typically have a drift region along the c-axis.Also, note that our model would give r n ≈ v sat I /v sat A 2 as follows: whereas ( 13) and ( 14) would give r n ≈ v sat I /v sat A in the high field saturation condition.
In Fig. 2, we compare the magnitude of the drift velocity as a function of the angle between the driving field and the c-axis depending on the driving force options.In the calculation, the magnitude of the driving field is set to a sufficiently large value (3 MV/cm) such that the drift velocity is near the saturation velocity.We can see that the drift velocity is almost insensitive to the field direction when ( 13) is employed.On the other hand, the present anisotropic driving force model and ( 14) give consistent limiting values when θ = 0 and 90 • while they give slightly different results when the angle is in between these values.

D. Anisotropic Impact Ionization Model
For the given current density vector J = J x , J y and the driving forces F I and F A , we compute the anisotropic impact ionization coefficient as follows: with where a I = 2.10 × 10 7 /cm, a A = 1.76 × 10 8 /cm, b I = 1.70 × 10 7 V/cm, and b A = 3.30 × 10 7 V/cm for 4H-SiC [9].Note that the present model performs the log-linear interpolation of the impact ionization coefficients in the xand y-directions based on the direction of the current density vector, which is different from [20], where the individual model parameters are linearly interpolated based on the direction of the current density vector.Fig. 3 compares the impact ionization coefficient as a function of the angle between the driving field and the c-axis depending on the driving force options.We also plot the coefficient from the Hatakeyama model [9].For the cases of ( 13) and ( 14), we interpolate the individual parameters Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.according to [20].Our model and the Hatakeyama model give very similar results, whereas (13) largely underestimates the coefficient along the c-axis.Finally, (14) gives the correct limiting values when θ = 0 and 90 • , but it gives the maximum coefficient when θ = 77 • .

III. SIMULATION RESULTS
As an application example, we simulate the I D -V G and the breakdown characteristics of a 4H-SiC vertical power MOS-FET structure, as shown in Fig. 4. We implemented the present anisotropic driving force model as well as ( 13) and ( 14) in our in-house device simulator [26] and performed DD simulations at room temperature (T = 300 K) with the mobility and impact ionization models described in Section II.In addition, we consider relevant models such as the incomplete ionization of dopants [5], Fermi-Dirac statistics, and the SRH and Auger recombination models [2].We ignore the self-heating effect to focus on the influence of the driving force options, although it affects power device characteristics significantly.
Figs. 5 and 6 show the simulated MOSFET I D -V G and breakdown characteristics depending on the driving force options.Our model and (14) give similar results, whereas (13) overestimates the drain current for the I D -V G curve and the breakdown voltage because (13) underestimates the driving force in the drift region.To investigate the differences between the driving force models more closely, we plot the electron velocity, the mobility anisotropy factor, and the electron impact ionization rate along the x-direction at y = 0 µm when V D = 800 V and V G = 0 V in Figs.7-9.While the electron velocity obtained from the present model and ( 14) correctly follows v sat A , the velocity from (13) follows v sat I .In addition, our model gives the mobility anisotropy factor close to v sat I /v sat A 2 , whereas (13) and ( 14) predict the factor close to v sat I /v sat A .These results are consistent with the discussions in Section II-C.Since (13) underestimates the magnitude of the driving force in the drift region by a factor of 0.6, it predicts much smaller impact ionization rates in the drift region, as shown in Fig. 9, which explains why (13) overestimates the breakdown voltage.
Finally, we compare the convergence properties of the three options.To obtain the I D -V G curves in Fig. 5, all three models require 57 bias steps with similar convergence rates.To obtain the breakdown characteristics in Fig. 6, the present model requires 224 bias steps, including 98 failed ones.On the other  hand, (13) requires 180 bias steps (including 76 failed steps), and ( 14) requires 223 bias steps (including 97 failed steps).Therefore, we can conclude that (13) gives slightly better convergence rates while the present model and ( 14) show similar convergence behaviors.

IV. CONCLUSION
We presented a new anisotropic driving force model suitable for the TCAD simulation of SiC devices derived from the underlying relationship between mobility, driving force, and carrier heating.Based on the anisotropic driving force model, we developed the anisotropic mobility and impact ionization models for the DD simulation.We compared the present model with the existing models in the commercial simulator and highlighted the differences.We simulated the 4H-SiC vertical power MOSFET device and demonstrated its applicability.APPENDIX COMPUTATION OF r n , µ A , AND µ I We solve the DD model on an unstructured grid using the control volume method.While we assemble the continuity equation, we compute µ A|I for each element node in parallel as the computation only involves local quantities.More specifically, we employ a 1-D root finding algorithm based on the Newton method to find r n by solving R(r n ) = 0 with Once r n is determined, we can obtain µ A|I .We also need to compute the derivative of µ A|I with respect to solution variables to ensure quadratic convergence.For example, the derivative of µ A with respect to a solution variable X can be obtained from ∂r n ∂ X with

Fig. 2 .
Fig. 2. Comparison of the magnitude of drift velocity as a function of the angle between the driving field and the c-axis depending on the driving force options.

Fig. 3 .
Fig. 3. Comparison of the impact ionization coefficient (α) as a function of the angle between the driving field and the c-axis depending on the driving force options.Also plotted is the coefficient from the Hatakeyama model [9].

Fig. 5 .
Fig. 5. Comparison of the simulated drain current versus gate voltage of the 4H-SiC MOSFET depending on the driving force options.

Fig. 6 .
Fig. 6.Comparison of the simulated breakdown characteristics of the 4H-SiC MOSFET depending on the driving force options.For the breakdown simulation, we set V G = 0 V and ramp the drain bias up to 800 V.Then, we ramp up the drain current to 10 −2 A/µm.

Fig. 7 .Fig. 8 .
Fig. 7. Comparison of the electron velocity along the x-direction at y = 0 µm when V D = 800 V and V G = 0 V.

Fig. 9 .
Fig. 9. Comparison of the electron impact ionization rate along the x-direction at y = 0 µm when V D = 800 V and V G = 0 V.