Acceleration of solving drift-diffusion equations enabled by estimation of initial value at nonequilibrium

1.   Introduction

Semiconductor device simulation plays a crucial role in fostering the next-generation device technology development, and is therefore widely used in the industry and academia ^[1,2,3]. Despite the recent trend of device miniaturization that requires simulators to account for quantum transport effects, many devices with larger dimensions (at the scale of 1 $\mu m$) and with more complex geometries are being designed and implemented for various applications. Among the semiclassical approaches developed for modeling charge carrier transport, the drift-diffusion (DD) equations are among the most popular ones because of its simplicity while being capable of explaining many essential characteristics of semiconductor devices. The nonlinear iterative approach, composed of inner iterations and outer iterations, is typically employed to solve the strongly coupled nonlinear DD equations of the semiconductor devices ^[4,5,6]. Each inner iteration solves a linearized version of the entire nonlinear algebraic system based on numerical methods, such as finite element method (FEM) and finite volume method (FVM) ^[7,8], consuming a lot of time to achieve a stable convergence. The selection of outer iteration method will affect the number of times the inner iterations are solved. Therefore, a high-speed outer iteration method is urgently required to reduce the iteration number and accelerate the simulation process of semiconductor devices. For the outer iteration method, an effective initial value is the key-point to improve simulation efficiency without divergence in the iteration process.

In the classical method, the self-consistent solution is obtained through bias ramping process in the outer iterations. The previous solution is generally adopted as the initial guess or as the basis for extrapolating the initial value to solve the equations under current bias voltage on devices. The bias ramping process increases bias from the equilibrium state and repeats a series of outer iterations, which is a time-consuming computational process ^[2]. By contrast, obtaining an effective initial value directly from a nonequilibrium state at the target bias voltage can bypass the tedious bias ramping process. However, it is a challenge to acquire such initial value due to the nonlinearity of the DD equations.

Nowadays, the prediction of initial value in semiconductor simulation has been approached through two new methods, consisting of the neural networks method ^[9] and analytical models method ^[10,11]. The neural networks method could efficiently generate initial guesses when provided with sufficient training data and device parameters. However, the process of preparing the training dataset is often time-consuming. Compared to the neural networks method that requires a time-consuming training process, physics-based analytical models can obtain a more unconstrained initial guess of device by solving simplified coupled nonlinear partial differential equations that describe the physics of the semiconductor devices. A straightforward initial guess, including both potential and carrier density distribution for simulating the P-N junction, which is a simple device composed of two semiconductor materials with different doping types, p-type and n-type, was proposed to accelerate device simulation based on physical perspective ^[10]. Another analytical models method based on the compact charge model was proposed to predict an initial guess for gate-all-around metal-oxide-semiconductor field-effect transistors ^[11]. Instead of the compact charge model derived from the density-gradient equation, the potential analytical model based on the Poisson equation is a better way to construct a new initial guess. Since the range of potential variation is relatively smaller compared to the carrier concentration, it is less likely to cause a divergence problem, and the carrier concentration distribution could be numerically obtained in our analytical models method automatically.

In this paper, a novel estimation of initial value at nonequilibrium for semiconductor device simulation, which caused by applied bias voltage, is proposed based on analytical models about electrical potential. To validate this method, a commonly encountered junction barrier Schottky (JBS) device has been simulated. Numerical experiments show that this method can significantly reduce the number of inner iterations required for convergence and possesses strong robustness in this application scenario. To extend this method to handle more complex devices, one can start by establishing a more extensive unified analytical model of electrical potential or even combining some neural networks methods, although it may bring additional training computational costs.

The organization of this article is as follows. In Section 2, the mathematical model employed in semiconductor device simulations and several types of boundary conditions is introduced. Section 3 discusses the numerical methods regarding the finite element discretization and iterative schemes. The potential analytical model of JBS devices before punch through is proposed in Section 4. Finally, numerical simulation results of the two-dimensional JBS device are presented in Section 5 to evaluate the effectiveness of this method.

2.   Mathematical model

2.1.   The DD equations

