Drift-Diffusion Device Equations

1. Theory

• Definition: The drift-diffusion model is the standard set of partial differential equations solved by device TCAD: Poisson’s equation for the potential and two continuity equations for electrons and holes, with drift-diffusion current relations (Bortolossi; del Alamo).
• Unknowns: Three fields, the potential and the two carrier densities or quasi-Fermi potentials, are solved on a mesh of the device with doping, materials and contacts (Bortolossi; DEVSIM).
• Assumptions: Carriers stay close to thermal equilibrium with the lattice, so their distributions are described by quasi-Fermi levels and the Einstein relation holds (del Alamo; background knowledge for the Einstein link).
• Closure: Physical models for mobility, recombination and band structure complete the system; they are covered on the following pages (background knowledge).

2. Working Principle

• Electrostatics: Charge from carriers and ionized dopants sets the potential through Poisson’s equation (Bortolossi; del Alamo).
• Transport: Carriers drift in the field and diffuse down density gradients; together these equal flow down the gradient of the quasi-Fermi potential (Bortolossi).
• Conservation: The divergence of each current balances generation and recombination and any change of stored carriers (del Alamo; Bortolossi).
• Boundaries: Ohmic contacts fix neutral, equilibrium densities and the applied voltage; insulating surfaces carry no current or flux (Bortolossi).
• Discretization: Finite-volume or finite-element methods turn the equations into coupled non-linear algebraic equations solved by Newton iteration (DEVSIM; Bortolossi; background knowledge for Newton).

3. Applications

• Device design: MOSFETs, diodes, BJTs, solar cells and power devices before fabrication (Bortolossi; background knowledge).
• Compact-model support: Generating I-V and C-V data and variability for model extraction (background knowledge; Process Corners and Statistical Models page).
• Physics insight: Mapping fields, carrier densities and quasi-Fermi levels inside devices (del Alamo; Bortolossi).
• Open-source tools: DEVSIM for scripted DC, AC, noise and transient analyses in 1-D, 2-D and 3-D (DEVSIM).

Device Equations

Name	Formula	Related Formulae	Variables/Constants
Poisson’s equation	$-\nabla \cdot \left(ϵ\nabla \psi \right)=q\left(p-n+N_D^{+}-N_A^{-}\right)$	$\vec{E}=-\nabla \psi ,\nabla \cdot \vec{D}=\rho$	$ϵ:permittivity$ $\psi :electrostaticpotential$ $q:elementarycharge$ $n,p:electronandholedensities$ $N_D^{+},N_A^{-}:ionizeddonorandacceptordensities$ $\vec{E}:field$ $\begin{array}{r}\vec{D}:displacement;Bortolossi:the\\ stationarysystemwithnetdopingD=N_D^{+}-\\ N_A^{-},writtenas\nabla \cdot (-\nabla \phi )-q(p-n)=qDin\\ theextract(seetheQ/A);delAlamo:\\ Gaus{s}^{\prime }slaw\nabla \cdot E=(q/ϵ)(p-n+N_D^{+}-N_A^{-});\\ DEVSIM:theintegralformϵ\oint \nabla \psi \cdot d\vec{s}\\ assembledasanedgefluxina\\ finite-volumescheme;1-Dformonthe\\ BasicSemiconductorEquationspage\end{array}$ $(nocomputedvalue)$
Continuity equations	$\frac{\partial n}{\partial t}=\frac{1}{q}\nabla \cdot {\vec{J}}_n+G-R,\frac{\partial p}{\partial t}=-\frac{1}{q}\nabla \cdot {\vec{J}}_p+G-R$	$\nabla \cdot {\vec{J}}_n=qR,\nabla \cdot {\vec{J}}_p=-qR(\text{stationary, }G=0)$	${\vec{J}}_n,{\vec{J}}_p:electronandholecurrentdensities$ $G:generationrate$ $R:recombinationrate$ $\begin{array}{r}t:time;Bortolossi:stationarycontinuity\\ \nabla \cdot (q{\mu }_{n}n\nabla \phi -q{D}_n\nabla n)=-qRandthehole\\ analogueinsilicon,withn=p=0and{J}_n\\ =J_p=0inoxide;delAlamo:carrier\\ dynamicsdn/dt=dp/dt=G-Randtotal\\ currentJ=J_e+J_h\end{array}$ $(nocomputedvalue)$
Drift-diffusion current densities	${\vec{J}}_n=qn{\mu }_n\vec{E}+q{D}_n\nabla n,{\vec{J}}_p=qp{\mu }_p\vec{E}-q{D}_p\nabla p$	$D_{n,p}=V_{th}{\mu }_{n,p},V_{th}=\frac{k_{B}T}{q}$	${\mu }_n,{\mu }_p:mobilities$ $D_n,D_p:diffusioncoefficients$ $V_{th}:thermalvoltage$ $k_B:Boltzmannconstant$ $\begin{array}{r}T:temperature;Bortolossi:theseforms\\ andtheEinsteinrelation,withlow-field\\ mobilitiesof1417and470.5\mathrm{c}{\mathrm{m}}^2/(\mathrm{V}\mathrm{s})in\\ thetestcase;delAlamo:J_{e}and{J}_h\\ (electronsignasextracted,seetheQ/A)\end{array}$ $(nocomputedvalue)$

