**The drift-diffusion current is stable because of the Bernoulli function, not because of the physics.** The Scharfetter–Gummel discretisation writes the inter-node current with the Bernoulli weight $B(x)=x/(e^x-1)$, and a central-difference approximation replaces that weight by its linear truncation $B(x)\approx 1-x/2$, which goes negative once a single cell drops more than $2\,kT/q = 51.7$ mV — a cell Péclet number of 2. Solve the textbook boundary layer both ways and the consequence is stark: central differencing drives the interior carrier density to -0.70 — a negative concentration, which is physically impossible — while Scharfetter–Gummel is exact at every mesh spacing. The 1968 scheme is not a numerical nicety; it is the reason a device simulator returns positive densities at all, and it is invisible in every glossy Id–Vg plot the tool produces. **A TCAD solve is a Newton iteration, so the solver is as much the model as the physics is.** The coupled Poisson–continuity system is nonlinear because the carrier densities depend exponentially on potential, and it is solved by Newton–Raphson, which squares its residual each step near the solution. Started from a sane guess the surface-potential residual falls ${3.4\times 10^{-1}}$ to ${2.7\times 10^{-2}}$ to ${2.3\times 10^{-3}}$ to ${1.9\times 10^{-5}}$ to ${1.3\times 10^{-9}}$ to ${8.3\times 10^{-17}}$ — quadratic convergence in 5 Newton steps to machine precision. Started far from the solution the same exponential overflows and the step must be damped or taken in Gummel's decoupled order instead, or the solve simply diverges and returns nothing. Whether an answer comes out, and which answer, is a property of the initial guess and the damping, which is why two engineers running the same deck can disagree — a numerical fact, not a physical one. **The mesh is a modeling decision that quietly changes the answer.** In strong inversion the electron sheet sits within a few extended Debye lengths of the surface, and $L_D=\sqrt{\varepsilon\, kT/(q^2 n)}$ is only 1.293 nm at an inversion density of $10^{19}$ cm$^{-3}$, against 12.929 nm in the lightly doped bulk. Over a 30 nm body that means the first grid cell must sit within about 0.5 nm of the interface or the inversion charge — and therefore the drive current — is simply wrong, and only once the layer is resolved does the error fall as the square of the spacing. A coarse mesh does not merely add noise; below two grid points per Debye length it changes the threshold voltage a designer reads off the curve. The grid is a knob with no physics in it, and it is set by the engineer, not the transistor. **Every mobility model is a curve fit, and swapping one moves the drive current with no new physics at all.** Low-field mobility, the Caughey–Thomas field dependence $\mu(E)=\mu_0/[1+(\mu_0 E/v_\text{sat})^\beta]^{1/\beta}$ with a saturation velocity of $v_\text{sat}=10^7$ cm/s, and the Lombardi surface-roughness model are three fits stacked on top of each other. At a 20 nm gate the lateral field is high enough that a constant-mobility model overpredicts the carrier velocity by 98% relative to the velocity-saturated form — an enormous swing in $I_\text{on}$ produced entirely by which empirical curve the engineer selected from a menu. The transistor did not change; the model did. **Heavy-doping band-gap narrowing is a fitted correction that multiplies every injection and generation current.** The Slotboom–de Graaff form shrinks the gap by 60 meV at $10^{19}$ cm$^{-3}$ and 92 meV at $10^{20}$ cm$^{-3}$, and because the effective intrinsic density enters squared, $n_i^2$ is enhanced by $e^{\Delta E_g/kT}$ — a factor of 10.4$\times$ and 35.4$\times$ respectively. A parameter fit to one process node therefore rescales bipolar gain, junction leakage and source-drain injection by more than an order of magnitude, which is why an uncalibrated deck applied to a new doping recipe can be confidently, quietly wrong. **Reliability and leakage predictions are only as good as the lifetimes and cross-sections you fit.** Off-state generation current through Shockley–Read–Hall traps scales as $1/\tau$, so a factor-of-2 error in the fitted carrier lifetime is a factor-of-2 error in predicted $I_\text{off}$ — linear and direct. NBTI and hot-carrier aging are worse: their trap-generation kinetics are empirical power laws whose exponents are fit to stress data, so a ten-year extrapolation is an extrapolation of a fit, not a derivation from physics. The standard model set below is a stack of such fits, and each row is a place where a number was chosen to match silicon. | Model | Role in the solve | Fitted handle | |-------|-------------------|----------------| | Drift-diffusion | Carrier transport, the default | mobilities $\mu_n,\mu_p$ | | Caughey–Thomas | Velocity saturation at high field | $v_\text{sat}$, $\beta$ | | Lombardi surface | Mobility degradation at the interface | roughness and Coulomb terms | | Slotboom BGN | Band-gap narrowing at heavy doping | $\Delta E_g$ prefactor | | SRH + Auger | Recombination and off-state leakage | lifetimes $\tau_n,\tau_p$ | | van Overstraeten | Impact ionisation, breakdown | ionisation coefficients | **Three dimensions is not three-halves the work; it is a different cost class.** A sparse Newton solve on $N$ unknowns costs roughly $N^{1.5}$ in memory and up to $N^2$ in time for a direct factorisation, so moving from a 2D cross-section of order $10^4$ nodes to a full 3D FinFET or gate-all-around structure of order $10^6$ nodes is 100$\times$ the unknowns but a 1,000x jump in factorisation memory and far more in solve time. That single scaling law is why production flows still lean on 2D splits, symmetry, and overnight runs on HPC clusters, and why a 3D reliability sweep is a capital-planning decision rather than a coffee break. **Calibration is the whole game: a deck predicts nothing until it matches measured silicon.** A modern device deck exposes on the order of 40 adjustable parameters across the transport, mobility, band and recombination models, and the flow is always the same — run split-lot wafers, measure threshold voltage, off-current, on-current and subthreshold slope, then tune parameters until the simulated curves sit within about ±5% of the data. Only inside that calibrated envelope is the ±10–15% predictive accuracy earned; push the geometry or the doping outside it and the deck reverts to a plausible-looking extrapolation of a fit. Done well it is what lets a foundry compress a development cycle by 30–50% against pure wafer experiments and cut the number of costly split lots by a comparable margin; foundry and IDM PDK teams at TSMC, Intel and imec keep this calibration alive precisely because it is the difference between a virtual fab and a physics-flavoured guess. ```flowchart Process recipe -> [Process TCAD] -> structure (doping, geometry) | v [Device TCAD] -> Id-Vg, Id-Vd, Vt, SS, leakage | [Compact-model extraction] -> SPICE parameters | v [Circuit simulation] -> ring oscillator, SRAM timing ^ | |__________ calibrate to split-lot silicon (+/-5%) <_________| ``` Read TCAD through a *calibration* lens rather than a *first-principles* lens, and every hard problem on this page becomes the same problem: the Bernoulli-stabilised scheme, the Newton damping, the Debye-resolved mesh, the mobility and band-gap and lifetime fits, and the cost of the third dimension are all knobs that the engineer sets and silicon adjudicates, not truths the physics hands over for free. A TCAD deck is a hypothesis about a transistor that has been argued into agreement with measured wafers; its power is real, but it is the power of a well-calibrated instrument, and the moment it is used outside the window it was fit to, it predicts with exactly the confidence of an extrapolated fit and no more.