The DD equations, simplified from the Boltzmann transport equation, is a semi-classical model for the mathematical description and numerical simulation for semiconductor devices. It expresses the potential and the carrier concentration relationships with the Poisson equation and continuity equations respectively. The equations read as follows ^[12]:

with

where $\psi$ is the electrostatic potential, $\epsilon$[CV$^{-1}$cm$^{-1}$] is the dielectric constant, $q$[C] is the fundamental electron charge, $n$ and $p$ are the electron and hole concentrations inside the semiconductor, and $C = N_D-N_A$ [cm$^{-3}$] is the doping profile. $C$ is assumed to be a given datum of the problem in terms of the donor and acceptor concentrations $N_D$ and $N_A$. Source terms $R_n, R_p$ [cm$^{-3} s^{-1}$] can be interpreted as the net recombination/generation rate of carriers per unit time and volume. In addition, assuming the temperature $T$[K] of the crystal is constant, carrier mobilities $\mu_n, \mu_p$ [cm$^2$ V$^{-1}$ s$^{-1}$] and carrier diffusion coefficients $D_n, D_p$ [cm$^2$ s$^{-1}$] follow Einstein's relations,

where $k_B$ [VCK$^{-1}$] is the Boltzmann constant. The carrier mobilities are treated as constant and calculated with the low field mobility model ^[2], Eq (2.4):

where $\mu_{n, 300}$, $\mu_{p, 300}$, $\gamma_n$, and $\gamma_p$ are low-field mobility parameters.

2.2.   Pretreatment of DD equations

First, the dimensionless electrostatic potential is treated with the following scaling for computational simplicity:

By taking Eqs (2.3) and (2.5) into DD Eqs (2.1) and (2.2), the scaled DD equations could be summarized as follows:

Second, the following scaled Slotboom variables in Eq (2.7) are introduced ^[13], which has been proven to be very useful and was adapted by many researchers ^[14,15]:

aiming at eliminating the cross terms, which is the advective terms in the flux density and transforming the continuity equations in Eq (2.6) into a set of self-adjoint second-order elliptic partial differential equations with exponential coefficients $exp⁡(-\psi)$ and $exp⁡(\psi)$:

2.3.   Generation-recombination terms of continuity equations

Carrier generation-recombination is a process in which the semiconductor material attempts to return to equilibrium after being disturbed from it. In this study, the generation-recombination mechanisms mainly include Shockley-Read-Hall (SRH) and Auger recombination ^[2]. The SRH model was introduced in 1952 ^[12,16] to describe the statistics of the recombination and generation of holes and electrons in semiconductors occurring through the mechanism of trapping. This model is an important ingredient of simulation models for semiconductor devices, especially the wide bandgap semiconductors such as gallium nitride (GaN) devices discussed in this paper. The net SRH recombination rate is written as

where $n_{ie}$ is the effective intrinsic carrier concentration of the semiconductor, $E_t$ is the trap energy level involved, $E_i$ is the intrinsic Fermi level, and $\tau_n$ and $\tau_p$ are the electron and hole lifetimes, which are related to the doping concentration:

where $\tau_{n0}$, $\tau_{p0}$, $N_n^{SRH}$, and $N_p^{SRH}$ are appropriate constants relating to materials.

For simplicity, in this paper, $E_t$ is set as $E_i$, which corresponds to the most efficient recombination center.

Auger generation-recombination mechanisms occurs through a three-particle transition whereby a mobile carrier is either captured or emitted ^[4]. Auger recombination is commonly modeled by the expression:

where $Aug_n$ and $Aug_p$ are appropriate constants relating to materials.

2.4.   Boundary conditions

Multiple boundary conditions are presented for kinds of semiconductor devices. In this paper, two main types of boundary conditions are taken into account: The nonhomogeneous Dirichlet conditions for ideal Ohmic contacts and Schottky contacts, and the homogeneous Neumann conditions for boundaries without contacts.

Ohmic contacts: Denoted with $\Gamma_O$, on which external voltages $V_{ext}$ are applied to electrically drive the device. The boundary conditions for the electrostatic potential $\psi$ and concentrations of carriers $n, p$ are all Dirichlet conditions:

Shottky contacts: Denoted with $\Gamma_S$, on which the applied external voltages are consistent with those on the Ohmic contact in JBS. According to the theory of hot electron emission, the boundary conditions are:

where $\chi$ is the electron affinity of the semiconductor material, $E_g$ is the bandgap, $N_c$ is the conduction band density of states, $N_v$ is the valence band density of states, and $\Phi_B$ is the barrier height at the metal-semiconductor interface in $eV$. The parameter of the working function, $W_f$, is used to specify a Schottky contact, which is defined as:

The above equation is formally similar to the metal work function. However, the barrier height is not strictly equal to the difference between the electron affinity and the metal work function, owing to the presence of interface states. In this paper, the parameter $W_f$ is set to $5.3938\; eV$ to have a Schottky barrier of $1.0838\; eV$.

The specific relevant parameters used in this paper are listed in Table 1.

Boundaries without contacts: All the external boundaries of the device without contacts, are denoted with $\Gamma_N$. These boundaries are treated as ideal Neumann boundary conditions:

3.   Numerical framework

Approximate solutions of DD equations are computed using numerical techniques as it is impractical to get analytical solutions. Various discretization methods such as finite difference method (FDM), FVM, and FEM are available to obtain numerical solutions. Many modern devices involve geometrically intricate structures. Therefore, FVM and FEM, which allow for unstructured meshes, have drawn more attention in recent years ^[18,19,20].

3.1.   FEM discretization for the DD equations

Let $\Omega \subset \mathbb{R}^2$ denote the whole semiconductor domain, $H^1 (\Omega) = \{u:\Omega\rightarrow \mathbb{R}|u, \frac{\partial u}{\partial x}, \frac{\partial u}{\partial y}\in L_2(\Omega) \}$ be the Sobolev space of weakly differentiable functions, and $H_D^1 (\Omega) = \{v\in H^1 (\Omega)|v = 0\ on\ \Gamma_O\cup \Gamma_S \}$. The variational formulas of the continuity equations in Eq (2.8) are used to find $\Phi_n\in H_D^1(\Omega)$ and $\Phi_p\in H_D^1(\Omega)$, satisfying

Let us assume that $\Omega$ is polygonal (as is often the case in practice), so that the domain could be tiled with a set of triangles $\bigtriangleup_k, k = 1, \cdots, K$, defining a triangulation $T_h$. The points where triangle vertices meet are called nodes, denoted as $q_i, i = 1, \cdots, N$. The test function $v$ is chosen in the piecewise linear finite element space $V_h\subset H_D^1(\Omega)$, which is denoted as $v_h$. The Slotboom variables $\Phi_n$ and $\Phi_p$ are discretized by $\Phi_{nh} = \varSigma_i\Phi_{nh}(q_i)\phi_i$ and $\Phi_{ph} = \varSigma_i\Phi_{ph} (q_i)\phi_i$, respectively. In these equations, $\phi_i$ denotes the linear Lagrangian basis function at $q_i$, then the discrete forms of the left-hand side of Eq (3.1) become,

where $\psi_*$ is the corresponding potential values at Gaussian points within a triangular region, extracted by interpolating the potential values at the vertices of $\bigtriangleup_k$.

For the Poisson's equation in Eq (2.8), by treating with Slotboom variables, the nonlinear equation in terms of the difference $\delta(x)$ between the available solution $\psi$ and the exact solution ^[21] can be expanded to

By adopting operations similar to those in the classical Newton method, neglecting terms of second and higher orders, and substituting Eq (3.3) into Eq (3.2), the following linear differential equation for $\delta$ is obtained:

The discrete variational form of the above equation is given by:

The complete coefficient matrices of Eqs (3.2) and (3.5) are constructed by the Galerkin finite element method ^[22].

3.2.   Outer iteration schemes

3.2.1.   Outer iteration scheme based on the analytical models

The process of constructing the potential analytical model-based iteration is given as follows. First, the analytical 2D potential distribution $\psi_{analytical}$ is obtained through the analytical model under a certain reverse bias voltage, $V_r$, for a specified doping concentration and device size, then the transitional Slotboom variables $\Phi_{n, trans}, \Phi_{p, trans}$ are obtained by numerically solving the decoupled continuity equations without generation terms using the potential distribution obtained in the previous step. With one step of Gummel iteration considering the generation terms in Eqs (2.9)–(2.11), the initial guess could be achieved. The combined $\psi_{initial}$, $\Phi_{n, initial}$, and $\Phi_{p, initial}$ are used as the initial guess under this certain reverse bias. Finally, the classical outer iteration scheme is employed to further iteratively solve the DD equations. The overall iteration procedure is shown in the flowchart in Figure 1.

