Discretization and Numerical Solution

1. Theory

• Definition: Device TCAD turns the drift-diffusion equations into algebraic equations on a mesh, most commonly by the box (Voronoi finite-volume) method with Scharfetter-Gummel fluxes, and solves the resulting non-linear system by Gummel or Newton iteration (Farrell; Chen and Bagci; Jerome and Kerkhoven).
• The flux problem: Carrier densities vary exponentially with potential, so linear interpolation between nodes gives spurious currents unless the mesh is extremely fine; the Scharfetter-Gummel flux integrates the current equation exactly along each edge (Farrell; computed).
• Non-linear solution: Gummel’s method solves Poisson and the two continuity equations in turn; Newton’s method solves them simultaneously with the full Jacobian (Chen and Bagci; Jerome and Kerkhoven).
• Accuracy: Errors fall with mesh size at a rate set by the scheme, about second order for SG currents in practice (Farrell; Chen and Bagci).

2. Working Principle

• Control volumes: Each node owns a Voronoi cell, and fluxes cross its faces, so charge and current are conserved exactly at the discrete level (Farrell; background knowledge for exact conservation).
• Exponential fitting: Assuming constant current and field along an edge gives an exponential density profile and the Bernoulli-weighted flux (Farrell).
• Decoupling: Gummel freezes the quasi-Fermi potentials while solving Poisson, then updates carriers with the new field, which converges well when coupling between the equations is weak (Chen and Bagci; background knowledge for the coupling condition).
• Full coupling: Newton linearizes all three equations together, converging quadratically near the solution, with damping to stay robust from poor initial guesses (Jerome and Kerkhoven; background knowledge for quadratic convergence).

3. Applications

• Commercial and open TCAD: Box-method SG discretization in DEVSIM and similar tools (Drift-Diffusion Device Equations page; Farrell).
• Heterostructures and degenerate devices: Generalized SG schemes for Fermi-Dirac statistics in lasers, LEDs and highly doped contacts (Farrell).
• Large 3-D problems: Iterative linear solvers with preconditioners for millions of unknowns (Chen and Bagci; Jerome and Kerkhoven).
• Alternative discretizations: Discontinuous Galerkin methods reaching the same accuracy on coarser meshes (Chen and Bagci).

Finite-Volume Discretization

Name	Formula	Related Formulae	Variables/Constants
Box (Voronoi finite-volume) discretization	${\sum }_{L\in N(K)}{\sigma }_{KL}ϵ\left({\psi }_K-{\psi }_L\right)=|K|q\left(p_K-n_K+C_K\right)$	${\sum }_{L\in N(K)}{\sigma }_{KL}{j}_{n,KL}=|K|q{R}_K,{\sigma }_{KL}=\frac{|\partial K\cap \partial L|}{|x_K-x_L|}$	$\begin{array}{r}K,L:neighbouringmeshnodes,N(K):\\ neighboursofK\end{array}$ $\begin{array}{r}{\sigma }_{KL}:edgecoefficient(sharedVoronoiface\\ areaoveredgelength),|K|:Voronoicell\\ volume\end{array}$ $\psi :potential$ $C:netdoping$ $j_{n,KL}:electronfluxonedgeKL$ $\begin{array}{r}R:recombination;Farrell:Voronoicells\\ withoneunknownpernode,integrationof\\ thevanRoosbroeckequationsovereach\\ cellwiththedivergencetheorem,giving\\ sumsoverneighbourswithedge\\ coefficientsandfluxes(edge-coefficient\\ definitionandsignarrangement,\\ backgroundknowledge);DEVSIMusesthe\\ sameedge-fluxassembly(Drift-Diffusion\\ DeviceEquationspage)\end{array}$ $(nocomputedvalue)$
Central-difference flux and its failure	$j_{KL}^{CD}=\frac{qD}{h}\left[\left(n_L-n_K\right)-\frac{n_K+n_L}{2}\frac{{\psi }_L-{\psi }_K}{U_T}\right]$	$\left|\frac{{\psi }_L-{\psi }_K}{U_T}\right|<2(\text{needed for stability})$	$\begin{array}{r}{j}^{CD}:naivecentral-differenceelectronflux\\ (signconventionwithfluxtowardK\\ positive,backgroundknowledge)\end{array}$ $D:diffusivity$ $h:edgelength$ $\begin{array}{r}{U}_T:thermalvoltage;thepotentialstepper\\ edgemuststaybelowabout2U_{T}orthe\\ schemeoscillatesanddensitiescanturn\\ negative(backgroundknowledge)\end{array}$ $\begin{array}{r}(computed:atequilibriumwitha5U_{T}step,\\ n_L=e^5{n}_K=148n_K,thecentralfluxisnot\\ zerobut-226inunitsofqD{n}_K/h,1.5times\\ thediffusionterm)\end{array}$

