Chip simulation
**Device Physics & Mathematical Modeling**
1. Fundamental Mathematical Structure
Semiconductor modeling is built on coupled nonlinear partial differential equations spanning multiple scales:
| Scale | Methods | Typical Equations |
|:------|:--------|:------------------|
| Quantum (< 1 nm) | DFT, Schrödinger | $H\psi = E\psi$ |
| Atomistic (1–100 nm) | MD, Kinetic Monte Carlo | Newton's equations, master equations |
| Continuum (nm–mm) | Drift-diffusion, FEM | PDEs (Poisson, continuity, heat) |
| Circuit | SPICE | ODEs, compact models |
Multiscale Hierarchy
The mathematics forms a hierarchy of models through successive averaging:
$$
\boxed{\text{Schrödinger} \xrightarrow{\text{averaging}} \text{Boltzmann} \xrightarrow{\text{moments}} \text{Drift-Diffusion} \xrightarrow{\text{fitting}} \text{Compact Models}}
$$
2. Process Physics & Models
2.1 Oxidation: Deal-Grove Model
Thermal oxidation of silicon follows linear-parabolic kinetics :
$$
\frac{dx_{ox}}{dt} = \frac{B}{A + 2x_{ox}}
$$
where:
- $x_{ox}$ = oxide thickness
- $B/A$ = linear rate constant (surface-reaction limited)
- $B$ = parabolic rate constant (diffusion limited)
Limiting Cases:
- Thin oxide (reaction-limited):
$$
x_{ox} \approx \frac{B}{A} \cdot t
$$
- Thick oxide (diffusion-limited):
$$
x_{ox} \approx \sqrt{B \cdot t}
$$
Physical Mechanism:
1. O₂ transport from gas to oxide surface
2. O₂ diffusion through growing SiO₂ layer
3. Reaction at Si/SiO₂ interface: $\text{Si} + \text{O}_2 \rightarrow \text{SiO}_2$
> Note: This is a Stefan problem (moving boundary PDE).
2.2 Diffusion: Fick's Laws
Dopant redistribution follows Fick's second law :
$$
\frac{\partial C}{\partial t} =
abla \cdot \left( D(C, T)
abla C \right)
$$
For constant $D$ in 1D:
$$
\frac{\partial C}{\partial t} = D \frac{\partial^2 C}{\partial x^2}
$$
Analytical Solutions (1D, constant D):
- Constant surface concentration (infinite source):
$$
C(x,t) = C_s \cdot \text{erfc}\left( \frac{x}{2\sqrt{Dt}} \right)
$$
- Limited source (e.g., implant drive-in):
$$
C(x,t) = \frac{Q}{\sqrt{\pi D t}} \exp\left( -\frac{x^2}{4Dt} \right)
$$
where $Q$ = dose (atoms/cm²)
Complications at High Concentrations:
- Concentration-dependent diffusivity: $D = D(C)$
- Electric field effects: Charged point defects create internal fields
- Vacancy/interstitial mechanisms: Different diffusion pathways
$$
\frac{\partial C}{\partial t} = \frac{\partial}{\partial x}\left[ D(C) \frac{\partial C}{\partial x} \right] + \mu C \frac{\partial \phi}{\partial x}
$$
2.3 Ion Implantation: Range Theory
The implanted dopant profile is approximately Gaussian :
$$
C(x) = \frac{\Phi}{\sqrt{2\pi} \Delta R_p} \exp\left( -\frac{(x - R_p)^2}{2 (\Delta R_p)^2} \right)
$$
where:
- $\Phi$ = implant dose (ions/cm²)
- $R_p$ = projected range (mean depth)
- $\Delta R_p$ = straggle (standard deviation)
LSS Theory (Lindhard-Scharff-Schiøtt) predicts stopping power:
$$
-\frac{dE}{dx} = N \left[ S_n(E) + S_e(E) \right]
$$
where:
- $S_n(E)$ = nuclear stopping power (dominant at low energy)
- $S_e(E)$ = electronic stopping power (dominant at high energy)
- $N$ = target atomic density
For asymmetric profiles , the Pearson IV distribution is used:
$$
C(x) = \frac{\Phi \cdot K}{\Delta R_p} \left[ 1 + \left( \frac{x - R_p}{a} \right)^2 \right]^{-m} \exp\left[ -
u \arctan\left( \frac{x - R_p}{a} \right) \right]
$$
> Modern approach: Monte Carlo codes (SRIM/TRIM) for accurate profiles including channeling effects.
2.4 Lithography: Optical Imaging
Aerial image formation follows Hopkins' partially coherent imaging theory :
$$
I(\mathbf{r}) = \iint TCC(f, f') \cdot \tilde{M}(f) \cdot \tilde{M}^*(f') \cdot e^{2\pi i (f - f') \cdot \mathbf{r}} \, df \, df'
$$
where:
- $TCC$ = Transmission Cross-Coefficient
- $\tilde{M}(f)$ = mask spectrum (Fourier transform of mask pattern)
- $\mathbf{r}$ = position in image plane
Fundamental Limits:
- Rayleigh resolution criterion:
$$
CD_{\min} = k_1 \frac{\lambda}{NA}
$$
- Depth of focus:
$$
DOF = k_2 \frac{\lambda}{NA^2}
$$
where:
- $\lambda$ = wavelength (193 nm for ArF, 13.5 nm for EUV)
- $NA$ = numerical aperture
- $k_1, k_2$ = process-dependent factors
Resist Modeling — Dill Equations:
$$
\frac{\partial M}{\partial t} = -C \cdot I(z) \cdot M
$$
$$
\frac{dI}{dz} = -(\alpha M + \beta) I
$$
where $M$ = photoactive compound concentration.
2.5 Etching & Deposition: Surface Evolution
Topography evolution is modeled with the level set method :
$$
\frac{\partial \phi}{\partial t} + V |
abla \phi| = 0
$$
where:
- $\phi(\mathbf{r}, t) = 0$ defines the surface
- $V$ = local velocity (etch rate or deposition rate)
For anisotropic etching:
$$
V = V(\theta, \phi, \text{ion flux}, \text{chemistry})
$$
CVD in High Aspect Ratio Features:
Knudsen diffusion limits step coverage:
$$
\frac{\partial C}{\partial t} = D_K
abla^2 C - k_s C \cdot \delta_{\text{surface}}
$$
where:
- $D_K = \frac{d}{3}\sqrt{\frac{8k_BT}{\pi m}}$ (Knudsen diffusivity)
- $d$ = feature width
- $k_s$ = surface reaction rate
ALD (Atomic Layer Deposition):
Self-limiting surface reactions follow Langmuir kinetics:
$$
\theta = \frac{K \cdot P}{1 + K \cdot P}
$$
where $\theta$ = surface coverage, $P$ = precursor partial pressure.
3. Device Physics: Semiconductor Equations
The core mathematical framework for device simulation consists of three coupled PDEs :
3.1 Poisson's Equation (Electrostatics)
$$
abla \cdot (\varepsilon
abla \psi) = -q \left( p - n + N_D^+ - N_A^- \right)
$$
where:
- $\psi$ = electrostatic potential
- $n, p$ = electron and hole concentrations
- $N_D^+, N_A^-$ = ionized donor and acceptor concentrations
3.2 Continuity Equations (Carrier Conservation)
Electrons:
$$
\frac{\partial n}{\partial t} = \frac{1}{q}
abla \cdot \mathbf{J}_n + G - R
$$
Holes:
$$
\frac{\partial p}{\partial t} = -\frac{1}{q}
abla \cdot \mathbf{J}_p + G - R
$$
where:
- $G$ = generation rate
- $R$ = recombination rate
3.3 Current Density Equations (Transport)
Drift-Diffusion Model:
$$
\mathbf{J}_n = q \mu_n n \mathbf{E} + q D_n
abla n
$$
$$
\mathbf{J}_p = q \mu_p p \mathbf{E} - q D_p
abla p
$$
Einstein Relation:
$$
\frac{D_n}{\mu_n} = \frac{D_p}{\mu_p} = \frac{k_B T}{q} = V_T
$$
3.4 Recombination Models
Shockley-Read-Hall (SRH) Recombination:
$$
R_{SRH} = \frac{np - n_i^2}{\tau_p (n + n_1) + \tau_n (p + p_1)}
$$
Auger Recombination:
$$
R_{Auger} = C_n n (np - n_i^2) + C_p p (np - n_i^2)
$$
Radiative Recombination:
$$
R_{rad} = B (np - n_i^2)
$$
3.5 MOSFET Physics
Threshold Voltage:
$$
V_T = V_{FB} + 2\phi_B + \frac{\sqrt{2 \varepsilon_{Si} q N_A (2\phi_B)}}{C_{ox}}
$$
where:
- $V_{FB}$ = flat-band voltage
- $\phi_B = \frac{k_BT}{q} \ln\left(\frac{N_A}{n_i}\right)$ = bulk potential
- $C_{ox} = \frac{\varepsilon_{ox}}{t_{ox}}$ = oxide capacitance
Drain Current (Gradual Channel Approximation):
- Linear region ($V_{DS} < V_{GS} - V_T$):
$$
I_D = \frac{W}{L} \mu_n C_{ox} \left[ (V_{GS} - V_T) V_{DS} - \frac{V_{DS}^2}{2} \right]
$$
- Saturation region ($V_{DS} \geq V_{GS} - V_T$):
$$
I_D = \frac{W}{2L} \mu_n C_{ox} (V_{GS} - V_T)^2
$$
4. Quantum Effects at Nanoscale
For modern devices with gate lengths $L_g < 10$ nm, classical models fail.
4.1 Quantum Confinement
In thin silicon channels, carrier energy becomes quantized :
$$
E_n = \frac{\hbar^2 \pi^2 n^2}{2 m^* t_{Si}^2}
$$
where:
- $n$ = quantum number (1, 2, 3, ...)
- $m^*$ = effective mass
- $t_{Si}$ = silicon body thickness
Effects:
- Increased threshold voltage
- Modified density of states: $g_{2D}(E) = \frac{m^*}{\pi \hbar^2}$ (step function)
4.2 Quantum Tunneling
Gate Leakage (Direct Tunneling):
WKB approximation:
$$
T \approx \exp\left( -2 \int_0^{t_{ox}} \kappa(x) \, dx \right)
$$
where $\kappa = \sqrt{\frac{2m^*(\Phi_B - E)}{\hbar^2}}$
Source-Drain Tunneling:
Limits OFF-state current in ultra-short channels.
Band-to-Band Tunneling:
Enables Tunnel FETs (TFETs):
$$
I_{BTBT} \propto \exp\left( -\frac{4\sqrt{2m^*} E_g^{3/2}}{3q\hbar |\mathbf{E}|} \right)
$$
4.3 Ballistic Transport
When channel length $L < \lambda_{mfp}$ (mean free path), the Landauer formalism applies:
$$
I = \frac{2q}{h} \int T(E) \left[ f_S(E) - f_D(E) \right] dE
$$
where:
- $T(E)$ = transmission probability
- $f_S, f_D$ = source and drain Fermi functions
Ballistic Conductance Quantum:
$$
G_0 = \frac{2q^2}{h} \approx 77.5 \, \mu\text{S}
$$
4.4 NEGF Formalism
The Non-Equilibrium Green's Function method is the gold standard for quantum transport:
$$
G^R = \left[ EI - H - \Sigma_1 - \Sigma_2 \right]^{-1}
$$
where:
- $H$ = device Hamiltonian
- $\Sigma_1, \Sigma_2$ = contact self-energies
- $G^R$ = retarded Green's function
Observables:
- Electron density: $n(\mathbf{r}) = -\frac{1}{\pi} \text{Im}[G^<(\mathbf{r}, \mathbf{r}; E)]$
- Current: $I = \frac{q}{h} \text{Tr}[\Gamma_1 G^R \Gamma_2 G^A]$
5. Numerical Methods
5.1 Discretization: Scharfetter-Gummel Scheme
The drift-diffusion current requires special treatment to avoid numerical instability:
$$
J_{n,i+1/2} = \frac{q D_n}{h} \left[ n_{i+1} B\left( -\frac{\Delta \psi}{V_T} \right) - n_i B\left( \frac{\Delta \psi}{V_T} \right) \right]
$$
where the Bernoulli function is:
$$
B(x) = \frac{x}{e^x - 1}
$$
Properties:
- $B(0) = 1$
- $B(x) \to 0$ as $x \to \infty$
- $B(-x) = x + B(x)$
5.2 Solution Strategies
Gummel Iteration (Decoupled):
1. Solve Poisson for $\psi$ (fixed $n$, $p$)
2. Solve electron continuity for $n$ (fixed $\psi$, $p$)
3. Solve hole continuity for $p$ (fixed $\psi$, $n$)
4. Repeat until convergence
Newton-Raphson (Fully Coupled):
Solve the Jacobian system:
$$
\begin{pmatrix}
\frac{\partial F_\psi}{\partial \psi} & \frac{\partial F_\psi}{\partial n} & \frac{\partial F_\psi}{\partial p} \\
\frac{\partial F_n}{\partial \psi} & \frac{\partial F_n}{\partial n} & \frac{\partial F_n}{\partial p} \\
\frac{\partial F_p}{\partial \psi} & \frac{\partial F_p}{\partial n} & \frac{\partial F_p}{\partial p}
\end{pmatrix}
\begin{pmatrix}
\delta \psi \\
\delta n \\
\delta p
\end{pmatrix}
= -
\begin{pmatrix}
F_\psi \\
F_n \\
F_p
\end{pmatrix}
$$
5.3 Time Integration
Stiffness Problem:
Time scales span ~15 orders of magnitude:
| Process | Time Scale |
|:--------|:-----------|
| Carrier relaxation | ~ps |
| Thermal response | ~μs–ms |
| Dopant diffusion | min–hours |
Solution: Use implicit methods (Backward Euler, BDF).
5.4 Mesh Requirements
Debye Length Constraint:
The mesh must resolve the Debye length:
$$
\lambda_D = \sqrt{\frac{\varepsilon k_B T}{q^2 n}}
$$
For $n = 10^{18}$ cm⁻³: $\lambda_D \approx 4$ nm
Adaptive Mesh Refinement:
- Refine near junctions, interfaces, corners
- Coarsen in bulk regions
- Use Delaunay triangulation for quality
6. Compact Models for Circuit Simulation
For SPICE-level simulation, physics is abstracted into algebraic/empirical equations.
Industry Standard Models
| Model | Device | Key Features |
|:------|:-------|:-------------|
| BSIM4 | Planar MOSFET | ~300 parameters, channel length modulation |
| BSIM-CMG | FinFET | Tri-gate geometry, quantum effects |
| BSIM-GAA | Nanosheet | Stacked channels, sheet width |
| PSP | Bulk MOSFET | Surface-potential-based |
Key Physics Captured
- Short-channel effects: DIBL, $V_T$ roll-off
- Quantum corrections: Inversion layer quantization
- Mobility degradation: Surface scattering, velocity saturation
- Parasitic effects: Series resistance, overlap capacitance
- Variability: Statistical mismatch models
Threshold Voltage Variability (Pelgrom's Law)
$$
\sigma_{V_T} = \frac{A_{VT}}{\sqrt{W \cdot L}}
$$
where $A_{VT}$ is a technology-dependent constant.
7. TCAD Co-Simulation Workflow
The complete semiconductor design flow:
```svg
```
Key Challenge: Propagating variability through the entire chain:
- Line Edge Roughness (LER)
- Random Dopant Fluctuation (RDF)
- Work function variation
- Thickness variations
8. Mathematical Frontiers
8.1 Machine Learning + Physics
- Physics-Informed Neural Networks (PINNs):
$$
\mathcal{L} = \mathcal{L}_{data} + \lambda \mathcal{L}_{physics}
$$
where $\mathcal{L}_{physics}$ enforces PDE residuals.
- Surrogate models for expensive TCAD simulations
- Inverse design and topology optimization
- Defect prediction in manufacturing
8.2 Stochastic Modeling
Random Dopant Fluctuation:
$$
\sigma_{V_T} \propto \frac{t_{ox}}{\sqrt{W \cdot L \cdot N_A}}
$$
Approaches:
- Atomistic Monte Carlo (place individual dopants)
- Statistical impedance field method
- Compact model statistical extensions
8.3 Multiphysics Coupling
Electro-Thermal Self-Heating:
$$
\rho C_p \frac{\partial T}{\partial t} =
abla \cdot (\kappa
abla T) + \mathbf{J} \cdot \mathbf{E}
$$
Stress Effects on Mobility (Piezoresistance):
$$
\frac{\Delta \mu}{\mu_0} = \pi_L \sigma_L + \pi_T \sigma_T
$$
Electromigration in Interconnects:
$$
\mathbf{J}_{atoms} = \frac{D C}{k_B T} \left( Z^* q \mathbf{E} - \Omega
abla \sigma \right)
$$
8.4 Atomistic-Continuum Bridging
Strategies:
- Coarse-graining from MD/DFT
- Density gradient quantum corrections:
$$
V_{QM} = \frac{\gamma \hbar^2}{12 m^*} \frac{
abla^2 \sqrt{n}}{\sqrt{n}}
$$
- Hybrid methods: atomistic core + continuum far-field
The mathematics of semiconductor manufacturing and device physics encompasses:
$$
\boxed{
\begin{aligned}
&\text{Process:} && \text{Stefan problems, diffusion PDEs, reaction kinetics} \\
&\text{Device:} && \text{Coupled Poisson + continuity equations} \\
&\text{Quantum:} && \text{Schrödinger, NEGF, tunneling} \\
&\text{Numerical:} && \text{FEM/FDM, Scharfetter-Gummel, Newton iteration} \\
&\text{Circuit:} && \text{Compact models (BSIM), variability statistics}
\end{aligned}
}
$$
Each level trades accuracy for computational tractability . The art lies in knowing when each approximation breaks down—and modern scaling is pushing us toward the quantum limit where classical continuum models become inadequate.