Solution of the Semiconductor-Device Equations by the Numerov Process

Article history: Received: 30 June, 2020 Accepted: 12 October, 2020 Online: 16 December, 2020


Introduction
A number of physical problems are described by second-order differential equations of the form (primes indicate derivatives with respect to x). When a uniform grid with grid size h is used, such a form can be tackled numerically, thanks to the absence of the first derivative, with an O(h 4 )-accurate discretization method generally known as Numerov Process (NP); this in contrast to O(h 2 ) of the finite-difference method (the derivation of the scheme is shown in the Appendix). Linearization of (1) yields − z = Q(x) z + P(x) .
Scope of this paper is to extend the applicability of NP to the system of equations that describe the transport of electric charge in solid-state devices. This paper is an extension of work originally presented at the SISPAD 2019 Conference [1]; besides the Appendix, more details about NP are given, e.g., in [2] and references therein.
In some cases, (2) can be reduced to equations that appear in the model for semiconductors; for instance, letting Q = 0, and naming z = u the electrostatic potential normalized to the thermal voltage k B T/q, and P = q/(ε k B T ) , with the charge per unit volume, transforms (2) into the Poisson equation, which is one of the equations of the standard semiconductor-device model (in the above definitions, symbols q, ε, k B , T denote the elementary charge, material permittivity, Boltzmann constant, and lattice temperature, respectively). If, instead, one lets P = 0, the form of the timeindependent Schrödinger equation is recovered, which enters the semiconductor-device model when nanometer-scale devices are considered; in this case, z = w is the time-independent wave function, whereas the form of Q reads Q = 2 m (E − V)/ 2 , where E = const is the total energy of the electron and V(x) the potential energy; in turn, m is the electron effective mass and the reduced Planck constant.
The elimination of the odd-order derivatives in the Taylor expansion (see Appendix), as required by NP, prescribes the use of a uniform grid; the better accuracy of the method largely compensates for this aspect. An application of NP to the analysis of a semiconductor device was shown in [2], where the coupled solution of the Schrödinger-Poisson system was dealt with; in that case, no manipulation of the equations was necessary because both have already the form (2).
Object of this paper is to show that, by a suitable choice of the unknowns, NP can be extended to the whole set of equations constituting the mathematical model of semiconductor devices. To this purpose it is convenient to briefly outline the recent progresses in modeling, following the semiconductor-device roadmap (now "International Roadmap for Devices and Systems -IRDS TM " [3]) and the corresponding evolution of the mathematical models.
The standard model of semiconductor devices (drift-diffusion model) is made of the Poisson equation and, for each type of mobile charges (electrons and holes), of the continuity and transport equations for the mobile-charge density; the transport equation has two terms: the first one is of the drift type, the second one of the diffusive type. For submicron devices the model is extended by adding another pair of continuity and transport equations, whose unknown is the energy density of the carriers (in this case, the model is termed hydrodynamic).
When the lateral extension of the device channel shrinks to the nanometer scale, the device can be described as a cylindrical structure. This makes it possible to separate the longitudinal and transversal directions; then, the calculation of the charge density is accomplished by solving the Schrödinger equation over each cross section, whereas the analysis along the longitudinal direction is still carried out in classical terms. An example is found in the solution of a cylindrical nanowire shown in [2]. Other examples of nanometersize structures made of compound materials are illustrated in [4]; a recent analysis of the Schrödinger-Poisson system of equations in InGaAs-based HEMTs is shown in [5].
The n + -n-n + structure considered in this paper is suitable for checking the extension of NP to the solution of the whole transport model of semiconductors. The structure is one of the benchmarks used in the literature (see, e.g., [6] and references therein): besides being one-dimensional, like a nanowire or a gate-all-around transistor, it has the advantage of being unipolar, namely, for its description it suffices to use only one type of carriers (electrons in this case), thus halving the number of continuity-transport pairs of equations to be solved.
The investigations on numerical methods for solving the semiconductor equations have gone in parallel with the strong development of the integrated-circuit technology of the last 50 years. At present, the most successful solution techniques for the driftdiffusion or hydrodynamic models are incorporated in huge deviceanalysis systems, whose update and maintenance is a concern of the semiconductor Companies in cooperation with specialized Software Houses. The advantage of NP is that it is applicable to the existing solution techniques without excessive changes in the organization of the software, and with a limited increase in the computational cost; for this reason, also considering the improvement in accuracy, the application of NP to the semiconductor equation is worth considering. The present investigation about the applicability of NP to the numerical solution of the semiconductor transport equations considers the drift-diffusion model for the sake of simplicity, since the hydrodynamic model and the other models of higher order, derived from the Boltzmann Transport Equation, have the same structure.
The paper is organized as follows: the model equations are shown in Section 2 along with the transformation that gives them an NP-suitable form; the NP-based discretization in the onedimensional case is carried out in Section 3. A method for calculating the electric potential at each node independently of the coupling with the transport equation is outlined in Section 4, followed by a stability analysis of the iterative solution (Section 5). The results are shown in Section 6, and the extension of the method to the two-dimensional case, still over a uniform grid, is given in Section 7. Finally, the conclusions are drawn in Section 8, while a summary of the derivation of NP is provided in the Appendix.