Scharfetter-Gummel Fluxes

Name	Formula	Related Formulae	Variables/Constants
Scharfetter-Gummel flux	$j_{n,KL}=\frac{q{\mu }_n{U}_T}{h}\left[B\left(\frac{{\psi }_L-{\psi }_K}{U_T}\right)n_L-B\left(\frac{{\psi }_K-{\psi }_L}{U_T}\right)n_K\right]$	$B\left(x\right)=\frac{x}{e^x-1},B\left(0\right)=1,B\left(-x\right)=B\left(x\right)+x$	$B:Bernoullifunction$ ${\mu }_n:electronmobility$ $\begin{array}{r}h:edgelength;Farrell:theclassical\\ fluxj=q\mu U_T[B(-\Delta {\psi }_{KL}/U_T)n_L-B(\Delta {\psi }_{KL}/U_T)n_K]with\\ B(x)=x/(e^x-1),exactforBoltzmann\\ statisticsunderconstantcurrentalong\\ theedge,andthedesiredpropertiesof\\ stability,maximumprinciple,boundary\\ layersandthermodynamicconsistency\\ (divisionbyhandsignconventionwritten\\ herewith\Delta {\psi }_{KL}={\psi }_K-{\psi }_L,seetheQ/A;B\\ identities,backgroundknowledge)\end{array}$ $\begin{array}{r}(computed:B(5)=0.0339,B(-5)=5.034\\ andB(\pm 0.01)=0.995and1.005;at\\ equilibrium{n}_L=e^5{n}_{K}thefluxiszeroto\\ 1\times 10^{-15})\end{array}$
Scharfetter-Gummel generalizations for Fermi-Dirac statistics	$g\left(x\right)=x{\left(F^{-1}\right)}^{\prime }\left(x\right)(\text{generalized SG})$	$\gamma \left(\eta \right)=\frac{F\left(\eta \right)}{e^{\eta }}(\text{inverse-activity scheme})$	$\begin{array}{r}F:statisticsfunction(Fermi-Dirac\\ integralorBlakemoreapproximation)\end{array}$ $\eta :reducedquasi-Fermienergy$ $g:diffusionenhancement$ $\begin{array}{r}\gamma :inverseactivitycoefficient;Farrell:\\ thegeneralSGschemeofKopruckiand\\ Gartner(2013),exactforBlakemore\\ statisticsandsolvedbyafixed-point\\ equationtomachineprecision;diffusion\\ averaging(Bessemoulin-Chatard2012;\\ Kopruckietal.2014)consistentonlyfor\\ aspecificlogarithmicaverage;the\\ inverse-activityscheme(Fuhrmann2015)\\ consistentwithequilibriumforany\\ average;allreducetoclassicalSGfor\\ largenegative\eta ;onapinbenchmarkwith\\ N_A=4.35\times 10^{17}and{N}_D=4.20\times 10^{18}\mathrm{c}{\mathrm{m}}^{-3},about{h}^2\\ convergence,thegeneralSGbestfor\\ currentanddiffusionenhancementgoodfor\\ thefield;codeddfermi\end{array}$ $\begin{array}{r}(computed:second-orderconvergencemeans\\ halvingthemeshsizecutsthecurrent\\ erroraboutfourtimes)\end{array}$

