Physics-based semiconductor modeling converts conservation laws, material relations, geometry, and operating conditions into a boundary-value or initial-boundary-value problem whose solution is useful only when the equations, numerical approximation, and comparison with experiment are each tested separately.

<svg viewBox="0 0 760 470" xmlns="http://www.w3.org/2000/svg" font-family="-apple-system,BlinkMacSystemFont,Segoe UI,Roboto,sans-serif">
  <rect x="0" y="0" width="760" height="470" rx="18" fill="#0d1117"/>
  <text x="380" y="28" fill="#e6edf3" font-size="21" font-weight="700" text-anchor="middle">Physics Based Modeling And Differential Equations Technical Microarchitecture</text>
  <text x="380" y="48" fill="#8b98a5" font-size="12" text-anchor="middle">Detailed Domain Pipeline, Architectural Blocks &amp; Engineering Performance Optimization (ID 10702)</text>
  
  
    <!-- FAB CROSS SECTION (2 Main Panels) -->
    <g transform="translate(25, 75)">
      <rect width="345" height="325" fill="#161b22" stroke="#30363d" stroke-width="1.5" rx="8"/>
      <text x="172.5" y="25" fill="#34d399" font-size="13" font-weight="700" text-anchor="middle">1. Physical Layer Cross-Section</text>
      <rect x="20" y="240" width="305" height="60" fill="#0d1117" stroke="#30363d" rx="4"/>
      <text x="172.5" y="275" fill="#8b98a5" font-size="11" font-weight="600" text-anchor="middle">Silicon Substrate / Base Crystal Wafers</text>
      <rect x="20" y="170" width="305" height="60" fill="#1f2937" stroke="#047857" rx="4"/>
      <text x="172.5" y="205" fill="#6ee7b7" font-size="11" font-weight="600" text-anchor="middle">Dielectric Oxide &amp; Isolation Barriers</text>
      <rect x="20" y="100" width="305" height="60" fill="#0f172a" stroke="#34d399" stroke-width="1.5" rx="4"/>
      <text x="172.5" y="135" fill="#e6edf3" font-size="12" font-weight="700" text-anchor="middle">Active Junctions &amp; Nanometer Channel</text>
      <rect x="40" y="45" width="70" height="45" fill="#10b981" rx="3"/>
      <text x="75" y="72" fill="#ffffff" font-size="10" font-weight="700" text-anchor="middle">Source</text>
      <rect x="137.5" y="45" width="70" height="45" fill="#fbbf24" rx="3"/>
      <text x="172.5" y="72" fill="#0d1117" font-size="10" font-weight="700" text-anchor="middle">Gate</text>
      <rect x="235" y="45" width="70" height="45" fill="#10b981" rx="3"/>
      <text x="270" y="72" fill="#ffffff" font-size="10" font-weight="700" text-anchor="middle">Drain</text>
    </g>
    <g transform="translate(390, 75)">
      <rect width="345" height="325" fill="#161b22" stroke="#30363d" stroke-width="1.5" rx="8"/>
      <text x="172.5" y="25" fill="#34d399" font-size="13" font-weight="700" text-anchor="middle">2. Process &amp; Materials Specs</text>
      <rect x="15" y="45" width="315" height="75" fill="#0d1117" stroke="#30363d" rx="5"/>
      <text x="30" y="70" fill="#6ee7b7" font-size="11" font-weight="700">Deposition &amp; Etch Selectivity:</text>
      <text x="30" y="95" fill="#8b98a5" font-size="10">&gt; 50:1 Target Selectivity, Sub-nm Uniformity Control</text>
      <rect x="15" y="135" width="315" height="75" fill="#0d1117" stroke="#30363d" rx="5"/>
      <text x="30" y="160" fill="#6ee7b7" font-size="11" font-weight="700">Thermal &amp; Stress Budget:</text>
      <text x="30" y="185" fill="#8b98a5" font-size="10">Rapid Thermal Anneal (RTA) &lt; 1050°C, Stress Migration Low</text>
      <rect x="15" y="225" width="315" height="80" fill="#0d1117" stroke="#047857" rx="5"/>
      <text x="30" y="250" fill="#3fb950" font-size="11" font-weight="700">Yield &amp; Defect Metric:</text>
      <text x="30" y="275" fill="#e6edf3" font-size="10">Critical Dimension (CD) Variation &lt; 1.2%, D0 Defect &lt; 0.05/cm²</text>
    </g>
  
  <!-- Key insight bar -->
  <rect x="25" y="415" width="710" height="22" rx="3" fill="#0b1220" stroke="#233043" stroke-width="0.8"/>
  <text x="380" y="430" fill="#fbbf24" font-size="9" font-weight="700" text-anchor="middle">Key Insight: Optimal Physics Based Modeling And Differential Equations architecture balances performance throughput, systemic latency, and physical constraints.</text>
  
  <text x="380" y="460" fill="#6b7684" font-size="11" text-anchor="middle">Technical specification &amp; verification reference for Physics Based Modeling And Differential Equations (Row ID 10702)</text>