Model Equations
The semiclassical, drift-diffusion equations describing the transport of carriers in solids is made of the Poisson equation and, for each type of charges (namely, electrons and holes), of a pair of equations made of the continuity and of the transport equations. Here a onedimensional device of the unipolar type is analyzed, whose carriers are electrons. This situation is found, e.g., in n + -n-n + structures like that shown in Figures 1 and 2, or in memory devices of the phasechange type (e.g., [7,8]); in these devices, the Poisson equation reads with n the electron concentration and N D (x) the concentration of ionized donors within the n + -n-n + structure (when phase-change materials are considered, N D is replaced with the concentration n T of empty traps). The continuity equation of the electrons, and the corresponding drift-diffusion equation for the electron-current density read, respectively, with U = U(n), D n being the net-recombination rate and the electron-diffusion coefficient, respectively (here a constant value of D n is considered). Combining (4) yields The Wronskian determinant of the homogeneous equation associated to (5) reads, by Abel's identity, W = const × exp(u) 0.
The system of equations to be solved is then made of (3) and (5), whose boundary conditions are the values of u and n at the two ends of the (finite) integration interval. If the discretization of the interval www.astesj.com produces N internal nodes, the boundary conditions are u 0 , u N+1 and n 0 , n N+1 , respectively; in the analysis of solid-state devices it is assumed that the contacts force a condition of electric neutrality, so that n 0 , n N+1 are such that P 0 = P N+1 = 0 (compare with (3)).
In a decoupled solution scheme, one assumes that n is given and solves (3) for u, then solves (5) for n using the form of u thus found; the procedure is iterated until convergence is reached. One notes that NP is not applicable to (5) as it stands, due to the presence of the first derivative; to give the equation a more suitable form one defines this transforming (5) into Transformation (6) does not consist in a reduction to a self-adjoint form; in fact, the expression of the latter would be transformation (6) is not either an "exponential fitting" of the type commonly adopted for solving the semiconductor equations: the exponential fitting would in fact yield where the current density J n is approximated with a different constant along each interval between two nodes (in the field of semiconductor-device modeling the scheme based on (9) is also known as Scharfetter-Gummel method [9]). The homogeneous equation corresponding to (7) reads g +c g = 0; in the iterative procedure by which the system made of (3) and (7) is solved within the finite integration interval, coefficient c is obtained from the electron concentration n and the normalized electrostatic potential u calculated at the previous iteration. As neither n nor u have poles, all points within the integration interval are ordinary.