Non-Linear Solvers and Convergence

Name	Formula	Related Formulae	Variables/Constants
Gummel (decoupled) iteration	$F^{\prime }\left({\phi }^t;{\delta }_{\phi }^{t+1}\right)=-F\left({\phi }^t,{\phi }_e^t,{\phi }_h^t\right),{\phi }^{t+1}={\phi }^t+{\delta }_{\phi }^{t+1}$	$F=\nabla \cdot \left(ϵ\nabla \phi \right)+q\left(C+n_i{e}^{({\phi }_h-\phi )/V_T}-n_i{e}^{(\phi -{\phi }_e)/V_T}\right)$	$\begin{array}{r}F:non-linearPoissonresidualwithfrozen\\ quasi-Fermipotentials{\phi }_e,{\phi }_h,F^{\prime }:its\\ Frechetderivative,{\delta }_{\phi }:potentialupdate\end{array}$ $\begin{array}{r}t:iterationindex;ChenandBagci:Gummel\\ decouplesthenon-linearPoissonequation\\ fromthedrift-diffusionequations,which\\ arethensolvedaslinearsystemsfor{n}_e\\ and{n}_{h}withthenewfield,needsless\\ memoryandcomputationperiterationand\\ islesssensitivetotheinitialguess\\ thanNewtonbutconvergesmoreslowly,\\ typicallyinfewerthan50iterations,to\\ atoleranceof1\times 10^{-7};Jeromeand\\ Kerkhoven:Gumme{l}^{\prime }smapastheoldest\\ approach,aGauss-Seidelsuccessive\\ substitution,withconvergenceanalysedin\\ theSlotboom-variableform\end{array}$ $\begin{array}{r}(computed:alinearcontractionfactorof\\ 0.5needs50iterationstoreach1\times 10^{-15})\end{array}$
Newton (coupled) iteration	$J\left({\mathbf{x}}^k\right)\Delta \mathbf{x}=-\mathbf{F}\left({\mathbf{x}}^k\right),{\mathbf{x}}^{k+1}={\mathbf{x}}^k+{\lambda }_k\Delta \mathbf{x}$	$\Vert {\mathbf{x}}^{k+1}-{\mathbf{x}}^{\ast }\Vert \le C{\Vert {\mathbf{x}}^k-{\mathbf{x}}^{\ast }\Vert }^2(\text{near the solution})$	$\mathbf{x}:allnodalunknowns(\psi ,n,p)$ $\begin{array}{r}\mathbf{F}:residualsofthethreediscrete\\ equations\end{array}$ $J:Jacobian$ ${\lambda }_k:dampingfactor(0<{\lambda }_k\le 1)$ $\begin{array}{r}{\mathbf{x}}^{\ast }:solution(forms,background\\ knowledge);JeromeandKerkhoven:adamped\\ NewtonouteriterationwiththeJacobian\\ fromlinearizingthePDEsystem,linear\\ systemssolveddirectlyoriteratively\\ withJacobian-basedpreconditioners,and\\ evidencethatNewtononthediscrete\\ fixed-pointmapistosomeextent\\ superior;ChenandBagci:Newtonasthe\\ fullycoupledalternativetoGummel\end{array}$ $\begin{array}{r}(computed:quadraticconvergencetakesan\\ errorof1\times 10^{-2}to1\times 10^{-4},1\times 10^{-8}and1\times 10^{-16}in\\ threesteps)\end{array}$
Mesh convergence and problem size	${\Vert u-u_h\Vert }_{H^1}\le C{h}^{\theta },0<\theta \le 1$	${ϵ}_I\propto h^2(\text{SG current error, observed})$	$u_h:discretesolutiononmeshsizeh$ $\theta :convergenceorder$ $\begin{array}{r}{ϵ}_I:currenterror;JeromeandKerkhoven:H^1\\ errorboundsoforder{h}^{\theta }for\\ piecewise-linearfiniteelements,withan\\ O(h)andanO(h^2)term,fordimensionsup\\ to3;ChenandBagci:onaMOSFET,SG\\ finitevolumesreach1\times 10^{-2}relativeI-V\\ errorath=1\mathrm{nm}with545342elementsand\\ 1499646unknowns,Galerkinleast-squares\\ finiteelementsneedh=0.3\mathrm{nm}and\\ discontinuousGalerkinh=3\mathrm{nm}(errorsof\\ 3.05\times 10^{-3}and5.28\times 10^{-3}),andDGandSG\\ electrondensitiesdifferby0.78percent\\ onaplasmonicphotoconductiveantenna;\\ Farrell:about{h}^{2}currentconvergence\end{array}$ $\begin{array}{r}(computed:1499646unknownsover545342\\ elementsis2.75unknownsperelement)\end{array}$