Quasi-Fermi Formulation

Name	Formula	Related Formulae	Variables/Constants
Carrier densities and quasi-Fermi potentials	$n=n_i\exp \left(\frac{\psi -{\phi }_n}{V_{th}}\right),p=n_i\exp \left(\frac{{\phi }_p-\psi }{V_{th}}\right)$	$np=n_i^2\exp \left(\frac{{\phi }_p-{\phi }_n}{V_{th}}\right),n=N_c{F}_{1/2}\left(\frac{E_{Fn}-E_c}{k_{B}T}\right)(\text{Fermi-Dirac})$	$n_i:intrinsicdensity$ ${\phi }_n,{\phi }_p:electronandholequasi-Fermipotentials$ $\begin{array}{r}{N}_c:conduction-bandeffectivedensityof\\ states\end{array}$ $F_{1/2}:Fermi-Diracintegral$ $E_{Fn}:electronquasi-Fermilevel$ $\begin{array}{r}{E}_c:conduction-bandedge;Bortolossi:the\\ Boltzmannrelationswith\psi themidgap\\ potentialand{n}_i=\sqrt{N_c{N}_v}exp(-E_g/(2k_{B}T)),E_g=\\ 1.11\mathrm{eV};delAlamo:theFermi-Diracand\\ Maxwell-Boltzmannforms,np=n_i^{2}exp((E_{fe}-\\ E_{fh})/k_{B}T),andquasi-Fermilevelsvalidwhen\\ distributionsstaynearthermal\\ equilibrium(energy-levelform;potentials\\ hereintheconvention{E}_F=-q\phi )\end{array}$ $(nocomputedvalue)$
Quasi-Fermi form of the currents	${\vec{J}}_n=-qn{\mu }_n\nabla {\phi }_n,{\vec{J}}_p=-qp{\mu }_p\nabla {\phi }_p$	$u_n=n_i{e}^{-{\phi }_n/V_{th}},{\vec{J}}_n=q{D}_n{e}^{\psi /V_{th}}\nabla u_n,u_p=n_i{e}^{{\phi }_p/V_{th}},{\vec{J}}_p=-q{D}_p{e}^{-\psi /V_{th}}\nabla u_p$	$\begin{array}{r}{u}_n,u_p:Slotboomvariables;Bortolossi:\\ bothforms(equations1.36to1.41);a\\ flatquasi-Fermipotentialmeanszero\\ currentofthatcarrier(background\\ knowledge)\end{array}$ $(nocomputedvalue)$