3.2.2.   Classical outer iteration schemes for solving the DD equations

In classical nonlinear iterative schemes, the DD equations are mainly solved by the inner-outer iterative methods ^[3,5]. Among them, the Gummel or Newton iteration is used as the outer iteration, and the conventional iterative methods are used as the inner iteration for solving the linearized systems. In order to better compare our method with classical nonlinear iterative methods, three outer iterative schemes based on Gummel or Newton iterations were adopted here.

Scheme $\mathbf{G_1}$: This scheme is based on the Gummel iteration, and the flowchart is shown in Figure 2(a). Unlike our method, the classical outer iteration requires obtaining the initial guess from previous data. In Scheme $G_1$, the currently calculated solution is used as the initial guess at the next bias. Scheme $G_1$ adopts a simple damping method, which truncates corrections that exceed a maximum allowable magnitude. Scheme $G_1$ also limits the number of linearized Poisson solutions per Gummel iteration to one, which leads to under-relaxation of the potential update. This 'single-Poisson' solution mode extends the usefulness of Gummel's method to higher currents.

Scheme $\mathbf{G_2}$: This scheme is also based on the Gummel iteration, but the difference from Scheme $G_1$ is that the Poisson equation is solved repeatedly until convergence is reached as shown in Figure 2(b). This is done because the nonlinear Poisson solver acts as a 'preconditioner' of Gummel's fixed point iteration. A more commonly used initial guess strategy, called linear extrapolation (LE), is adopted in Scheme $G_2$. In Scheme $G_2$, an initial guess for a new bias is extrapolated from two previous solutions, that is

in which $\psi_{i+2, initial}$, $\Phi_{n, i+2, initial}$, and $\Phi_{p, i+2, initial}$ are initial guesses for the new bias and the other variables with subscripts $i$ and $i+1$ are the convergent solutions of the two previous iterations, with $V_i$ and $V_{i+1}$ representing their corresponding bias voltages.

In Scheme $G_2$, an exponential decline search damping is adopted by introducing variable $k$ and a 'residual' function $g(k)$. Derive the following equations from Eq (3.4):

Based on Eqs (3.7) and (3.8), the 'residual' function $g(k)$ is expressed as Eq (3.9):

in which $\psi_i$ represents the initial guess and $\delta \psi_i$ represents the correction for $\psi_i$. A particular $k_i$ is accepted if

in which, $\alpha$ is a small positive parameter, set as $10^{-4}$.

The damping process starts from $g(1)$. If $g(1)$ is sufficiently small, then the full Newton step is accepted, otherwise $k_i$ will be further reduced to $k_{i+1}$ by $k_{i+1} = k_i/2$, then test Eq (3.10) again.

Scheme $\mathbf{N}$: This scheme is based on the Newton iteration. The Newton-like method usually has the local quadratic convergency property if a suitable initial guess is given. According to Eq (2.8), let

In each Newton iteration, the following Newton equations based on Eq (3.11) are solved:

where

The same initial guess strategy, LE in Eq (3.6), as in Scheme $G_2$, is adopted in Scheme $N$. The damping method is based on Eqs (3.12) and (3.13), but it has been scaled to ensure that the corrections of different variables are within the same magnitude. So, Eq (3.12) is transformed into

in which scale factor $A$ is calculated by $max(1, \psi_{mmax}^{i+1})$, where $mmax$ is the node, at which $|\psi_m^{i+1}-\psi_m^{i}|$ has its maximum value, and $m$ is the node identifier. For scale factor $B$ and $C$, they are calculated by $max(\Phi_n^i)$ and $max(\Phi_p^i)$ respectively. By representing Eq (3.14) in the form of $\mathbf{M} \cdot \boldsymbol{\delta} = -\mathbf{F}$, the 'residual' function for the Newton iteration would be expressed as Eq (3.15):