Application of NP
In this Section, the Numerov Process is applied to (3) and (7) over a grid made of N, uniformly-spaced internal nodes with element size h. The Poisson equation (3) is transformed into the N × N algebraic system i = 1, 2, . . . , N, while (7) transforms into the N × N system The matrix of the algebraic system (10) is symmetric, whereas that of (11) is non symmetric due to terms c i−1 and c i+1 at the left hand side. Although in the latter system the unknowns are g i = n i exp(−u i /2), terms c i and s i still depend on the original unknown n; also, each g i depends, through exp(−u i /2), on the (arbitrary) zero of the electric potential; finally, due the presence of the exponentials of the electric potential, the matrix of the coefficients of (11) may become stiff. For the above reasons it is convenient to multiply both sides of (11) by exp(u i /2). In this way, the following replacements are carried out in (11): www.astesj.com 1416 The new form of (11), where the original unknown appears, is preferable because the exponents bearing differences like, e.g., u i − u i−1 , instead of a nodal potential alone, are expected to be relatively small. To complete the calculation it is necessary to express the derivative u that appears in (7); this is also done using NP, so that the accuracy of the scheme is preserved, and yields When i = 2, . . . N − 1, only the indices of internal nodes appear in (12); if, in contrast, i = 1 or i = N, then one boundary condition is present (u 0 and P 0 = 0 or, respectively, u N+1 and P N+1 = 0); the last possible case occurs when i = 0 or i = N + 1, in which the values of the electric potential and charge in the interior of the contacts appear in (12); remembering that the contacts are equipotential and neutral, it is u −1 = u 0 , P −1 = 0 and u N+2 = u N+1 , P N+2 = 0.

Decoupling the Model Equations
The model equations are solved by iterating between (10) and (11). In general, equations like these are tackled with algebraic solvers able to provide the values of the unknowns in a given sequence; for instance, considering the solution of (10) obtained through the A = L U decomposition, the ith nodal value of the normalized electric potential u is obtained only after calculating it at nodes 1 through i − 1 (or at nodes N through i + 1). Here a different approach is adopted, in which u i is obtained as soon as necessary, without the need of calculating its other nodal values. Such a result is achievable for the discretized form (10) of the Poisson equation by solving it with the scheme of [10, p. 769], that provides each nodal value u i independently of the others; using such a scheme after indicating with C i the right hand side of (10) yields where The differences that appear in (11) after carrying out the replacements outlined above, are easily evaluated from (13)(14)(15): In practice, this scheme decouples (10) from (11); it may appear that this result is due to the discretization: in fact, one can also eliminate the Poisson equation prior to the discretization, by recasting it in integral form [10, p. 781-784]. The calculation of the overall number of operations also shows that the method based on (13)(14)(15) requires fewer operations than the A = L U decomposition; specifically, for a matrix of order N the number of operations required by the A = L U decomposition is 6 (N − 1) multiplications and 3 (N − 1) additions, whereas the method based on (13-15) requires N − 1 multiplications and 4 (N − 1) additions.
Another difference between the standard discretization schemes like, e.g., finite differences, and the present one is the following: here, the discretized functions and their derivatives of all orders "belong" to the nodes; as a consequence, no hypothesis is necessary about the form of the discretized functions over each element. In the finite-difference scheme the functions and their derivatives of even order belong to the nodes, whereas the derivatives of odd order belong to the elements.

Stability
Observing that the equations of the model are non linear, it is necessary to implement an iterative solution. Specifically, considering for instance the equilibrium case, the normalized charge concentrations P i−1 , P i , P i+1 at the right hand side of (10) depend on the electric potential u; the dependence is exponential in a non-degenerate semiconductor, or via a Fermi integral in a degenerate one [10]. In both cases, the derivatives dP i /du are negative irrespective of the type of carriers (electrons or holes) that is considered: in fact, electrons (holes) contribute negatively (positively) to the charge density, and the electron (hole) concentration increases (decreases) when u increases. This reasoning also applies in the non-equilibrium case, because the functional dependence of the concentration on the electric potential is the same. In conclusion, linearizing (10) with respect to u adds weight to the main diagonal of the system matrix, resulting in an improved convergence rate; this behavior, typical of the linearized Poisson equation in semiconductors, is due to the dependence of the carrier concentrations on the electrostatic potential and is common to a number of discretization schemes. In contrast, as shown below, the properties of the algebraic system obtained by discretizing the continuity and transport equations with the scheme shown in this paper are quite different from those obtained with the standard scheme.
Coming now to (11), when the expressions of c i−1 , c i , c i+1 that appear in (7) are replaced in (11), one finds an algebraic system whose right hand side, in the ith row, is made of three summands: the structure of the first summand, A i = −g i−1 + 2 g i − g i+1 , is the same as that of Poisson's equation (10). The form of the other two terms is www.astesj.com    The coefficients defined by (16) are non negative, therefore they add weight to all diagonals of the system matrix (11); due to factor 10, the weight added to the main diagonal is generally dominant, unless node i is near an extremum of u. As for (17), the presence in it of the normalized charge density P, which may have either sign, makes the analysis more complicated; on the other hand, (17) is of order 2 in h, whereas, due to a cancellation that takes place when (12) is replaced in (16), B i is of the same order as A i , namely, order 0 in h. In essence, the structure of the algebraic system (11) deriving from the present discretization scheme of the continuity and transport equations is similar to that deriving from the Poisson equation. In contrast, the standard exponential-fitting scheme for the continuity and transport equations, based on (9), approximates U as a piecewise-constant function at each node, and linearizes u between each pair of nodes; with these premises, the scheme yields whose coefficients K k j depend on the normalized electric potential through the Bernoulli function B: Besides the approximations mentioned above, it is easily found that the main diagonal of the algebraic system (18) is not necessarily dominant [10, p. 781].