Boundary Conditions and Scales

Name	Formula	Related Formulae	Variables/Constants
Ohmic contact boundary conditions	$n_D=\frac{N+\sqrt{N^2+4n_i^2}}{2},p_D=\frac{-N+\sqrt{N^2+4n_i^2}}{2},N=N_D^{+}-N_A^{-}$	${\phi }_n={\phi }_p=V_{ext},{\psi }_D=V_{ext}+V_{th}\ln \frac{n_D}{n_i}$	$n_D,p_D:carrierdensitiesatthecontact$ $N:netdoping$ $V_{ext}:appliedcontactvoltage$ $\begin{array}{r}{\psi }_D:contactpotential;Bortolossi:\\ Dirichletconditionsfromcharge\\ neutralityandequilibrium{p}_D{n}_D=n_i^2,these\\ solutions,quasi-Fermipotentialsequalto\\ theappliedvoltage,andthecontact\\ potential\end{array}$ $\begin{array}{r}(computed:with{n}_i=9.65\times 10^9\mathrm{c}{\mathrm{m}}^{-3}(assumed),\\ a1\times 10^{18}\mathrm{c}{\mathrm{m}}^{-3}n-contactsitsat+0.477\mathrm{V}and\\ a1\times 10^{17}\mathrm{c}{\mathrm{m}}^{-3}p-contactat-0.418\mathrm{V}relative\\ totheappliedvoltage,sothejunction\\ betweenthemhasabuilt-inpotentialof\\ 0.895\mathrm{V})\end{array}$
Insulating and artificial boundaries	$\vec{D}\cdot \hat{n}=0,{\vec{J}}_n\cdot \hat{n}=0,{\vec{J}}_p\cdot \hat{n}=0$	${ϵ}_1\frac{\partial \psi }{\partial n}{|}_1-{ϵ}_2\frac{\partial \psi }{\partial n}{|}_2=Q_{int}(\text{interface})$	$\hat{n}:outwardnormal$ $\begin{array}{r}{Q}_{int}:interfacechargeperarea(interface\\ condition,backgroundknowledge);\\ Bortolossi:homogeneousNeumannconditions\\ oninsulatingandartificialboundaries\\ andzerocarriersandcurrentsinoxide\\ regions\end{array}$ $(nocomputedvalue)$
Characteristic length and time scales	$L_D=\sqrt{\frac{ϵV_{th}}{qN}},{\tau }_d=\frac{ϵ}{\sigma }=\frac{ϵ}{qn{\mu }_n}$	$\psi \to \frac{\psi }{V_{th}},n,p\to \frac{n,p}{n_i}(\text{scaling})$	$L_D:extrinsicDebyelength$ $N:doping$ ${\tau }_d:dielectricrelaxationtime$ $\begin{array}{r}\sigma :conductivity(definitionsandscaling,\\ backgroundknowledge);ameshmustresolve\\ L_{D}nearjunctions,anddensitiesspanmany\\ decades,whichiswhysimulatorsuse\\ scaledandlogarithmicorSlotboom\\ variables(backgroundknowledge);\\ Bortolossi:a0.05by0.05by0.1\mu \mathrm{m}\\ junctionmeshedwith4933verticesand\\ 24576elements,biasedat0.3and1.0\mathrm{V}\end{array}$ $\begin{array}{r}(computed:L_D=12.9\mathrm{nm}at1\times 10^{17}and4.1\mathrm{nm}\\ at1\times 10^{18}\mathrm{c}{\mathrm{m}}^{-3},and41.6\mu \mathrm{m}forintrinsic\\ silicon;{\tau }_d=81\mathrm{fs}at1\times 10^{17}\mathrm{c}{\mathrm{m}}^{-3}with800\\ \mathrm{c}{\mathrm{m}}^2/(\mathrm{V}\mathrm{s});thecarrierdensityina1\times 10^{18}\\ n-regionspansfrom1\times 10^{18}downto93\mathrm{c}{\mathrm{m}}^{-3}\\ forholes,16decades)\end{array}$