Important Questions

Ques: Why can’t simple central differences be used for the current?

Because carrier density changes exponentially between nodes. A central-difference flux treats the density as linear along an edge and averages it in the drift term. At equilibrium with a potential step of 5 thermal voltages across one edge, the density ratio is $e^5=148$, and the central flux comes out at -226 in units of $qD{n}_K/h$, 1.5 times the diffusion term, instead of zero (computed): a large spurious current. The scheme is only stable when the potential step per edge stays below about $2U_T$, 52 mV at room temperature (background knowledge), which would need impractically fine meshes across junctions with built-in potentials near 0.9 $\mathrm{V}$. The Scharfetter-Gummel flux gives exactly zero for the same case (computed), which is why it became the standard. Reference: P. Farrell (with T. Koprucki, J. Fuhrmann, N. Rotundo and H. Doan), ‘Scharfetter-Gummel schemes for non-Boltzmann statistics’, ECMI 2016 talk, Weierstrass Institute Berlin (slides; signs hard to read in the extract) (fetch-tool summary).

Ques: What is the Scharfetter-Gummel flux and why is it exact?

It is the exact solution of the one-dimensional current equation along an edge with constant current and field. Farrell writes the flux as $q\mu U_T[B(-\Delta {\psi }_{KL}/U_T)n_L-B(\Delta {\psi }_{KL}/U_T)n_K]$ with the Bernoulli function $B(x)=x/(e^x-1)$, derived by assuming constant current, solving the resulting ODE for $n(x)$ with the end values $n_K$ and $n_L$. For small potential steps, $B(\pm 0.01)=0.995$ and 1.005 (computed), so the flux reduces to central differencing; for large steps, $B(5)=0.0339$ and $B(-5)=5.034$ (computed), so the flux takes its density almost entirely from the upwind node. It therefore blends diffusion-dominated and drift-dominated transport automatically and is consistent with thermodynamic equilibrium, one of the properties Farrell lists. Reference: P. Farrell (with T. Koprucki, J. Fuhrmann, N. Rotundo and H. Doan), ‘Scharfetter-Gummel schemes for non-Boltzmann statistics’, ECMI 2016 talk, Weierstrass Institute Berlin (slides; signs hard to read in the extract) (fetch-tool summary).

Ques: Is the sign convention of the Scharfetter-Gummel flux in the slides correct?