Results
As an example of application, the solution scheme based on NP is applied here to the n + -n-n + structure whose lateral section and dopant profiles are shown in Figures 1 and 2, respectively. The device length is L = 10 µm, the concentration of the light, constant dopant (black line in Figure 2) is N B = 10 15 cm −3 , while the two profiles (red lines) are given by N L = N 0 exp(−x 2 /x 2 0 ) and , with N 0 = 10 17 cm −3 and x 0 = 1.165 µm, respectively. The total dopant distribution to be used in (3) is N D (x) = N B + N L + N R . The device is uniform in the direction normal to the field; the equations of the models have been solved at different applied voltages V, ranging from equilibrium (0 V) to 1.6 V. The form of the electric field E is shown in Figure 3 for V = 0 V and V = 1 V.
Each solution was repeated on different grids, starting with a reference grid having N = 2000 nodes and successively decreasing it down to N = 100. Figure 4 compares the electric field E calculated using the 2000-node grid with that calculated using coarser grids. As all grids are obtained by refining the 100-node one, for each node of the latter there exists a node of the denser grids having the same abscissa. Consider, e.g., the V = 0 curve of Figure 4, for which the 500-node grid yields about 4 × 10 −9 ; to obtain this value one considers two vectors, namely, the electric field E c at each node of the 500-node grid and the electric field E d at the corresponding abscissae of the 2000-node grid; from this vector one then extracts the maximum relative difference max g |(E d − E c )/E d |, which is the strictest metric for the problem in hand. Figures 5 and 6 show the same type of comparison carried out for the electron concentration n and the electric potential ϕ, respectively. The worst case (about 9%) occurs for the electron concentration at the largest bias; this was to be expected, given the exponential dependence of the concentration on the electric potential.
Another part of the investigation aimed at determining the number of calls to the solver (i.e., iterations) that are necessary to bring the error below a prescribed value, given the number of grid nodes and the applied bias. The error is defined as where i is the node index, k is the iteration index, and q stands for ϕ or n. The number of calls is plotted in Figure 7 for a 1500-node grid, in Figure 8 for a 500-node grid and, finally, in Figure 9 for a 100-node grid, respectively; the bias values are the same in all cases. The curves show that the prescribed error is reached smoothly for all values of the applied bias, with no substantial difference from one grid to another; thus, the method provides a significant gain as far as the computational cost is concerned. www.astesj.com In [1], a different electronic device has been used to compare the error of the method presented here with that of the exponentialfitting one (the Scharfetter-Gummel method), in which the discretization scheme for the carrier concentration is based on (9). Like in the example above, after obtaining the reference solution over a dense grid (N = 2000), more solutions were calculated by making N to progressively decrease; after each solution, the maximum difference was calculated with respect to the reference solution. The improvement of the present method with respect to the standard one was found to be about one order of magnitude in all cases.