Important Questions

Ques: What equations does a drift-diffusion device simulator solve?

Poisson’s equation and the electron and hole continuity equations, with drift-diffusion currents. Bortolossi writes the stationary system for potential, electrons and holes on a 3-D domain, with currents $qn{\mu }_n\vec{E}+q{D}_n\nabla n$ and $qp{\mu }_p\vec{E}-q{D}_p\nabla p$ and the Einstein relation $D=V_{th}\mu$. Del Alamo’s five equations of semiconductor physics are the same set: Gauss’s law, two current equations, total current, and carrier dynamics $dn/dt=dp/dt=G-R$. DEVSIM assembles such equations by the finite-volume method, integrating fluxes over control-volume surfaces. With three unknowns per mesh node, Bortolossi’s test junction of 4933 vertices already has about 14,800 unknowns (computed). Everything else in TCAD, mobility, recombination and band structure, enters as coefficients of these equations. Reference: A. Bortolossi, ‘3D Finite Element Drift-Diffusion Simulation of Semiconductor Devices’, Laurea Magistrale thesis in Mathematical Engineering, Politecnico di Milano (academic year 2013-2014) (fetch-tool summary); J. A. del Alamo, ‘6.720J/3.43J Integrated Microelectronic Devices’, Lecture 8: Carrier Drift and Diffusion (cont.), Carrier Flow, MIT OpenCourseWare (21 February 2007) (equations partly garbled in the extract) (fetch-tool summary); DEVSIM LLC, ‘Introduction’, DEVSIM documentation (J. E. Sanchez; Apache License 2.0) (fetch-tool summary).

Ques: Why do simulators use quasi-Fermi potentials?

Because they give the current in its simplest form and handle the huge range of carrier densities. Bortolossi shows that the drift-diffusion current equals $-qn{\mu }_n\nabla {\phi }_n$: drift and diffusion combine into flow down the gradient of the electron quasi-Fermi potential, so a flat ${\phi }_n$ means zero electron current (background knowledge for the corollary). In the Slotboom variables $u_n=n_i{e}^{-{\phi }_n/V_{th}}$ the current becomes $q{D}_n{e}^{\psi /V_{th}}\nabla u_n$, which leads to a self-adjoint, numerically convenient form. Carrier densities in one device span about 16 decades, from $1\times 10^{18}$ $\mathrm{c}{\mathrm{m}}^{-3}$ majority electrons to 93 $\mathrm{c}{\mathrm{m}}^{-3}$ minority holes in the same region (computed), whereas quasi-Fermi potentials vary by at most the applied voltage. Del Alamo notes that quasi-Fermi levels are valid as long as distributions stay near thermal equilibrium. Reference: A. Bortolossi, ‘3D Finite Element Drift-Diffusion Simulation of Semiconductor Devices’, Laurea Magistrale thesis in Mathematical Engineering, Politecnico di Milano (academic year 2013-2014) (fetch-tool summary); J. A. del Alamo, ‘6.720J/3.43J Integrated Microelectronic Devices’, Lecture 8: Carrier Drift and Diffusion (cont.), Carrier Flow, MIT OpenCourseWare (21 February 2007) (equations partly garbled in the extract) (fetch-tool summary).

Ques: What conditions does an ohmic contact impose?