It is correct if $\Delta {\psi }_{KL}$ means ${\psi }_K-{\psi }_L$. The extracted slide writes $q\mu U_T[B(-\Delta {\psi }_{KL}/U_T)n_L-B(\Delta {\psi }_{KL}/U_T)n_K]$ but warns that its signs are hard to read. At equilibrium, Boltzmann statistics require $n_L=n_K{e}^{({\psi }_L-{\psi }_K)/U_T}$, and with the identity $B(-x)=B(x)+x$, equivalently $B(-x)=e^{x}B(x)$ (background knowledge), the flux vanishes exactly when the argument multiplying $n_L$ is $({\psi }_L-{\psi }_K)/U_T$ (computed: residual below $1\times 10^{-15}$ at a 5 $U_T$ step). With $\Delta {\psi }_{KL}={\psi }_K-{\psi }_L$ the slide’s form satisfies this, while the opposite reading would produce a large equilibrium current. Codes differ in whether they define $\Delta \psi$ as ${\psi }_K-{\psi }_L$ or the reverse, so the equilibrium test is the reliable check (my reading). A useful companion test is a uniform resistor, whose simulated current must match Ohm’s law at any mesh spacing. Reference: P. Farrell (with T. Koprucki, J. Fuhrmann, N. Rotundo and H. Doan), ‘Scharfetter-Gummel schemes for non-Boltzmann statistics’, ECMI 2016 talk, Weierstrass Institute Berlin (slides; signs hard to read in the extract) (fetch-tool summary).

Ques: How is the Scharfetter-Gummel scheme extended to Fermi-Dirac statistics?

In three main ways, compared by Farrell. The generalized SG scheme of Koprucki and Gartner uses $g(x)=x(F^{-1}{)}^{\prime }(x)$ and solves for the flux with a fixed-point equation to machine precision; it is exact for Blakemore statistics. Diffusion averaging replaces g by an edge average built from logarithms of the inverse statistics function and is consistent with equilibrium only for that particular average. Fuhrmann’s inverse-activity scheme uses $\gamma (\eta )=F(\eta )/e^{\eta }$ and is consistent for any average. All three reduce to classical SG in the non-degenerate limit of large negative $\eta$. On a pin benchmark with doping up to $4.20\times 10^{18}$ $\mathrm{c}{\mathrm{m}}^{-3}$, all converged at about $h^2$; the generalized SG gave the best current and diffusion enhancement the best field. Degenerate regions such as contacts, lasers and LEDs need these schemes. Reference: P. Farrell (with T. Koprucki, J. Fuhrmann, N. Rotundo and H. Doan), ‘Scharfetter-Gummel schemes for non-Boltzmann statistics’, ECMI 2016 talk, Weierstrass Institute Berlin (slides; signs hard to read in the extract) (fetch-tool summary).

Ques: When is Gummel iteration preferred over Newton’s method?

When memory or a poor initial guess is the limitation, and coupling is weak. Chen and Bagci state that Gummel needs less memory and computation per iteration and is less sensitive to the initial guess than Newton, but converges more slowly, typically in fewer than 50 iterations to a tolerance of $1\times 10^{-7}$ for their devices. Jerome and Kerkhoven describe it as the oldest approach, a Gauss-Seidel successive substitution, whose convergence has been analysed in Slotboom variables. Gummel’s linear convergence slows when the equations are strongly coupled, as at high injection or high current, where Newton is preferred (background knowledge). A linear contraction factor of 0.5 needs 50 iterations to reach $1\times 10^{-15}$ (computed). Simulators often start with Gummel iterations and switch to Newton once close to the solution (background knowledge). Reference: L. Chen and H. Bagci, ‘Steady-state Simulation of Semiconductor Devices using Discontinuous Galerkin Methods’, arXiv:1910.05977 (fetch-tool summary); J. W. Jerome and T. Kerkhoven, ‘A Finite Element Approximation Theory for the Drift Diffusion Semiconductor Model’ (received 1989, accepted 1990; preprint, journal not stated) (fetch-tool summary).

Ques: How fast does Newton’s method converge, and why is it damped?