and the remaining processing of the damping method is consistent with Scheme $G_2$.

4.   Analytical model about potential distribution of JBS with reverse voltage

The cross-sectional schematic diagram of the JBS device is shown in Figure 3, where the width and junction depth of the P+ implantation region are denoted as $w$ and $X_j$, respectively. The three boundary conditions described by Eqs (2.12)–(2.15) are represented by blue, red, and orange boundary lines in The drift region with a depth of $W_D$, which is a section designed to handle high voltages with low implantaion, is divided into two parts: the grid area and the yellow area. The grid region between the two P+ regions is part of the un-implanted n-drift region, with a width of $S$ and a depth of $d$. This grid region located beneath the anode Schottky contact is the most critical area of the JBS device. The parameter $d$ represents the maximum depth influenced by the P-N junctions on both sides. The drift region outside this rectangular area is considered unaffected by the lateral P-N junctions and has the same potential variation rate as the 1D Schottky junction diffusion region.

Setting the anode as the zero-potential reference, the potential distribution at a reverse bias voltage $V_r$ ($V_r < \frac{qN_dW_D^2}{2\varepsilon})$ is derived under several assumptions, which can be listed as follows:

$\bullet A_1$: The change rate of the longitudinal electric field is much greater than that of the transverse field in this grid region.

$\bullet A_2$: The potential along the P-N interface is zero.

$\bullet A_3$: The implantation regions on both sides of JBS only affect the area where $y < d$, while the change rate of the electric field in the region where $y > d$ is the same as that of a Schottky barrier diode (SBD).

With the boundary conditions as in Section 2.4, the potential distribution in the grid area can be concluded as follows ^[23,24]:

The expressions for $k_1$ and $k_2$ are provided in Eq (4.4).

The drift region outside the grid shaped area can be treated under the assumption that this region is not affected by the transverse PN junction and complete depletion. Also, the continuity of the potential needs to be maintained. So, the approach is to use a quadratic function with a second-order derivative of $-\frac{qN_d}{2\epsilon}$ to fit the potential at the upper and lower ends. The relevant formulas are concluded as follows,

The complete expressions for the coefficient $k_3, k_{30}, k_4$, and $k_{40}$ are provided in Eq (4.4).

The potential distribution on the remaining region, such as drift regions that only exhibit one-dimensional effects and heavily doped regions, can be figured out with a physical explanation ^[2]. The $n^+$ and $p^+$ doping regions can be approximated to have uniformly distributed quasi-Fermi levels $E_{F_n}^N$ and $E_{F_p}^P$. Therefore, the potentials in $n^+$ and $p^+$ doping region can be denoted as constant $\psi_N$ and $\psi_P$, respectively, given by

in which, $V_d$, which could be calculated by Eq (4.1), represents the voltage at a distance of $d$ from the Schottky contact at the horizontal boundary of the grid area, $V_{ohm}\approx V_r$ represents the voltage at a distance of $X_j$ from the Ohmic contact, and $V_{cathode}$ represents the voltage of the cathode voltage at the bottom.

So, the Eqs (4.1)–(4.3) form the analytical model about potential, $\psi_{analytical}$.

5.   Numerical results

The parameters for simulating the GaN JBS regarding the dimensions and doping concentration are presented herein. The substrate has an n-type doping concentration of $1\times10^{19}\; cm^{-3}$. The drift region, with a thickness of $1.5\; \mu m$, exhibits an n-type background carrier concentration of $2\times10^{16}\; cm^{-3}$. Additionally, the $p^+$ implantation regions have a width of $1.4\; \mu m$ and a depth of $0.2\; \mu m$, featuring a p-type concentration of $3\times10^{17}\; cm^{-3}$. The channel region has a sub-micron dimension of $0.6\; \mu m$ wide and $0.2\; \mu m$ deep, contributing to the realization of superior JBS rectifier characteristics ^[25].