</svg>

A comprehensive reference for the physics and mathematics governing semiconductor fabrication processes.

Thermal Oxidation of Silicon

Deal-Grove Model

The foundational model for silicon oxidation describes oxide thickness growth through coupled transport and reaction.

Governing Equation:

Parameter Definitions:

$x$ — oxide thickness $A = \frac{2D_{ox}}{k_s}$ — linear rate constant parameter (related to surface reaction) $B = \frac{2D_{ox}C^}{N_1}$ — parabolic rate constant (related to diffusion) $D_{ox}$ — oxidant diffusivity through oxide $k_s$ — surface reaction rate constant $C^$ — equilibrium oxidant concentration at gas-oxide interface $N_1$ — number of oxidant molecules incorporated per unit volume of oxide $\tau$ — time shift accounting for initial oxide

Underlying Diffusion Physics

Steady-state diffusion through the oxide:

Boundary Conditions:

Gas-oxide interface (flux from gas phase):

Si-SiO₂ interface (surface reaction):

Steady-state flux through the oxide:

Limiting Growth Regimes

Dopant Diffusion

Fick's Laws of Diffusion

First Law (Flux Equation):

Second Law (Mass Conservation / Continuity):

For constant diffusivity in 1D:

Analytical Solutions

Constant Surface Concentration (Predeposition)

Initial condition: $C(x, 0) = 0$ Boundary condition: $C(0, t) = C_s$

where the complementary error function is:

Fixed Dose / Drive-in (Gaussian Distribution)

Initial condition: Delta function at surface with dose $Q$

Key Parameters:

$Q$ — total dose per unit area (atoms/cm²) $\sqrt{Dt}$ — diffusion length Peak concentration: $C_{max} = \frac{Q}{\sqrt{\pi Dt}}$

Concentration-Dependent Diffusion

At high doping concentrations, diffusivity becomes concentration-dependent:

Fair-Tsai Model for Diffusivity:

Parameter Definitions:

$D_i$ — intrinsic diffusivity (via neutral defects) $D^-$ — diffusivity via negatively charged defects $D^+$ — diffusivity via singly positive charged defects $D^{++}$ — diffusivity via doubly positive charged defects $n, p$ — electron and hole concentrations $n_i$ — intrinsic carrier concentration

Point Defect Coupled Diffusion

Modern TCAD uses coupled equations for dopants and point defects (vacancies $V$ and interstitials $I$):

Vacancy Continuity:

Interstitial Continuity:

Term Definitions:

$D_V, D_I$ — diffusion coefficients for vacancies and interstitials $k_{IV}$ — recombination rate constant for $V$-$I$ annihilation $G_V, G_I$ — generation rates $C_V^, C_I^$ — equilibrium concentrations $\tau_V, \tau_I$ — lifetimes at sinks (surfaces, dislocations)

Effective Dopant Diffusivity:

where $f_I$ and $f_V$ are the interstitial and vacancy fractions for the specific dopant species.

Ion Implantation

Range Distribution (LSS Theory)

The implanted dopant profile follows approximately a Gaussian distribution:

Parameters:

$\Phi$ — dose (ions/cm²) $R_p$ — projected range (mean implant depth) $\Delta R_p$ — straggle (standard deviation of range distribution)

Higher-Order Moments (Pearson IV Distribution):