Charge neutrality, thermal equilibrium and the applied voltage. Bortolossi sets $p_D{n}_D=n_i^2$ and $p_D-n_D+N=0$ at the contact, giving $n_D=(N+\sqrt{N^2+4n_i^2})/2$, sets both quasi-Fermi potentials equal to the applied voltage, and fixes the potential at $V_{ext}+V_{th}\ln (n_D/n_i)$. With an assumed $n_i$ of $9.65\times 10^9$ $\mathrm{c}{\mathrm{m}}^{-3}$, a $1\times 10^{18}$ $\mathrm{c}{\mathrm{m}}^{-3}$ n-type contact sits 0.477 $\mathrm{V}$ above its applied voltage and a $1\times 10^{17}$ $\mathrm{c}{\mathrm{m}}^{-3}$ p-type contact 0.418 $\mathrm{V}$ below it, so with both grounded the device holds a built-in potential of 0.895 $\mathrm{V}$ internally (computed). An ideal ohmic contact thus has infinite recombination velocity; Schottky contacts replace these Dirichlet conditions with thermionic-emission currents (background knowledge; Metal Semiconductor Contacts page). Reference: A. Bortolossi, ‘3D Finite Element Drift-Diffusion Simulation of Semiconductor Devices’, Laurea Magistrale thesis in Mathematical Engineering, Politecnico di Milano (academic year 2013-2014) (fetch-tool summary).

Ques: Does Poisson’s equation in the thesis omit the permittivity?

As extracted, yes, but only apparently. The extract of Bortolossi’s stationary system reads $\nabla \cdot (-\nabla \phi )-q(p-n)=qD$, which lacks $ϵ$ and has the charge terms on both sides, while del Alamo writes Gauss’s law as $\nabla \cdot E=(q/ϵ)(p-n+N_D^{+}-N_A^{-})$. Rearranged, the thesis form is $-{\nabla }^2\phi =q(p-n+D)$, the same equation with $ϵ$ set to 1 (computed). Either the permittivity was lost in extraction, as happens with inline symbols, or the system is written in scaled units; the thesis also models oxide, where $ϵ$ changes, so it must appear in the full text (my reading). The standard form $-\nabla \cdot (ϵ\nabla \psi )=q(p-n+N_D^{+}-N_A^{-})$ is certainly correct, and $ϵ$ must stay inside the divergence for multi-material devices. Reference: A. Bortolossi, ‘3D Finite Element Drift-Diffusion Simulation of Semiconductor Devices’, Laurea Magistrale thesis in Mathematical Engineering, Politecnico di Milano (academic year 2013-2014) (fetch-tool summary); J. A. del Alamo, ‘6.720J/3.43J Integrated Microelectronic Devices’, Lecture 8: Carrier Drift and Diffusion (cont.), Carrier Flow, MIT OpenCourseWare (21 February 2007) (equations partly garbled in the extract) (fetch-tool summary).

Ques: Are the signs of the electron current consistent between the sources?

Not in the extracts. Bortolossi gives ${\vec{J}}_n=qn{\mu }_n\vec{E}+q{D}_n\nabla n$, while the extract of del Alamo’s slide reads $J_e=-qn{\mu }_{e}E+q{D}_e\nabla n$. The conventional form has a plus sign on the drift term: electrons drift against the field, but their negative charge makes the current flow along it, and diffusion current flows up the electron gradient because electrons flow down it (background knowledge). The same lecture’s hole current, $qp{\mu }_{h}E-q{D}_h\nabla p$, agrees with Bortolossi, and the extract warns that its equations were reconstructed from garbled text, so the minus sign is almost certainly an extraction artefact (my reading). A sign error here would make electrons flow up the potential and is easy to catch with a resistor test: the simulated current must equal $V/R$ with $R=L/(qn{\mu }_{n}A)$. Reference: A. Bortolossi, ‘3D Finite Element Drift-Diffusion Simulation of Semiconductor Devices’, Laurea Magistrale thesis in Mathematical Engineering, Politecnico di Milano (academic year 2013-2014) (fetch-tool summary); J. A. del Alamo, ‘6.720J/3.43J Integrated Microelectronic Devices’, Lecture 8: Carrier Drift and Diffusion (cont.), Carrier Flow, MIT OpenCourseWare (21 February 2007) (equations partly garbled in the extract) (fetch-tool summary).

Ques: What length and time scales must a device simulation resolve?