In this test, the mesh consists of irregular triangular elements generated based on predefined density distribution functions, while adhering to Delaunay triangulation. Figure 4 presents a comparison between the initial guess obtained from the initial iteration scheme and the resulting converged solution. The experiments are conducted on a grid scale consisting of 10398 points, referred to as the $M_2$ grid scale. Since the analytical model is established before the device reaches punch-through voltage, the following reverse bias voltages are considered: $-5, -10$, and $-15\; V$. Figure 4(a)–(c) illustrates the 2D results of potential, electron concentration, and hole concentration for $V_r = -10\; V$. In these figures, the red circle represents the initial value distribution and the grid lines depict the final solution. To facilitate display, the number of initial guesses has been diluted by a factor of 10. As depicted in these figures, the initial distribution of potential and carrier concentration obtained through the improved initial guess strategy serves as an effective approximation of the final solution. Only the minority carrier concentration, specifically the hole concentration, exhibits noticeable variations in the channel region, such as the area between the $p^+$ implantation regions. Figure 4(d)–(f) displays the extracted data along the centerline traversing the channel region ($x = 0\; \mu m$). When the reverse bias voltage increases from $-5$ to $-15\; V$, the 2D initial potential and hole concentration values closely align with the final solution in Figure 4(d), and Figure 4(e), demonstrating a significant match. As shown in Figure 4(f), the deviations are confined to the hole concentration range of 0$\sim$1$\; \mu m$, with hole concentrations below $1\times10^{-14}\; cm^{-3}$.

Figure 5 presents a 2D contour plot that illustrates the potential difference between the initial guess and the final results. The plot demonstrates that near the lateral P-N junction interface, the maximum potential deviation is approximately $0.88$, $1.17$, and $1.3\; V$, corresponding to $V_r$ values of $-5$, $-10$, and $-15\; V$, respectively. While certain regions exhibit some deviations, the majority of regions exhibit deviations within $0.2\; V$, providing an effective initial potential value for subsequent nonlinear iterations.

The iterative scheme utilized in this study, based on the analytical model, eliminates the need for a ramping process and will be compared to the classical iteration method that starts from the equilibrium state. For simplicity, the former method will be referred to as the $AMmethod$, while the latter will be referred to as the $ESmethod$. Tables 2–4, respectively, display the number of inner iterations required for these two methods in different outer itereation schemes (Scheme $G_1$, Scheme $G_2$ and Scheme $N$) to achieve the same level of converence. It should be noted that due to the consistent system scale processed in each iteration, the calculation time is similar in each iteration of the same outer iteration scheme.

In the $ESmethod$ using Scheme $G_1$, a fixed bias increment of $0.5\; V$ was used. The experiment data in Table 2 revealed that for the $ESmethod$ using Scheme $G_1$, with each $5\; V$ increment in reverse bias voltage, the number of iterations increased by approximately 250. As the reverse bias voltage increased, the $ESmethod$ using Scheme $G_1$ consistently required more iterations to achieve a self-consistent solution, while the $AMmethod$ never exceeded 10 iterations.

For the $ESmethod$ using Scheme $G_2$, several bias steps have been used, ranging from $-0.5$ to $-5\; V$. Due to the use of the LE as the initial guess strategy for the $ESmethod$ using Scheme $G_2$, in addition to the need to first calculate the equilibrium state in Scheme $G_1$, an additional solution at $V_r = -0.2\; V$ was calculated as the basis for extrapolating initial guesses, which requires an additional 23 steps of solving the Poisson equation and 4 steps of solving continuity equations. As shown in Table 3, for different reverse biases, the lowest number of iterations occurs with the bias step, $\mathbf{V}_{step}$, at the maximum step size. However, excessive bias step size may lead to divergence, making it a trade-off issue. Comparing the data under $\mathbf{V}_{step} = -0.5\; V$ in Table 3 with the data in Table 2, it can be found that Scheme $G_2$ reduces the number of times the Poisson's equation is solved by more than half compared to Scheme $G_1$, and the number of times the continuity equations is solved is reduced by more than four fifths.

For the $ESmethod$ using Scheme $N$, only bias step of $-0.5$ and $-1\; V$ are adopted, and a lager bias step size will lead to divergence. This is because the Newton method relies on more accurate initial guesses than the Gummel method, and the initial guess obtained by the initial guess strategy, LE, still cannot meet the needs of the converged solution by Newton's iteration in this experiment. As shown in Table 4, whether $V_{step} = -0.5$ or $-1\; V$, the number of iterations required for the $ESmethod$ to reach a converged solution under different reverse biases is more than ten times that of the $AMmethod$.