$\gamma$ — skewness (asymmetry) $\beta$ — kurtosis (peakedness)

Stopping Power (Energy Loss)

The rate of energy loss as ions traverse the target:

Components:

$S_n(E)$ — nuclear stopping power (elastic collisions with target nuclei) $S_e(E)$ — electronic stopping power (inelastic interactions with electrons) $N$ — atomic density of target material (atoms/cm³)

LSS Electronic Stopping (Low Energy):

Nuclear Stopping: Uses screened Coulomb potentials with Thomas-Fermi or ZBL (Ziegler-Biersack-Littmark) universal screening functions.

Boltzmann Transport Equation

For rigorous treatment (typically solved via Monte Carlo methods):

Variables:

$f(\vec{r}, \vec{v}, t)$ — particle distribution function $\vec{F}$ — external force Right-hand side — collision integral

Damage Accumulation

Kinchin-Pease Model:

Parameters:

$N_d$ — number of displaced atoms $E_{damage}$ — energy available for displacement $E_d$ — displacement threshold energy ($\approx 15$ eV for silicon)

Chemical Vapor Deposition (CVD)

Coupled Transport Equations

Species Transport (Convection-Diffusion-Reaction):

Navier-Stokes Equations (Momentum):

Continuity Equation (Incompressible Flow):

Energy Equation:

Variable Definitions:

$C_i$ — concentration of species $i$ $\vec{u}$ — velocity vector $D_i$ — diffusion coefficient of species $i$ $R_i$ — net reaction rate for species $i$ $\rho$ — density $p$ — pressure $\mu$ — dynamic viscosity $c_p$ — specific heat at constant pressure $k$ — thermal conductivity $Q_{reaction}$ — heat of reaction

Surface Reaction Kinetics

Flux Balance at Wafer Surface:

Deposition Rate:

Parameters:

$h_m$ — mass transfer coefficient $k_s$ — surface reaction rate constant $C_b$ — bulk gas concentration $C_s$ — surface concentration

Limiting Cases:

Step Coverage — Knudsen Diffusion

In high-aspect-ratio features, molecular (Knudsen) flow dominates:

Parameters:

$d$ — characteristic feature dimension $k_B$ — Boltzmann constant $T$ — temperature $m$ — molecular mass

Thiele Modulus (Reaction-Diffusion Balance):

Interpretation:

$\phi \ll 1$ — Reaction-limited → Conformal deposition $\phi \gg 1$ — Diffusion-limited → Poor step coverage

Atomic Layer Deposition (ALD)

Surface Site Model

Precursor A Adsorption Kinetics:

Parameters:

$\theta_A$ — fractional surface coverage of precursor A $s_0$ — sticking coefficient $P_A$ — partial pressure of precursor A $m_A$ — molecular mass of precursor A $k_{des}$ — desorption rate constant

Growth Per Cycle (GPC)

Parameters:

$n_{sites}$ — surface site density (sites/cm²) $\Omega$ — atomic volume (volume per deposited atom) $\theta_A^{sat}$ — saturation coverage achieved during half-cycle

Plasma Etching

Plasma Fluid Equations

Electron Continuity:

Ion Continuity:

Drift-Diffusion Flux (Electrons):

Drift-Diffusion Flux (Ions):

Poisson's Equation (Self-Consistent Field):

Electron Energy Balance:

Sheath Physics

Bohm Criterion (Sheath Edge Condition):

Child-Langmuir Law (Collisionless Sheath Ion Current):

Parameters:

$u_i$ — ion velocity at sheath edge $u_B$ — Bohm velocity $T_e$ — electron temperature $M_i$ — ion mass $V_0$ — sheath voltage drop $d$ — sheath thickness

Surface Etch Kinetics

Ion-Enhanced Etching Rate:

Components:

$Y_i\Gamma_i$ — physical sputtering contribution $Y_n\Gamma_n(1-\theta)$ — spontaneous chemical etching $Y_{syn}\Gamma_i\theta$ — ion-enhanced (synergistic) etching

Yield Parameters:

$Y_i$ — physical sputtering yield $Y_n$ — spontaneous chemical etch yield $Y_{syn}$ — synergistic yield (ion-enhanced chemistry) $\Gamma_i, \Gamma_n$ — ion and neutral fluxes $\theta$ — fractional surface coverage of reactive species

Surface Coverage Dynamics:

Terms:

$s\Gamma_n(1-\theta)$ — adsorption onto empty sites $Y_{syn}\Gamma_i\theta$ — consumption by ion-enhanced reaction $k_v\theta$ — thermal desorption/volatilization

Lithography

Aerial Image Formation

Hopkins Formulation (Partially Coherent Imaging):

Parameters:

$TCC$ — Transmission Cross Coefficient (encapsulates partial coherence) $\tilde{M}(f,g)$ — Fourier transform of mask transmission function $f, g$ — spatial frequencies

Rayleigh Resolution Criterion:

Depth of Focus:

Parameters:

$k_1, k_2$ — process-dependent factors $\lambda$ — exposure wavelength $NA$ — numerical aperture

Photoresist Exposure — Dill Model

Intensity Attenuation with Photobleaching:

where the absorption coefficient depends on PAC concentration:

Photoactive Compound (PAC) Decomposition:

Dill Parameters:

Chemically Amplified Resists

Photoacid Generation:

Post-Exposure Bake — Acid Diffusion and Reaction:

Deprotection Reaction (Catalytic Amplification):

Parameters:

$[PAG]$ — photoacid generator concentration $D_{acid}$ — acid diffusion coefficient $k_{loss}$ — acid loss rate (neutralization, evaporation) $k_{cat}$ — catalytic deprotection rate constant

Development Rate — Mack Model

Parameters:

$R_{max}$ — maximum development rate (fully exposed) $R_{min}$ — minimum development rate (unexposed) $a$ — selectivity parameter $n$ — contrast parameter $M$ — normalized PAC concentration after exposure

Epitaxy

Burton-Cabrera-Frank (BCF) Theory

Adatom Diffusion on Terraces:

Parameters:

$n$ — adatom density on terrace $D_s$ — surface diffusion coefficient $F$ — deposition flux (atoms/cm²·s) $\tau$ — adatom lifetime before desorption

Step Velocity:

Steady-State Solution for Step Flow:

Parameters:

$\Omega$ — atomic volume $\lambda_s = \sqrt{D_s \tau}$ — surface diffusion length $l$ — terrace width

Rate Equations for Island Nucleation

Monomer (Single Adatom) Density:

Cluster of Size $j$:

Parameters:

$n_j$ — density of clusters containing $j$ atoms $\sigma_j$ — capture cross-section for clusters of size $j$

Chemical Mechanical Polishing (CMP)

Preston Equation

Parameters:

$MRR$ — material removal rate (nm/min) $K_p$ — Preston coefficient (material/process dependent) $P$ — applied pressure $V$ — relative velocity between pad and wafer

Contact Mechanics — Greenwood-Williamson Model

Real Contact Area:

Parameters:

$\eta$ — asperity density $A_n$ — nominal contact area $R_p$ — asperity radius $d$ — separation distance $\phi(z)$ — asperity height distribution

Slurry Hydrodynamics — Reynolds Equation

Parameters:

$h$ — film thickness $p$ — pressure $\mu$ — dynamic viscosity $U$ — sliding velocity

Thin Film Stress

Stoney Equation

Film Stress from Wafer Curvature:

Parameters:

$\sigma_f$ — film stress $E_s$ — substrate Young's modulus $u_s$ — substrate Poisson's ratio $h_s$ — substrate thickness $h_f$ — film thickness $R$ — radius of curvature

Thermal Stress

Parameters:

$E_f$ — film Young's modulus $u_f$ — film Poisson's ratio $\alpha_s, \alpha_f$ — thermal expansion coefficients (substrate, film) $\Delta T$ — temperature change from deposition

Electromigration (Reliability)

Black's Equation (Empirical MTTF)

Parameters:

$MTTF$ — mean time to failure $j$ — current density $n$ — current density exponent (typically 1-2) $E_a$ — activation energy $A$ — material/geometry constant

Drift-Diffusion Model

Parameters:

$C$ — atomic concentration $D$ — diffusion coefficient $Z^*$ — effective charge number (wind force parameter) $\rho$ — electrical resistivity $\vec{j}$ — current density vector

Stress Evolution — Korhonen Model

Parameters:

$\sigma$ — hydrostatic stress $D_a$ — atomic diffusivity $B$ — effective bulk modulus $\Omega$ — atomic volume

Numerical Solution Methods

Common Numerical Techniques

Discretization Examples

Explicit Forward Euler (1D Diffusion):

Stability Criterion:

Implicit Backward Euler:

Major TCAD Software Tools

Synopsys Sentaurus — comprehensive process and device simulation Silvaco ATHENA/ATLAS — process and device modeling COMSOL Multiphysics — general multiphysics platform SRIM/TRIM — ion implantation Monte Carlo PROLITH — lithography simulation

Processes and Governing Equations

A physics model should begin with the quantity of interest and a falsifiable claim. Predicting terminal current, junction temperature, wafer-scale dopant uniformity, trench profile, residual stress, or defect density requires different domains, state variables, and closures. The quantity of interest determines acceptable spatial and temporal resolution, experimental evidence, and uncertainty. State the operating envelope and decision before choosing equations. A model calibrated to one current-voltage curve does not automatically predict self-heating or breakdown, and a process model matching average depth does not automatically predict lateral profile. Model scope is a scientific claim, not a software feature list.

Conservation laws connect a control volume to a differential equation. For a conserved density $u$, accumulation equals net boundary flux plus volumetric source: $\frac{d}{dt}\int_Vu\,dV=-\int_{\partial V}\mathbf{J}\cdot\mathbf{n}\,dA+\int_VS\,dV$. Applying the divergence theorem gives $\partial_tu+\nabla\cdot\mathbf{J}=S$. Charge, particles, chemical species, mass, energy, and momentum share this structure even though their fluxes and sources differ. Writing the integral form first exposes units, signs, boundary exchange, and conservation checks. A discretization should reproduce the corresponding global balance.

Constitutive laws close conservation equations with material physics. A balance equation alone does not specify flux. Fourier conduction uses $\mathbf{q}=-k\nabla T$; Fickian diffusion uses $\mathbf{J}=-D\nabla c$; carrier drift-diffusion combines field-driven and concentration-driven terms; elasticity relates stress and strain; reaction models define sources. Coefficients may depend on temperature, field, concentration, crystal direction, stress, phase, and history. Every closure has a validity range. Treating mobility, thermal conductivity, diffusivity, or reaction rate as a universal constant can shift error into a fitted boundary condition.

Scale analysis determines which physics can be neglected. Compare characteristic time, length, field, velocity, and energy scales before solving. Debye length indicates electrostatic screening, diffusion length relates transport to recombination, mean free path tests continuum assumptions, thermal diffusion time tests quasi-steady heat flow, and dimensionless groups compare convection, diffusion, reaction, or inertia. A term small in the bulk may dominate in a thin interface. Nondimensionalization reveals singular perturbations and improves numerical scaling. Assumptions such as quasineutrality, isothermal operation, local equilibrium, or steady state should follow from these ratios.

~~~svg ~~~

PDE classification predicts information flow and numerical behavior. Elliptic equations such as steady Poisson problems communicate boundary influence across the domain. Parabolic equations such as diffusion and heat conduction smooth initial data while evolving in time. Hyperbolic equations propagate finite-speed waves and require attention to characteristics and upwinding. Coupled semiconductor systems can mix these types and become strongly nonlinear or degenerate. Classification guides the number and placement of boundary conditions, timestep restrictions, stabilization, and solver choice. Calling every spatial equation a diffusion equation hides important mathematical differences.

Initial, boundary, and interface conditions are part of the physical model. Dirichlet data prescribe a state, Neumann data prescribe flux, Robin data couple state to exchange, and dynamic boundaries carry their own storage. Semiconductor contacts may impose electrochemical potentials, recombination velocities, tunneling currents, or circuit relations. Material interfaces enforce appropriate flux continuity and may include sheet charge, thermal resistance, segregation, reaction, or mechanics. Initial conditions must satisfy constraints closely enough for the intended transient. Boundary data inferred from a tool setting often require a separate transfer model.