Quadratically near the solution, so it is damped far from it. Newton solves $J\Delta \mathbf{x}=-\mathbf{F}$ for all potentials and densities together, and near the solution the error roughly squares each step: from $1\times 10^{-2}$ to $1\times 10^{-4}$, $1\times 10^{-8}$ and $1\times 10^{-16}$ in three steps (computed). Far from the solution the exponential carrier terms make full steps overshoot, so a damping factor ${\lambda }_k\le 1$ shortens the update (background knowledge). Jerome and Kerkhoven describe such a damped Newton outer iteration with the Jacobian formed by linearizing the PDE system and linear systems solved directly or iteratively with Jacobian-based preconditioners. Each Newton step costs more than a Gummel step because the coupled Jacobian is three times larger in each dimension, as Chen and Bagci note in comparing memory needs. Reference: J. W. Jerome and T. Kerkhoven, ‘A Finite Element Approximation Theory for the Drift Diffusion Semiconductor Model’ (received 1989, accepted 1990; preprint, journal not stated) (fetch-tool summary); L. Chen and H. Bagci, ‘Steady-state Simulation of Semiconductor Devices using Discontinuous Galerkin Methods’, arXiv:1910.05977 (fetch-tool summary).

Ques: How fine must a TCAD mesh be?

Fine enough that the extracted currents stop changing, which depends on the scheme. Chen and Bagci report that for a MOSFET, SG finite volumes reach $1\times 10^{-2}$ relative I-V error at a 1 $\mathrm{nm}$ mesh, with 545342 elements and about 1.5 million unknowns, 2.75 per element (computed), while Galerkin least-squares finite elements need 0.3 $\mathrm{nm}$ and discontinuous Galerkin only 3 $\mathrm{nm}$. Farrell observes roughly $h^2$ current convergence for SG schemes, so halving the mesh cuts the error about four times (computed), and Jerome and Kerkhoven prove $h^{\theta }$ error bounds with $\theta$ up to 1 for linear elements. In practice the mesh must resolve the Debye length near junctions, a few nanometres at high doping (Drift-Diffusion Device Equations page), and is refined where doping, potential or current change fastest (background knowledge). Reference: L. Chen and H. Bagci, ‘Steady-state Simulation of Semiconductor Devices using Discontinuous Galerkin Methods’, arXiv:1910.05977 (fetch-tool summary); P. Farrell (with T. Koprucki, J. Fuhrmann, N. Rotundo and H. Doan), ‘Scharfetter-Gummel schemes for non-Boltzmann statistics’, ECMI 2016 talk, Weierstrass Institute Berlin (slides; signs hard to read in the extract) (fetch-tool summary); J. W. Jerome and T. Kerkhoven, ‘A Finite Element Approximation Theory for the Drift Diffusion Semiconductor Model’ (received 1989, accepted 1990; preprint, journal not stated) (fetch-tool summary); Drift-Diffusion Device Equations page of this vault.

References

1. P. Farrell, T. Koprucki, J. Fuhrmann, N. Rotundo and H. Doan, “Scharfetter-Gummel schemes for non-Boltzmann statistics,” ECMI 2016 talk, Weierstrass Institute for Applied Analysis and Stochastics, Berlin, 2016. Available: https://www.wias-berlin.de/people/peschka/files/ecmi2016_farrell.pdf
2. L. Chen and H. Bagci, “Steady-state Simulation of Semiconductor Devices using Discontinuous Galerkin Methods,” arXiv:1910.05977. Available: https://arxiv.org/abs/1910.05977
3. J. W. Jerome and T. Kerkhoven, “A Finite Element Approximation Theory for the Drift Diffusion Semiconductor Model,” preprint (received 1989, accepted 1990). Available: https://sites.math.northwestern.edu/~jwj/preprints/Kerkhoven/Kerkhoven.pdf

Related pages in this vault: Drift-Diffusion Device Equations, Thermal and Electro-Thermal Simulation, Band-Gap Narrowing and Quantum Corrections, Hydrodynamic and Energy-Balance Transport, Generation-Recombination Models, Basic Semiconductor Equations.