The Debye length in space and the dielectric relaxation time in time. The extrinsic Debye length $\sqrt{ϵV_{th}/(qN)}$ is 12.9 $\mathrm{nm}$ at $1\times 10^{17}$ $\mathrm{c}{\mathrm{m}}^{-3}$ and 4.1 $\mathrm{nm}$ at $1\times 10^{18}$ $\mathrm{c}{\mathrm{m}}^{-3}$, and 41.6 $\mu \mathrm{m}$ for intrinsic silicon (computed), so meshes must be finest at junctions and heavily doped regions, where potential and density change over a few Debye lengths (background knowledge). The dielectric relaxation time $ϵ/\sigma$ is 81 $\mathrm{fs}$ at $1\times 10^{17}$ $\mathrm{c}{\mathrm{m}}^{-3}$ with 800 $\mathrm{c}{\mathrm{m}}^2/(\mathrm{V}\mathrm{s})$ (computed), far shorter than circuit time steps, which is why majority carriers can be treated as quasi-neutral and why fully coupled implicit solvers are needed for transients (background knowledge). Bortolossi’s 0.05 by 0.05 by 0.1 $\mu \mathrm{m}$ junction used 24576 elements. Reference: A. Bortolossi, ‘3D Finite Element Drift-Diffusion Simulation of Semiconductor Devices’, Laurea Magistrale thesis in Mathematical Engineering, Politecnico di Milano (academic year 2013-2014) (fetch-tool summary).

Ques: When does the drift-diffusion model break down?

When carriers are far from equilibrium with the lattice. Del Alamo ties quasi-Fermi levels to the assumption that carrier energy distributions never depart too far from thermal equilibrium on practical time scales. In very short channels, high fields heat carriers, velocity overshoot and ballistic transport appear, and the local mobility and diffusion coefficients no longer describe transport; energy-balance, hydrodynamic or Monte Carlo methods are then used (background knowledge; Hydrodynamic and Energy-Balance Transport page, planned). In nanometre inversion layers and thin films, quantum confinement shifts the charge away from the interface, which classical Boltzmann statistics miss; density-gradient or Schrodinger-Poisson corrections are added (background knowledge). At high doping, Fermi-Dirac statistics replace the Boltzmann relations, as del Alamo’s $F_{1/2}$ form shows. Reference: J. A. del Alamo, ‘6.720J/3.43J Integrated Microelectronic Devices’, Lecture 8: Carrier Drift and Diffusion (cont.), Carrier Flow, MIT OpenCourseWare (21 February 2007) (equations partly garbled in the extract) (fetch-tool summary).

References

1. A. Bortolossi, “3D Finite Element Drift-Diffusion Simulation of Semiconductor Devices,” Laurea Magistrale thesis in Mathematical Engineering, Politecnico di Milano, 2014. Available: https://www.politesi.polimi.it/retrieve/a81cb05a-eff6-616b-e053-1605fe0a889a/2014_7_Bortolossi.pdf
2. J. A. del Alamo, “Carrier Drift and Diffusion (cont.), Carrier Flow,” Lecture 8, 6.720J/3.43J Integrated Microelectronic Devices, MIT OpenCourseWare, Spring 2007. Available: https://ocw.tau.edu.ng/courses/electrical-engineering-and-computer-science/6-720j-integrated-microelectronic-devices-spring-2007/lecture-notes/lecture8.pdf
3. DEVSIM LLC, “Introduction,” DEVSIM documentation (J. E. Sanchez). Available: https://devsim.org/introduction.html

Related pages in this vault: Hydrodynamic and Energy-Balance Transport, Discretization and Numerical Solution, Band-Gap Narrowing and Quantum Corrections, Thermal and Electro-Thermal Simulation, Basic Semiconductor Equations, Generation-Recombination Models, Metal Semiconductor Contacts, Electrostatics and Capacitance, Magnetism, Flux and Inductance.