Poisson’s equation links electrostatic potential to charge. A common semiconductor form is $-\nabla\cdot(\epsilon\nabla\psi)=q(p-n+N_D^+-N_A^-+\rho_t/q)$, where permittivity and charge populations can vary in space and state. Electric field $\mathbf{E}=-\nabla\psi$ drives transport and feeds back through carrier distributions. Interface charge, incomplete ionization, traps, polarization, and quantum corrections modify the source or closure. Gauge or reference potential must be fixed. Global charge and terminal displacement current provide useful consistency checks.

Carrier continuity expresses generation, recombination, and current divergence. Electron and hole equations take forms such as $\partial_t n=(1/q)\nabla\cdot\mathbf{J}_n+G-R$ and $\partial_t p=-(1/q)\nabla\cdot\mathbf{J}_p+G-R$, subject to the chosen sign convention. Integrating over the device connects contact currents, stored charge, and net generation-recombination. Van Roosbroeck’s 1950 formulation established the coupled electrostatic, drift, diffusion, and recombination structure that remains central to device simulation. Numerical residuals should be accompanied by terminal-current and charge-balance checks.

Drift-diffusion is a local near-equilibrium transport closure. For nondegenerate statistics, electron current may be written $\mathbf{J}_n=q\mu_nn\mathbf{E}+qD_n\nabla n$, with a corresponding hole expression and Einstein relation under its assumptions. Quasi-Fermi potentials often provide better variables because current follows their gradients and equilibrium is represented naturally. Mobility can depend on doping, field, temperature, interfaces, and stress. At nanoscale or high field, energy transport, hydrodynamic, Monte Carlo, Boltzmann, ballistic, or quantum models may be required. More advanced physics should be justified by a failed observable, not fashion.

~~~svg ~~~

Generation and recombination models encode distinct microscopic mechanisms. Shockley–Read–Hall recombination depends on trap energy, capture parameters, and carrier populations; Auger processes grow strongly at high injection; radiative recombination matters in direct-gap materials; impact ionization drives avalanche; optical generation depends on absorption and field distribution. Summing named rates is not enough if parameters are unidentifiable or double-counted. Temperature and field dependencies must be consistent. Lifetime measured in one structure may include surface effects that should not become a bulk constant in another.

Contact models often control the result more than the bulk equations. Ohmic contacts may impose carrier populations or quasi-Fermi levels, while Schottky contacts require barrier, image-force, thermionic, tunneling, and interface-state considerations. Series resistance, current crowding, contact geometry, and circuit loading can alter terminal data. Thermal boundaries at contacts also control self-heating. Calibrating bulk mobility against contact-limited current creates a nonportable parameter. Use geometry or temperature splits that distinguish contact from channel and validate internal profiles when possible.

Scharfetter–Gummel fluxes stabilize drift-dominated carrier transport. Naive centered differences can oscillate when electrostatic drift overwhelms diffusion. Scharfetter and Gummel integrated a local one-dimensional flux under assumptions across a mesh edge, producing a Bernoulli-function form that preserves equilibrium structure and behaves like upwinding in strong fields. The method’s success does not remove mesh, statistics, mobility, multidimensional, or interface issues. Verify sign conventions and limiting behavior for small potential differences. Flux continuity and positivity deserve explicit tests.

Nonlinear coupling can be solved segregated or monolithically. Gummel iteration alternates Poisson and carrier equations, often with damping, while Newton methods assemble a coupled Jacobian and can converge rapidly near a solution. Segregated schemes are modular but may stall under strong coupling; monolithic schemes are robust in some regimes but demand accurate derivatives, scaling, and linear algebra. Continuation in voltage, generation, temperature, or model complexity helps trace difficult branches. Convergence should be judged in scaled residuals and physical balances, not only update size.

Electrothermal modeling closes power and temperature feedback. Heat conduction with storage can be written $\rho c_p\partial_tT-\nabla\cdot(k\nabla T)=Q$, where $Q$ may include Joule heating, recombination, optical absorption, reactions, and thermoelectric terms. Temperature changes mobility, band structure, leakage, reaction rates, and boundary exchange, creating feedback and possibly thermal runaway. Thermal boundary resistance, package spreading, and pulsed duty cycle are often more uncertain than bulk conductivity. Validate temperature with a measurement model that represents spatial and temporal averaging.