These results clearly demonstrate that the $ESmethod$ necessitates a significantly greater number of iterations to produce converged results comparable to those obtained using the $AMmethod$.

In addition, the convergence behavior of the $AMmethod$ using Scheme $G_1$ is depicted in Figure 6. In Figure 6(a), it is evident that the differences in root mean square (RMS) values of the potential update $\delta$ become significantly pronounced after the sixth iteration. To enhance clarity in the presentation, Figure 6(b) was created to display the RMS values associated with the number of iterations with different bias voltages, starting from the seventh iteration. Specifically, at the ninth iteration, the corresponding RMS values are $2.321\times10^{-9}, 1.976\times10^{-7}$, and $7.816\times10^{-7}$ for $V_r = -5, -10, -15\; V$, respectively. The convergence rate is faster for lower reverse bias voltages, indicating a relatively quicker attainment of the same RMS condition through nonlinear iterations. This can be attributed to the increased nonlinearity of the DD equations with higher bias voltages. In Figure 6(c), the RMS curves for different grid scales, namely, $M_1 = 4674$, $M_2 = 10398$ and $M_3 = 20881$, at three bias voltages almost overlap. This suggests that the grid scale size, under the same grid generation method, has minimal impact on the convergence rate of nonlinear iterations in semiconductor simulation. The convergence rate of the solving process primarily depends on the inherent complexity of the problem itself and the employed nonlinear iteration strategy.

6.   Conclusions

In conclusion, a novel estimation of initial value at nonequilibrium for solving the DD equations is presented based on analytical models about electrical potential. The method starts by deriving an analytical potential distribution, taking into account the device dimensions, shape, and electrical parameters. Subsequently, the transient carrier concentration distribution is obtained by solving the decoupled continuity equations without the generation terms. By taking generation terms into account, a complete initial value under the specified bias is achieved.

Our method was verified on JBS devices under reverse bias voltage and compared with classical semiconductor device simulation methods that start from equilibrium state. The results indicate that regardless if the outer iteration is a Gummel iteration or Newton iteration with damping, our method saves more than ten times the number of iterations compared to classical methods.

Since one step of our method needs to ignore the generation and combination terms in solving the continuity equation. The proposed method is more applicable under reverse bias voltage. Specifically, it is most effective when there is minimal recombination of charge carriers. In the cases where the generation and recombination effect is significant, it is recommended to use the Scheme $G_1$, which limits the number of linearized Poisson solutions per Gummel iteration to one, to do the outer iteration.

For grid area of arbitrary shape, which is caused by different shapes of $p^+$ implantation region in JBS, it is convenient to continue using the method in the paper by adjusting the parameter $d$ to match different grid shapes. Since $d$ is a quantity strongly correlated with dimensions, it will be obtained using a fitting approach with experimental data. For more complex situations, such as nonuniform doping or more complex physical shapes, which can cause great difficulties in constructing analytical models. In this situation, neural networks, which require some time for training, can be used to solve the problem of getting an appropriate initial potential distribution, and then combined it with our nonlinear iteration scheme for a converged solution.

The proposed method offers a substantial reduction in simulation time for semiconductor device simulations, such as reverse leakage and reverse breakdown voltage, as it eliminates the need for a bias ramping process. It also eliminates the need to build a complex charge model from the beginning. Instead, it solely relies on the analytical model about electrical potential, leading to a significant reduction in the number of iterations required to achieve convergence for the desired bias voltage.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was supported by ShanghaiTech University Startup Fund 2017F0203-000-14, the National Natural Science Foundation of China (Grant No. 52131303), Natural Science Foundation of Shanghai (Grant No. 22ZR1442300), and in part by CAS Strategic Science and Technology Program under Grant No. XDA18000000.

The authors gratefully thank Professor Qifeng Liao at the Visual & Data Intelligence Center of ShanghaiTech University for his helpful instructions on finite element discretization.

Conflict of interest

The authors declare there is no conflict of interest.