Extension to Two Dimensions
It may be argued that a method for solving one-dimensional problems has a limited range of applicability. This is not so in the modeling of solid-state devices; in fact, when dealing with transport problems in structures like, e.g., carbon-nanotube transistors or nanowires, that are essentially cylindrical objects whose diameter is in the nanometer range, the restriction to the one-dimensional case is not too severe a constraint. In fact, for such devices the problem is typically separated by decoupling the longitudinal coordinate (the channel) from the transversal ones [12]; as no charge transport takes place in the directions transversal to the channel, separating the problem yields, for each transversal section, a system made of the Poisson and Schrödinger equations, which are further separated by exploiting the radial symmetry of the device. As shown in Section 1, the one-dimensional Poisson and Schrödinger equations are solvable by NP. In turn, the transport along the channel is described by equations of the form discussed here, also amenable to the NP scheme.
Clearly the NP method would gain in flexibility from the extension to a non-uniform grid. This issue is outside the scope of this paper; the interested reader may refer to a method for extending NP to a variable stepsize, still in one dimension, based on the k-step Cowell method, which has been tested on the Schrödinger equation (see [11] and references therein).
Another direction for evolving the NP scheme is the extension to more than one dimension using uniform grids of the tensor-product type. Such an extension, again applied to the case of the Schrödinger equation, is described in [13]. It is of interest to ascertain whether the approach of [13] can be extended to the class of equations investigated here. To this purpose, one must first determine the multi-dimensional form of (7); the latter is obtained as follows, this yielding The multi-dimensional form (22) shows that the equations describing charge transport in a semiconductor can be reduced to a single second-order equation of the elliptic type.
The notation becomes rather awkward even if one limits the extension to the two-dimensional case; taking a uniform, tensorproduct grid, and using the matrix notation of [13], namely, a 11 g j−1 i+1 + a 12 g j i+1 + a 13 g j+1 i+1 + a 21 g j−1 i + · · · (23) + a 23 g j+1 i + · · · + a 33 g j+1 a 11 a 12 a 13 a 21 a 22 a 23 a 31 a 32 a 33 in which the lower indices of g j i , . . . refer to the x axis, and the upper ones refer to the y axis, the two-dimensional expansion yields an interpolation of the form Calculating the Laplacian of both sides of (24), eliminating [∇ 2 (∇ 2 g)] j i after neglecting the derivatives of the 6th-order and, finally, replacing (∇ 2 g) j i from (22) and defining the matrices eventually yields the two-dimensional generalization of (11): The result expressed by (26) shows that the approach of [13] can indeed be extended to the second-order equation of the general form; like in the one-dimensional case, no special assumption is necessary on the form of the discretized functions inside each elements or along its edges. www.astesj.com

Conclusions
It has been shown in this paper that the standard charge-transport model in solid-state devices, made of the continuity and transport equations for the charge, can be given a form suitable for the application of the NP scheme, using a uniform grid in one or two dimensions. With respect to the standard discretization schemes used in the analysis of semiconductor devices, the increase in the computational cost is modest; like in the case examined in [2], where the solution of the Schrödinger equation was takled, the NP approach has provided a simple tool for improving the solution of the equations.
In the field of semiconductor-device modeling, the results presented here are important from two points of view, already mentioned in Section 1: first, the NP approach, besides being able to significantly improve the solution, is applicable without excessive changes in the organization of the existing software; second, the drift-diffusion model on which NP has been tested is in fact the simplest level of description of solid-state devices; more sophisticated models exist, made of higher-order moments of the Boltzmann transport equation. On the other hand, such moments can always be grouped in pairs, each pair having the same structure: specifically, one equation of the pair is an even-order moment (order 2 k, with k = 0, 1, . . .) whose general form, with a suitable meaning of symbols, reads [10] − div S = C .
The next equation of the pair is the odd-order moment of order 2 k + 1, whose form is S = a grad σ + σ ∇b .
Comparing (27, 28) with (4) shows that all pairs of moments of the Boltzmann transport equation have the same structure. It follows that the method worked out in this paper applies to any order of transport models; also, considering that in the dynamic case the term C in (27) embeds the time derivative ∂σ/∂t, the method is applicable also in the dynamic operating mode of solid-state devices.