~~~svg ~~~

Mechanical equations connect process history to stress and deformation. Small-strain equilibrium uses $\nabla\cdot\boldsymbol{\sigma}+\mathbf{b}=0$ with a constitutive law relating stress to elastic, thermal, plastic, creep, transformation, and eigenstrains. Thin films carry intrinsic stress; thermal expansion mismatch bends wafers; oxidation or phase change creates volume strain; stress alters mobility and band structure. Geometry may evolve enough to require nonlinear mechanics or contact. Curvature alone cannot uniquely identify a depth-varying stress field, so calibration needs appropriate observables.

Process transport couples species diffusion, reaction, and moving material state. Dopant diffusion may depend on charged point defects, clustering, activation, stress, and concentration. Oxidation couples oxidant transport with interface reaction and volume expansion. Deposition and etching combine gas or plasma transport, surface coverage, reaction, and profile motion. A reaction-diffusion equation $\partial_tc=\nabla\cdot(D\nabla c)+R(c,T,\ldots)$ is only the starting skeleton. Material labels, interfaces, and history variables determine which coefficients and sources apply at each location.

Moving boundaries require both a velocity law and conservative geometry update. Level-set, phase-field, arbitrary Lagrangian–Eulerian, volume-of-fluid, and front-tracking methods represent interfaces differently. Normal speed may follow incorporated flux, etch yield, curvature, stress, or local chemistry. Reinitialization, remeshing, and topology changes can add or remove material numerically. Verify a planar analytical case, volume conservation, symmetry, and grid convergence before trusting a complex trench or oxidation front. The interface condition is the physics; the geometric method is its numerical carrier.

Plasma process models introduce kinetic and electromagnetic closures. Electron-impact rates depend on the electron energy distribution, while charged-species transport couples to electric fields and sheath boundaries. Global models, fluid models, hybrid methods, particle-in-cell simulations, and Boltzmann solvers answer different questions. Surface ion energy-angle distributions may matter more than volume density for profile evolution. Quasineutral bulk assumptions fail in sheaths. Couple plasma outputs to feature models with resolved flux, energy, angle, time, and uncertainty rather than a generic power setting.

~~~svg ~~~

Nondimensional variables improve interpretation and solver conditioning. Choose reference scales so principal unknowns and residuals are comparable. The Péclet number compares convection with diffusion, Damköhler number reaction with transport, Fourier number transient diffusion time, Biot number internal with boundary thermal resistance, and Debye-scaled ratios electrostatic length scales. Distinct definitions apply at reactor, wafer, device, and feature scales. A nondimensional equation exposes dominant terms and parameters; numerical scaling then prevents one equation’s units from overwhelming a coupled Newton norm.

Weak and integral forms map naturally to different discretizations. Finite volume methods conserve flux locally by construction, finite elements derive a weak form and handle complex geometry and multiphysics flexibly, finite differences are efficient on structured domains, and spectral methods can converge rapidly for smooth solutions. Discontinuous Galerkin and mixed methods offer additional conservation or stability properties. Method choice should follow conservation, regularity, geometry, discontinuities, and quantities of interest. Software convenience is not a numerical analysis.

Mesh refinement must target the quantity of interest. Resolve depletion regions, interfaces, boundary layers, high-field zones, steep thermal gradients, and moving fronts. A small global element size does not guarantee adequate anisotropic resolution. Compare systematic refinements and estimate observed order where the solution is smooth; use goal-oriented indicators when terminal current or peak temperature matters more than a field norm. Mesh adaptation must not erase conservation or move material inconsistently. Report the mesh and convergence of the actual decision metric.

Transient accuracy requires more than a stable timestep. Explicit schemes may face diffusion or wave stability limits; implicit schemes allow larger steps but can smear fast events or converge to an inaccurate trajectory. Stiff reactions and widely separated carrier, thermal, and process timescales motivate adaptive implicit integration and consistent error control. Resolve input edges, circuit dynamics, trapping, nucleation, or pulses relevant to the quantity of interest. Timestep convergence should be independent of nonlinear tolerances. Event handling must preserve state across discontinuous boundary changes.

Code verification asks whether the equations were implemented correctly. Analytical solutions, method of manufactured solutions, symmetry, conservation identities, and independent implementations test coding and discretization. A manufactured source creates a chosen exact solution, allowing observed spatial or temporal convergence even when the manufactured field is not a physical experiment. NIST guidance distinguishes code verification from claims about physical reality. Unit tests should cover constitutive limits, Jacobians, boundary conditions, interface fluxes, and parameter units before system validation begins.

~~~svg ~~~

Solution verification estimates numerical error in the reported calculation. Mesh, timestep, iterative tolerance, domain truncation, quadrature, and stochastic sampling each contribute. Demonstrating residual convergence on one mesh is insufficient. Repeat the calculation at controlled resolutions, check conservation, inspect local indicators, and estimate uncertainty in the quantity of interest. Nonlinear multiple solutions and hysteresis require continuation and initialization studies. Solver failure is obvious; silent numerical diffusion or premature convergence is more dangerous.

Validation asks whether the model represents the intended physical system. Compare predictions with measurements not used to fit parameters, including measurement uncertainty and the instrument’s transfer function. Validation is conditional on quantity, operating range, geometry, and material state. NIST summarizes the distinction as verification checking mathematical implementation while validation checks physical representation. A model may be valid for terminal current but not internal temperature or breakdown. Predefine metrics and acceptance criteria, and preserve failed cases as evidence for model-form revision.

Calibration must respect identifiability and experimental design. Mobility, lifetime, contact resistance, interface charge, thermal boundary resistance, and geometry can compensate in a fit. Sensitivity matrices, profile likelihoods, posterior correlations, or Fisher information reveal which combinations the data constrain. Use experiments that perturb competing mechanisms differently and reserve independent validation conditions. Priors and regularization can stabilize inference but do not create information. Report parameter uncertainty and correlation rather than a single fitted deck.

Sensitivity and uncertainty should guide the next measurement. Local derivatives reveal nearby response; global variance methods reveal interactions across an operating envelope; adjoints efficiently differentiate many parameters for few outputs. Propagate uncertainty in geometry, material data, boundary conditions, calibration, numerical error, and model discrepancy to prediction intervals. Rank contributions to determine whether better metrology, a new experiment, finer mesh, or improved physics has the most value. A deterministic best-fit contour cannot support a risk-based manufacturing decision.

Reduced-order and learned models inherit the high-fidelity validity envelope. Surrogates, proper orthogonal decomposition, Gaussian processes, neural operators, and physics-informed networks can accelerate optimization or control. They must preserve key conservation and boundary behavior, quantify interpolation uncertainty, and detect departure from training conditions. Randomly splitting nearby simulation points overstates extrapolation skill. Compare against withheld geometries or regimes and fall back to the verified solver when outside the trusted domain. Speed is valuable only after the source model and data are credible.

~~~svg ~~~

start: Define quantity of interest operating envelope and decision
scale: Choose domain state variables scales and neglected physics
laws: Derive integral conservation balances
closure: Add constitutive source boundary and interface laws
regime: Nondimensionalize and classify the coupled PDE system
discretize: Select conservative stable spatial and temporal methods
verify: Verify units limits Jacobians balances and manufactured solutions
solution: Quantify mesh timestep iteration and domain error
calibrate: Estimate only identifiable uncertain parameters
validate: Predict independent measurements through an instrument model
accept: Are metrics within predefined uncertainty-aware limits?
deploy: Report validity envelope and propagate prediction uncertainty
revise: Classify residuals and replace the falsified assumption
start->scale->laws->closure->regime->discretize->verify->solution->calibrate->validate->accept
accept->deploy
accept->revise
revise->scale

A physics-based model earns trust by surviving a prediction it was not tuned to reproduce. After verification and calibration, specify the sign, magnitude, and uncertainty of a response to a new bias, geometry, temperature, material state, or process condition before measuring it. Success supports the claimed envelope; failure identifies a closure, boundary, parameter, or scale assumption to revise. Read physics-based modeling and differential equations through a conservation-and-credibility lens rather than an equation-collection lens.