darkfield inspection,metrology
Scatter light off defects for enhanced detection.
58 technical terms and definitions
Scatter light off defects for enhanced detection.
Manufacturing date marking.
# Semiconductor Manufacturing: Dielectric Deposition Process (DDP) Modeling ## Overview **DDP (Dielectric Deposition Process)** refers to the set of techniques used to deposit insulating films in semiconductor fabrication. Dielectric materials serve critical functions: - **Gate dielectrics** — $\text{SiO}_2$, high-$\kappa$ materials like $\text{HfO}_2$ - **Interlayer dielectrics (ILD)** — isolating metal interconnect layers - **Spacer dielectrics** — defining transistor gate dimensions - **Passivation layers** — protecting finished devices - **Hard masks** — etch selectivity during patterning ## Dielectric Deposition Methods ### Primary Techniques | Method | Full Name | Temperature Range | Typical Applications | |--------|-----------|-------------------|---------------------| | **PECVD** | Plasma-Enhanced CVD | $200-400°C$ | $\text{SiO}_2$, $\text{SiN}_x$ for ILD, passivation | | **LPCVD** | Low-Pressure CVD | $400-800°C$ | High-quality $\text{Si}_3\text{N}_4$, poly-Si | | **HDPCVD** | High-Density Plasma CVD | $300-450°C$ | Gap-fill for trenches and vias | | **ALD** | Atomic Layer Deposition | $150-350°C$ | Ultra-thin gate dielectrics ($\text{HfO}_2$, $\text{Al}_2\text{O}_3$) | | **Thermal Oxidation** | — | $800-1200°C$ | Gate oxide ($\text{SiO}_2$) | | **Spin-on** | SOG/SOD | $100-400°C$ | Planarization layers | ### Selection Criteria - **Conformality requirements** — ALD > LPCVD > PECVD - **Thermal budget** — PECVD/ALD for low-$T$, thermal oxidation for high-quality - **Throughput** — CVD methods faster than ALD - **Film quality** — Thermal > LPCVD > PECVD generally ## Physics of Dielectric Deposition Modeling ### Fundamental Transport Equations Modeling dielectric deposition requires solving coupled partial differential equations for mass, momentum, and energy transport. #### Mass Transport (Species Concentration) $$ \frac{\partial C}{\partial t} + \nabla \cdot (\mathbf{v}C) = D\nabla^2 C + R $$ Where: - $C$ — species concentration $[\text{mol/m}^3]$ - $\mathbf{v}$ — velocity field $[\text{m/s}]$ - $D$ — diffusion coefficient $[\text{m}^2/\text{s}]$ - $R$ — reaction rate $[\text{mol/m}^3 \cdot \text{s}]$ #### Energy Balance $$ \rho C_p \left(\frac{\partial T}{\partial t} + \mathbf{v} \cdot \nabla T\right) = k\nabla^2 T + Q $$ Where: - $\rho$ — density $[\text{kg/m}^3]$ - $C_p$ — specific heat capacity $[\text{J/kg} \cdot \text{K}]$ - $k$ — thermal conductivity $[\text{W/m} \cdot \text{K}]$ - $Q$ — heat generation rate $[\text{W/m}^3]$ #### Momentum Balance (Navier-Stokes) $$ \rho\left(\frac{\partial \mathbf{v}}{\partial t} + \mathbf{v} \cdot \nabla \mathbf{v}\right) = -\nabla p + \mu \nabla^2 \mathbf{v} + \rho \mathbf{g} $$ Where: - $p$ — pressure $[\text{Pa}]$ - $\mu$ — dynamic viscosity $[\text{Pa} \cdot \text{s}]$ - $\mathbf{g}$ — gravitational acceleration $[\text{m/s}^2]$ ### Surface Reaction Kinetics #### Arrhenius Rate Expression $$ k = A \exp\left(-\frac{E_a}{RT}\right) $$ Where: - $k$ — rate constant - $A$ — pre-exponential factor - $E_a$ — activation energy $[\text{J/mol}]$ - $R$ — gas constant $= 8.314 \, \text{J/mol} \cdot \text{K}$ - $T$ — temperature $[\text{K}]$ #### Langmuir Adsorption Isotherm (for ALD) $$ \theta = \frac{K \cdot p}{1 + K \cdot p} $$ Where: - $\theta$ — fractional surface coverage $(0 \leq \theta \leq 1)$ - $K$ — equilibrium adsorption constant - $p$ — partial pressure of adsorbate #### Sticking Coefficient $$ S = S_0 \cdot (1 - \theta)^n \cdot \exp\left(-\frac{E_a}{RT}\right) $$ Where: - $S$ — sticking coefficient (probability of adsorption) - $S_0$ — initial sticking coefficient - $n$ — reaction order ### Plasma Modeling (PECVD/HDPCVD) #### Electron Energy Distribution Function (EEDF) For non-Maxwellian plasmas, the Druyvesteyn distribution: $$ f(\varepsilon) = C \cdot \varepsilon^{1/2} \exp\left(-\left(\frac{\varepsilon}{\bar{\varepsilon}}\right)^2\right) $$ Where: - $\varepsilon$ — electron energy $[\text{eV}]$ - $\bar{\varepsilon}$ — mean electron energy - $C$ — normalization constant #### Ion Bombardment Energy $$ E_{ion} = e \cdot V_{sheath} + \frac{1}{2}m_{ion}v_{Bohm}^2 $$ Where: - $V_{sheath}$ — plasma sheath voltage - $v_{Bohm} = \sqrt{\frac{k_B T_e}{m_{ion}}}$ — Bohm velocity #### Radical Generation Rate $$ R_{radical} = n_e \cdot n_{gas} \cdot \langle \sigma v \rangle $$ Where: - $n_e$ — electron density $[\text{m}^{-3}]$ - $n_{gas}$ — neutral gas density - $\langle \sigma v \rangle$ — rate coefficient (energy-averaged cross-section × velocity) ## Feature-Scale Modeling ### Critical Phenomena in High Aspect Ratio Structures Modern semiconductor devices require filling trenches and vias with aspect ratios (AR) exceeding 50:1. #### Knudsen Number $$ Kn = \frac{\lambda}{d} $$ Where: - $\lambda$ — mean free path of gas molecules - $d$ — characteristic feature dimension | Regime | Knudsen Number | Transport Type | |--------|---------------|----------------| | Continuum | $Kn < 0.01$ | Viscous flow | | Slip | $0.01 < Kn < 0.1$ | Transition | | Transition | $0.1 < Kn < 10$ | Mixed | | Free molecular | $Kn > 10$ | Ballistic/Knudsen | #### Mean Free Path Calculation $$ \lambda = \frac{k_B T}{\sqrt{2} \pi d_m^2 p} $$ Where: - $d_m$ — molecular diameter $[\text{m}]$ - $p$ — pressure $[\text{Pa}]$ ### Step Coverage Model $$ SC = \frac{t_{sidewall}}{t_{top}} \times 100\% $$ For diffusion-limited deposition: $$ SC \approx \frac{1}{\sqrt{1 + AR^2}} $$ For reaction-limited deposition: $$ SC \approx 1 - \frac{S \cdot AR}{2} $$ Where: - $S$ — sticking coefficient - $AR$ — aspect ratio = depth/width ### Void Formation Criterion Void formation occurs when: $$ \frac{d(thickness_{sidewall})}{dz} > \frac{w(z)}{2 \cdot t_{total}} $$ Where: - $w(z)$ — feature width at depth $z$ - $t_{total}$ — total deposition time ## Film Properties to Model ### Structural Properties - **Thickness uniformity**: $$ U = \frac{t_{max} - t_{min}}{t_{max} + t_{min}} \times 100\% $$ - **Film stress** (Stoney equation): $$ \sigma_f = \frac{E_s t_s^2}{6(1-\nu_s)t_f} \cdot \frac{1}{R} $$ Where: - $E_s$, $\nu_s$ — substrate Young's modulus and Poisson ratio - $t_s$, $t_f$ — substrate and film thickness - $R$ — radius of curvature - **Density from refractive index** (Lorentz-Lorenz): $$ \frac{n^2 - 1}{n^2 + 2} = \frac{4\pi}{3} N \alpha $$ Where $N$ is molecular density and $\alpha$ is polarizability ### Electrical Properties - **Dielectric constant** (capacitance method): $$ \kappa = \frac{C \cdot t}{\varepsilon_0 \cdot A} $$ - **Breakdown field**: $$ E_{BD} = \frac{V_{BD}}{t} $$ - **Leakage current density** (Fowler-Nordheim tunneling): $$ J = \frac{q^3 E^2}{8\pi h \phi_B} \exp\left(-\frac{8\pi\sqrt{2m^*}\phi_B^{3/2}}{3qhE}\right) $$ Where: - $E$ — electric field - $\phi_B$ — barrier height - $m^*$ — effective electron mass ## Multiscale Modeling Hierarchy ### Scale Linking Framework ``` ┌─────────────────────────────────────────────────────────────────────┐ │ ATOMISTIC (Å-nm) MESOSCALE (nm-μm) CONTINUUM │ │ ───────────────── ────────────────── (μm-mm) │ │ ────────── │ │ • DFT calculations • Kinetic Monte Carlo • CFD │ │ • Molecular Dynamics • Level-set methods • FEM │ │ • Ab initio MD • Cellular automata • TCAD │ │ │ │ Outputs: Outputs: Outputs: │ │ • Binding energies • Film morphology • Flow │ │ • Reaction barriers • Growth rate • T, C │ │ • Diffusion coefficients • Surface roughness • Profiles │ └─────────────────────────────────────────────────────────────────────┘ ``` ### DFT Calculations Solve the Kohn-Sham equations: $$ \left[-\frac{\hbar^2}{2m}\nabla^2 + V_{eff}(\mathbf{r})\right]\psi_i(\mathbf{r}) = \varepsilon_i \psi_i(\mathbf{r}) $$ Where: $$ V_{eff} = V_{ext} + V_H + V_{xc} $$ - $V_{ext}$ — external potential (nuclei) - $V_H$ — Hartree potential (electron-electron) - $V_{xc}$ — exchange-correlation potential ### Kinetic Monte Carlo (kMC) Event selection probability: $$ P_i = \frac{k_i}{\sum_j k_j} $$ Time advancement: $$ \Delta t = -\frac{\ln(r)}{\sum_j k_j} $$ Where $r$ is a random number $\in (0,1]$ ## Specific Process Examples ### PECVD $\text{SiO}_2$ from TEOS #### Overall Reaction $$ \text{Si(OC}_2\text{H}_5\text{)}_4 + 12\text{O}^* \xrightarrow{\text{plasma}} \text{SiO}_2 + 8\text{CO}_2 + 10\text{H}_2\text{O} $$ #### Key Process Parameters | Parameter | Typical Range | Effect | |-----------|--------------|--------| | RF Power | $100-1000 \, \text{W}$ | ↑ Power → ↑ Density, ↓ Dep rate | | Pressure | $0.5-5 \, \text{Torr}$ | ↑ Pressure → ↑ Dep rate, ↓ Conformality | | Temperature | $300-400°C$ | ↑ Temp → ↑ Density, ↓ H content | | TEOS:O₂ ratio | $1:5$ to $1:20$ | Affects stoichiometry, quality | #### Deposition Rate Model $$ R_{dep} = k_0 \cdot p_{TEOS}^a \cdot p_{O_2}^b \cdot \exp\left(-\frac{E_a}{RT}\right) $$ Typical values: $a \approx 0.5$, $b \approx 0.3$, $E_a \approx 0.3 \, \text{eV}$ ### ALD High-$\kappa$ Dielectrics ($\text{HfO}_2$) #### Half-Reactions **Cycle A (Metal precursor):** $$ \text{Hf(N(CH}_3\text{)}_2\text{)}_4\text{(g)} + \text{*-OH} \rightarrow \text{*-O-Hf(N(CH}_3\text{)}_2\text{)}_3 + \text{HN(CH}_3\text{)}_2 $$ **Cycle B (Oxidizer):** $$ \text{*-O-Hf(N(CH}_3\text{)}_2\text{)}_3 + 2\text{H}_2\text{O} \rightarrow \text{*-O-Hf(OH)}_3 + 3\text{HN(CH}_3\text{)}_2 $$ #### Growth Per Cycle (GPC) $$ \text{GPC} = \frac{\theta_{sat} \cdot \rho_{site} \cdot M_{HfO_2}}{\rho_{HfO_2} \cdot N_A} $$ Typical GPC for $\text{HfO}_2$: $0.8-1.2 \, \text{Å/cycle}$ #### ALD Window ``` ┌────────────────────────────┐ GPC │ ┌──────────────┐ │ (Å/ │ /│ │\ │ cycle) │ / │ ALD │ \ │ │ / │ WINDOW │ \ │ │ / │ │ \ │ │/ │ │ \ │ └─────┴──────────────┴─────┴─┘ T_min T_max Temperature (°C) ``` Below $T_{min}$: Condensation, incomplete reactions Above $T_{max}$: Precursor decomposition, CVD-like behavior ### HDPCVD Gap Fill #### Deposition-Etch Competition Net deposition rate: $$ R_{net}(z) = R_{dep}(\theta) - R_{etch}(E_{ion}, \theta) $$ Where: - $R_{dep}(\theta)$ — angular-dependent deposition rate - $R_{etch}$ — ion-enhanced etch rate - $\theta$ — angle from surface normal #### Sputter Yield (Yamamura Formula) $$ Y(E, \theta) = Y_0(E) \cdot f(\theta) $$ Where: $$ f(\theta) = \cos^{-f}\theta \cdot \exp\left[-\Sigma(\cos^{-1}\theta - 1)\right] $$ ## Machine Learning Applications ### Virtual Metrology **Objective:** Predict film properties from in-situ sensor data without destructive measurement. $$ \hat{y} = f_{ML}(\mathbf{x}_{sensors}, \mathbf{x}_{recipe}) $$ Where: - $\hat{y}$ — predicted property (thickness, stress, etc.) - $\mathbf{x}_{sensors}$ — OES, pressure, RF power signals - $\mathbf{x}_{recipe}$ — setpoints and timing ### Gaussian Process Regression $$ y(\mathbf{x}) \sim \mathcal{GP}\left(m(\mathbf{x}), k(\mathbf{x}, \mathbf{x}')\right) $$ Posterior mean prediction: $$ \mu(\mathbf{x}^*) = \mathbf{k}^T(\mathbf{K} + \sigma_n^2\mathbf{I})^{-1}\mathbf{y} $$ Uncertainty quantification: $$ \sigma^2(\mathbf{x}^*) = k(\mathbf{x}^*, \mathbf{x}^*) - \mathbf{k}^T(\mathbf{K} + \sigma_n^2\mathbf{I})^{-1}\mathbf{k} $$ ### Bayesian Optimization for Recipe Development **Acquisition function** (Expected Improvement): $$ \text{EI}(\mathbf{x}) = \mathbb{E}\left[\max(f(\mathbf{x}) - f^+, 0)\right] $$ Where $f^+$ is the best observed value. ## Advanced Node Challenges (Sub-5nm) ### Critical Challenges | Challenge | Technical Details | Modeling Complexity | |-----------|------------------|---------------------| | **Ultra-high AR** | 3D NAND: 100+ layers, AR > 50:1 | Knudsen transport, ballistic modeling | | **Atomic precision** | Gate dielectrics: 1-2 nm | Monolayer-level control, quantum effects | | **Low-$\kappa$ integration** | $\kappa < 2.5$ porous films | Mechanical integrity, plasma damage | | **Selective deposition** | Area-selective ALD | Nucleation control, surface chemistry | | **Thermal budget** | BEOL: $< 400°C$ | Kinetic limitations, precursor chemistry | ### Equivalent Oxide Thickness (EOT) For high-$\kappa$ gate stacks: $$ \text{EOT} = t_{IL} + \frac{\kappa_{SiO_2}}{\kappa_{high-k}} \cdot t_{high-k} $$ Where: - $t_{IL}$ — interfacial layer thickness - $\kappa_{SiO_2} = 3.9$ - Typical high-$\kappa$: $\kappa_{HfO_2} \approx 20-25$ ### Low-$\kappa$ Dielectric Design Effective dielectric constant: $$ \kappa_{eff} = \kappa_{matrix} \cdot (1 - p) + \kappa_{air} \cdot p $$ Where $p$ is porosity fraction. Target for advanced nodes: $\kappa_{eff} < 2.0$ ## Tools and Software ### Commercial TCAD - **Synopsys Sentaurus Process** — full process simulation - **Silvaco Victory Process** — alternative TCAD suite - **Lam Research SEMulator3D** — 3D topography simulation ### Multiphysics Platforms - **COMSOL Multiphysics** — coupled PDE solving - **Ansys Fluent** — CFD for reactor design - **Ansys CFX** — alternative CFD solver ### Specialized Tools - **CHEMKIN** (Ansys) — gas-phase reaction kinetics - **Reaction Design** — combustion and plasma chemistry - **Custom Monte Carlo codes** — feature-scale simulation ### Open Source Options - **OpenFOAM** — CFD framework - **LAMMPS** — molecular dynamics - **Quantum ESPRESSO** — DFT calculations - **SPARTA** — DSMC for rarefied gas dynamics ## Summary Dielectric deposition modeling in semiconductor manufacturing integrates: 1. **Transport phenomena** — mass, momentum, energy conservation 2. **Reaction kinetics** — surface and gas-phase chemistry 3. **Plasma physics** — for PECVD/HDPCVD processes 4. **Feature-scale physics** — conformality, void formation 5. **Multiscale approaches** — atomistic to continuum 6. **Machine learning** — for optimization and virtual metrology The goal is predicting and optimizing film properties based on process parameters while accounting for the extreme topography of modern semiconductor devices.
Separate temporarily bonded wafers.
Measure trap levels.
High-aspect-ratio etching for TSV.
Spatial distribution of defects.
Automated optical or e-beam inspection to find particles scratches defects.
High-resolution follow-up of detected defects.
Trace defects to source.
Remove molding flash.
Film thickness grown per unit time.
Model film growth.
Range of wafer height where image stays in focus.
Acceptable focus range.
Moisture absorber.
Systematically vary parameters to optimize.
Lowest detectable amount.
Chemical process to remove soluble resist revealing pattern.
# 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} = \nabla \cdot \left( D(C, T) \nabla 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[ -\nu \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 |\nabla \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 \nabla^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) $$ \nabla \cdot (\varepsilon \nabla \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} \nabla \cdot \mathbf{J}_n + G - R $$ Holes: $$ \frac{\partial p}{\partial t} = -\frac{1}{q} \nabla \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 \nabla n $$ $$ \mathbf{J}_p = q \mu_p p \mathbf{E} - q D_p \nabla 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: ```text ┌─────────────────────────────────────────────────────────────┐ │ ┌───────────────┐ ┌───────────────┐ ┌───────────────┐ │ │ │ Process │──▶│ Device │──▶│ Parameter │ │ │ │ Simulation │ │ Simulation │ │ Extraction │ │ │ │ (Sentaurus) │ │ (Sentaurus) │ │ (BSIM Fit) │ │ │ └───────────────┘ └───────────────┘ └───────────────┘ │ │ │ │ │ │ │ ▼ ▼ ▼ │ │ ┌───────────────┐ ┌───────────────┐ ┌───────────────┐ │ │ │• Implantation │ │• I-V, C-V │ │• BSIM params │ │ │ │• Diffusion │ │• Breakdown │ │• Corner extr. │ │ │ │• Oxidation │ │• Hot carrier │ │• Variability │ │ │ │• Etching │ │• Noise │ │ statistics │ │ │ └───────────────┘ └───────────────┘ └───────────────┘ │ │ │ │ │ ▼ │ │ ┌───────────────┐ │ │ │ Circuit │ │ │ │ Simulation │ │ │ │(SPICE,Spectre)│ │ │ └───────────────┘ │ └─────────────────────────────────────────────────────────────┘ ``` 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} = \nabla \cdot (\kappa \nabla 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 \nabla \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{\nabla^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.
# Device Physics, TCAD, and Mathematical Modeling 1. Physical Foundation 1.1 Band Theory and Electronic Structure - Energy bands arise from the periodic potential of the crystal lattice - Conduction band (empty states available for electron transport) - Valence band (filled states; holes represent missing electrons) - Bandgap $E_g$ separates these bands (Si: ~1.12 eV at 300K) - Effective mass approximation - Electrons and holes behave as quasi-particles with modified mass - Electron effective mass: $m_n^*$ - Hole effective mass: $m_p^*$ - Carrier statistics follow Fermi-Dirac distribution: $$ f(E) = \frac{1}{1 + \exp\left(\frac{E - E_F}{k_B T}\right)} $$ - Carrier concentrations in non-degenerate semiconductors: $$ n = N_C \exp\left(-\frac{E_C - E_F}{k_B T}\right) $$ $$ p = N_V \exp\left(-\frac{E_F - E_V}{k_B T}\right) $$ Where: - $N_C$, $N_V$ = effective density of states in conduction/valence bands - $E_C$, $E_V$ = conduction/valence band edges - $E_F$ = Fermi level 1.2 Carrier Transport Mechanisms | Mechanism | Driving Force | Current Density | |-----------|---------------|-----------------| | Drift | Electric field $\mathbf{E}$ | $\mathbf{J} = qn\mu\mathbf{E}$ | | Diffusion | Concentration gradient | $\mathbf{J} = qD\nabla n$ | | Thermionic emission | Thermal energy over barrier | Exponential in $\phi_B/k_BT$ | | Tunneling | Quantum penetration | Exponential in barrier | - Einstein relation connects mobility and diffusivity: $$ D = \frac{k_B T}{q} \mu $$ 1.3 Generation and Recombination - Thermal equilibrium condition: $$ np = n_i^2 $$ - Three primary recombination mechanisms: 1. Shockley-Read-Hall (SRH) — trap-assisted 2. Auger — three-particle process (dominant at high injection) 3. Radiative — photon emission (important in direct bandgap materials) 2. Mathematical Hierarchy 2.1 Quantum Mechanical Level (Most Fundamental) Time-Independent Schrödinger Equation $$ \left[-\frac{\hbar^2}{2m^*}\nabla^2 + V(\mathbf{r})\right]\psi = E\psi $$ Where: - $\hbar$ = reduced Planck constant - $m^*$ = effective mass - $V(\mathbf{r})$ = potential energy - $\psi$ = wavefunction - $E$ = energy eigenvalue Non-Equilibrium Green's Function (NEGF) For open quantum systems (nanoscale devices, tunneling): $$ G^R = [EI - H - \Sigma]^{-1} $$ - $G^R$ = retarded Green's function - $H$ = device Hamiltonian - $\Sigma$ = self-energy (encodes contact coupling) Applications: - Tunnel FETs - Ultra-scaled MOSFETs ($L_g < 10$ nm) - Quantum well devices - Resonant tunneling diodes 2.2 Boltzmann Transport Level Boltzmann Transport Equation (BTE) $$ \frac{\partial f}{\partial t} + \mathbf{v} \cdot \nabla_{\mathbf{r}} f + \frac{\mathbf{F}}{\hbar} \cdot \nabla_{\mathbf{k}} f = \left(\frac{\partial f}{\partial t}\right)_{\text{coll}} $$ Where: - $f(\mathbf{r}, \mathbf{k}, t)$ = distribution function in phase space - $\mathbf{v}$ = group velocity - $\mathbf{F}$ = external force - RHS = collision integral Solution Methods: - Monte Carlo (stochastic particle tracking) - Spherical Harmonics Expansion (SHE) - Moments methods → leads to drift-diffusion, hydrodynamic Captures: - Hot carrier effects - Velocity overshoot - Non-equilibrium distributions - Ballistic transport 2.3 Hydrodynamic / Energy Balance Level Derived from moments of BTE with carrier temperature as variable: $$ \frac{\partial (nw)}{\partial t} + \nabla \cdot \mathbf{S} = \mathbf{J} \cdot \mathbf{E} - \frac{n(w - w_0)}{\tau_w} $$ - $w$ = carrier energy density - $\mathbf{S}$ = energy flux - $\tau_w$ = energy relaxation time - $w_0$ = equilibrium energy density Key feature: Carrier temperature $T_n \neq$ lattice temperature $T_L$ 2.4 Drift-Diffusion Level (The Workhorse) The most widely used TCAD formulation — three coupled PDEs: Poisson's Equation (Electrostatics) $$ \nabla \cdot (\varepsilon \nabla \psi) = -\rho = -q(p - n + N_D^+ - N_A^-) $$ - $\psi$ = electrostatic potential - $\varepsilon$ = permittivity - $\rho$ = charge density - $N_D^+$, $N_A^-$ = ionized donor/acceptor concentrations Electron Continuity Equation $$ \frac{\partial n}{\partial t} = \frac{1}{q}\nabla \cdot \mathbf{J}_n + G_n - R_n $$ Hole Continuity Equation $$ \frac{\partial p}{\partial t} = -\frac{1}{q}\nabla \cdot \mathbf{J}_p + G_p - R_p $$ Current Density Equations Standard form: $$ \mathbf{J}_n = q\mu_n n \mathbf{E} + qD_n \nabla n $$ $$ \mathbf{J}_p = q\mu_p p \mathbf{E} - qD_p \nabla p $$ Quasi-Fermi level formulation: $$ \mathbf{J}_n = q\mu_n n \nabla E_{F,n} $$ $$ \mathbf{J}_p = q\mu_p p \nabla E_{F,p} $$ System characteristics: - Coupled, nonlinear, elliptic-parabolic PDEs - Carrier concentrations vary exponentially with potential - Spans 10+ orders of magnitude across junctions 3. Numerical Methods 3.1 Spatial Discretization Finite Difference Method (FDM) - Simple implementation - Limited to structured (rectangular) grids - Box integration for conservation Finite Element Method (FEM) - Handles complex geometries - Basis function expansion - Weak (variational) formulation Finite Volume Method (FVM) - Ensures local conservation - Natural for semiconductor equations - Control volume integration 3.2 Scharfetter-Gummel Discretization Critical for numerical stability — handles exponential carrier variations: $$ J_{n,i+\frac{1}{2}} = \frac{qD_n}{h}\left[n_i B\left(\frac{\psi_i - \psi_{i+1}}{V_T}\right) - n_{i+1} B\left(\frac{\psi_{i+1} - \psi_i}{V_T}\right)\right] $$ Where the Bernoulli function is: $$ B(x) = \frac{x}{e^x - 1} $$ Properties: - Reduces to central difference for small $\Delta\psi$ - Reduces to upwind for large $\Delta\psi$ - Prevents spurious oscillations - Thermal voltage: $V_T = k_B T / q \approx 26$ mV at 300K 3.3 Mesh Generation - 2D: Delaunay triangulation - 3D: Tetrahedral meshing Adaptive refinement criteria: - Junction regions (high field gradients) - Oxide interfaces - Contact regions - High current density areas Quality metrics: - Aspect ratio - Orthogonality (important for FVM) - Delaunay property (circumsphere criterion) 3.4 Nonlinear Solvers Gummel Iteration (Decoupled) repeat: 1. Solve Poisson equation → ψ 2. Solve electron continuity → n 3. Solve hole continuity → p until convergence Pros: - Simple implementation - Robust for moderate bias - Each subproblem is smaller Cons: - Poor convergence at high injection - Slow for strongly coupled systems Newton-Raphson (Fully Coupled) Solve the linearized system: $$ \mathbf{J} \cdot \delta\mathbf{x} = -\mathbf{F}(\mathbf{x}) $$ Where: - $\mathbf{J}$ = Jacobian matrix $\partial \mathbf{F}/\partial \mathbf{x}$ - $\mathbf{F}$ = residual vector - $\delta\mathbf{x}$ = update vector Pros: - Quadratic convergence near solution - Handles strong coupling Cons: - Requires good initial guess - Expensive Jacobian assembly - Larger linear systems Hybrid Methods - Start with Gummel to get close - Switch to Newton for fast final convergence 3.5 Linear Solvers For large, sparse, ill-conditioned Jacobian systems: | Method | Type | Characteristics | |--------|------|-----------------| | LU (PARDISO, UMFPACK) | Direct | Robust, memory-intensive | | GMRES | Iterative | Krylov subspace, needs preconditioning | | BiCGSTAB | Iterative | Non-symmetric systems | | Multigrid | Iterative | Optimal for Poisson-like equations | 4. Physical Models in TCAD 4.1 Mobility Models Matthiessen's Rule Combines independent scattering mechanisms: $$ \frac{1}{\mu} = \frac{1}{\mu_{\text{lattice}}} + \frac{1}{\mu_{\text{impurity}}} + \frac{1}{\mu_{\text{surface}}} + \cdots $$ Lattice Scattering $$ \mu_L = \mu_0 \left(\frac{T}{300}\right)^{-\alpha} $$ - Si electrons: $\alpha \approx 2.4$ - Si holes: $\alpha \approx 2.2$ Ionized Impurity Scattering Brooks-Herring model: $$ \mu_I \propto \frac{T^{3/2}}{N_I \cdot \ln(1 + b^2) - b^2/(1+b^2)} $$ High-Field Saturation (Caughey-Thomas) $$ \mu(E) = \frac{\mu_0}{\left[1 + \left(\frac{\mu_0 E}{v_{\text{sat}}}\right)^\beta\right]^{1/\beta}} $$ - $v_{\text{sat}}$ = saturation velocity (~$10^7$ cm/s for Si) - $\beta$ = fitting parameter (~2 for electrons, ~1 for holes) 4.2 Recombination Models Shockley-Read-Hall (SRH) $$ R_{\text{SRH}} = \frac{np - n_i^2}{\tau_p(n + n_1) + \tau_n(p + p_1)} $$ Where: - $\tau_n$, $\tau_p$ = carrier lifetimes - $n_1 = n_i \exp[(E_t - E_i)/k_BT]$ - $p_1 = n_i \exp[(E_i - E_t)/k_BT]$ - $E_t$ = trap energy level Auger Recombination $$ R_{\text{Auger}} = (C_n n + C_p p)(np - n_i^2) $$ - $C_n$, $C_p$ = Auger coefficients (~$10^{-31}$ cm$^6$/s for Si) - Dominant at high carrier densities ($>10^{18}$ cm$^{-3}$) Radiative Recombination $$ R_{\text{rad}} = B(np - n_i^2) $$ - $B$ = radiative coefficient - Important in direct bandgap materials (GaAs, InP) 4.3 Band-to-Band Tunneling For tunnel FETs, Zener diodes: $$ G_{\text{BTBT}} = A \cdot E^2 \exp\left(-\frac{B}{E}\right) $$ - $A$, $B$ = material-dependent parameters - $E$ = electric field magnitude 4.4 Quantum Corrections Density Gradient Method Adds quantum potential to classical equations: $$ V_Q = -\frac{\hbar^2}{6m^*} \frac{\nabla^2\sqrt{n}}{\sqrt{n}} $$ Or equivalently, the quantum potential term: $$ \Lambda_n = \frac{\hbar^2}{12 m_n^* k_B T} \nabla^2 \ln(n) $$ Applications: - Inversion layer quantization in MOSFETs - Thin body SOI devices - FinFETs, nanowires 1D Schrödinger-Poisson For stronger quantum confinement: 1. Solve 1D Schrödinger in confinement direction → subbands $E_i$, $\psi_i$ 2. Calculate 2D density of states 3. Compute carrier density from subband occupation 4. Solve 2D Poisson with quantum charge 5. Iterate to self-consistency 4.5 Bandgap Narrowing At high doping ($N > 10^{17}$ cm$^{-3}$): $$ \Delta E_g = A \cdot N^{1/3} + B \cdot \ln\left(\frac{N}{N_{\text{ref}}}\right) $$ Effect: Increases $n_i^2$ → affects recombination and device characteristics 4.6 Interface Models - Interface trap density: $D_{it}(E)$ — states per cm$^2$·eV - Oxide charges: - Fixed oxide charge $Q_f$ - Mobile ionic charge $Q_m$ - Oxide trapped charge $Q_{ot}$ - Interface trapped charge $Q_{it}$ 5. Process TCAD 5.1 Ion Implantation Monte Carlo Method - Track individual ion trajectories - Binary collision approximation - Accurate for low doses, complex geometries Analytical Profiles Gaussian: $$ N(x) = \frac{\Phi}{\sqrt{2\pi}\Delta R_p} \exp\left[-\frac{(x - R_p)^2}{2\Delta R_p^2}\right] $$ - $\Phi$ = dose (ions/cm$^2$) - $R_p$ = projected range - $\Delta R_p$ = straggle Pearson IV: Adds skewness and kurtosis for better accuracy 5.2 Diffusion Fick's First Law: $$ \mathbf{J} = -D \nabla C $$ Fick's Second Law: $$ \frac{\partial C}{\partial t} = \nabla \cdot (D \nabla C) $$ Concentration-dependent diffusion: $$ D = D_i \left(\frac{n}{n_i}\right)^2 + D_v + D_x \left(\frac{n}{n_i}\right) $$ (Accounts for charged point defects) 5.3 Oxidation Deal-Grove Model: $$ x_{ox}^2 + A \cdot x_{ox} = B(t + \tau) $$ - $x_{ox}$ = oxide thickness - $A$, $B$ = temperature-dependent parameters - Linear regime: $x_{ox} \approx (B/A) \cdot t$ (thin oxide) - Parabolic regime: $x_{ox} \approx \sqrt{B \cdot t}$ (thick oxide) 5.4 Etching and Deposition Level-set method for surface evolution: $$ \frac{\partial \phi}{\partial t} + v_n |\nabla \phi| = 0 $$ - $\phi$ = level-set function (zero contour = surface) - $v_n$ = normal velocity (etch/deposition rate) 6. Multiphysics and Advanced Topics 6.1 Electrothermal Coupling Heat equation: $$ \rho c_p \frac{\partial T}{\partial t} = \nabla \cdot (\kappa \nabla T) + H $$ Heat generation: $$ H = \mathbf{J} \cdot \mathbf{E} + (R - G)(E_g + 3k_BT) $$ - First term: Joule heating - Second term: recombination heating Thermoelectric effects: - Seebeck effect - Peltier effect - Thomson effect 6.2 Electromechanical Coupling Strain effects on mobility: $$ \mu_{\text{strained}} = \mu_0 (1 + \Pi \cdot \sigma) $$ - $\Pi$ = piezoresistance coefficient - $\sigma$ = mechanical stress Applications: Strained Si, SiGe channels 6.3 Statistical Variability Sources of random variation: - Random Dopant Fluctuations (RDF) — discrete dopant positions - Line Edge Roughness (LER) — gate patterning variation - Metal Gate Granularity (MGG) — work function variation - Oxide Thickness Variation (OTV) Simulation approach: - Monte Carlo sampling over device instances - Statistical TCAD → threshold voltage distributions 6.4 Reliability Modeling Bias Temperature Instability (BTI): - Defect generation at Si/SiO$_2$ interface - Reaction-diffusion models Hot Carrier Injection (HCI): - High-energy carriers damage interface - Coupled with energy transport 6.5 Noise Modeling Noise sources: - Thermal noise: $S_I = 4k_BT/R$ - Shot noise: $S_I = 2qI$ - 1/f noise (flicker): $S_I \propto I^2/(f \cdot N)$ Impedance field method for spatial correlation 7. Computational Architecture 7.1 Model Hierarchy Comparison | Level | Physics | Math | Cost | Accuracy | |-------|---------|------|------|----------| | NEGF | Quantum coherence | $G = [E-H-\Sigma]^{-1}$ | $$$$$ | Highest | | Monte Carlo | Distribution function | Stochastic DEs | $$$$ | High | | Hydrodynamic | Carrier temperature | Hyperbolic-parabolic PDEs | $$$ | Good | | Drift-Diffusion | Continuum transport | Elliptic-parabolic PDEs | $$ | Moderate | | Compact Models | Empirical | Algebraic | $ | Calibrated | 7.2 Software Architecture ```text ┌─────────────────────────────────────────┐ │ User Interface (GUI) │ ├─────────────────────────────────────────┤ │ Structure Definition │ │ (Geometry, Mesh, Materials) │ ├─────────────────────────────────────────┤ │ Physical Models │ │ (Mobility, Recombination, Quantum) │ ├─────────────────────────────────────────┤ │ Numerical Engine │ │ (Discretization, Solvers, Linear Alg) │ ├─────────────────────────────────────────┤ │ Post-Processing │ │ (Visualization, Parameter Extraction) │ └─────────────────────────────────────────┘ ``` 7.3 TCAD ↔ Compact Model Flow ```text ┌──────────┐ calibrate ┌──────────────┐ │ TCAD │ ──────────────► │ Compact Model│ │(Physics) │ │ (BSIM,PSP) │ └──────────┘ └──────────────┘ │ │ │ validate │ enable ▼ ▼ ┌──────────┐ ┌──────────────┐ │ Silicon │ │ Circuit │ │ Data │ │ Simulation │ └──────────┘ └──────────────┘ ``` Equations: Fundamental Constants | Symbol | Name | Value | |--------|------|-------| | $q$ | Elementary charge | $1.602 \times 10^{-19}$ C | | $k_B$ | Boltzmann constant | $1.381 \times 10^{-23}$ J/K | | $\hbar$ | Reduced Planck | $1.055 \times 10^{-34}$ J·s | | $\varepsilon_0$ | Vacuum permittivity | $8.854 \times 10^{-12}$ F/m | | $V_T$ | Thermal voltage (300K) | 25.9 mV | Silicon Properties (300K) | Property | Value | |----------|-------| | Bandgap $E_g$ | 1.12 eV | | Intrinsic carrier density $n_i$ | $1.0 \times 10^{10}$ cm$^{-3}$ | | Electron mobility $\mu_n$ | 1450 cm$^2$/V·s | | Hole mobility $\mu_p$ | 500 cm$^2$/V·s | | Electron saturation velocity | $1.0 \times 10^7$ cm/s | | Relative permittivity $\varepsilon_r$ | 11.7 |
Wafer containing actual devices.
Mechanical displacement gauge.
Adhesive at die edge.
Adhesives for die attachment.
Bond line thickness.
Voids in die attach layer.
Attach die to substrate or package.
Attach individual die to substrate or package.
Die breaking during bonding.
Number of chips that fit on one wafer.
Die per wafer quantifies how many complete chips fit on wafer determining manufacturing efficiency.
Lateral die displacement.
Non-parallel die placement.
Connections between chiplets.
# Semiconductor Die Dicing and Singulation: Mathematical Modeling Framework ## 1. Introduction and Overview Die dicing (singulation) is the process of separating individual semiconductor dies from a processed wafer. The mathematical modeling of this process spans multiple physics domains: - **Mechanics** - Cutting forces, stress distribution - **Thermodynamics** - Heat generation, thermal diffusion - **Fracture Mechanics** - Crack propagation, chipping - **Plasma Chemistry** - Etch rates, selectivity ### 1.1 Primary Dicing Methods | Method | Mechanism | Key Physics | |--------|-----------|-------------| | Blade Dicing | Mechanical abrasion | Contact mechanics, grinding | | Laser Dicing | Thermal ablation | Heat transfer, phase change | | Stealth Dicing | Internal modification | Nonlinear optics, fracture | | Plasma Dicing | Chemical etching | Plasma physics, reaction kinetics | ## 2. Wafer Geometry and Die Yield Calculations ### 2.1 Dies Per Wafer (DPW) - Basic Formula The fundamental yield calculation determines how many dies can be extracted from a circular wafer: $$ \text{DPW} = \frac{\pi \cdot \left(\frac{D_{wafer}}{2}\right)^2}{A_{die} + A_{scribe}} \times U $$ **Where:** - $D_{wafer}$ = Wafer diameter (mm) - Standard sizes: 150mm, 200mm, 300mm - $A_{die}$ = Die area (mm²) - $A_{die} = L_{die} \times W_{die}$ - $A_{scribe}$ = Scribe lane area per die (mm²) - $U$ = Utilization factor - Typical range: $0.75 \leq U \leq 0.95$ ### 2.2 Edge Loss Correction (Gross Die Formula) Accounts for partial dies at wafer edge: $$ N_{gross} = \frac{\pi \cdot r^2}{A_{die}} - \frac{\pi \cdot D_{wafer}}{\sqrt{2} \cdot S_{die}} $$ **Where:** - $r$ = Wafer radius (mm) - $S_{die}$ = Die size (characteristic dimension) ### 2.3 Rectangular Die Packing For rectangular dies on circular wafer: $$ N_{dies} = \left\lfloor \frac{D_{eff}}{L + k_x} \right\rfloor \times \left\lfloor \frac{D_{eff}}{W + k_y} \right\rfloor $$ **Where:** - $D_{eff}$ = Effective wafer diameter after edge exclusion - $L, W$ = Die length and width - $k_x, k_y$ = Kerf widths in x and y directions - $\lfloor \cdot \rfloor$ = Floor function ### 2.4 Die Yield Models #### 2.4.1 Poisson Yield Model (Simple) $$ Y = e^{-D_0 \cdot A_{die}} $$ **Where:** - $Y$ = Yield (fraction of good dies) - $D_0$ = Defect density (defects/cm²) - $A_{die}$ = Die area (cm²) #### 2.4.2 Murphy's Yield Model $$ Y = \left( \frac{1 - e^{-D_0 \cdot A_{die}}}{D_0 \cdot A_{die}} \right)^2 $$ #### 2.4.3 Seeds Yield Model $$ Y = e^{-\sqrt{D_0 \cdot A_{die}}} $$ #### 2.4.4 Comprehensive Die Yield Formula $$ Y = \left( \frac{W_a}{S + 2\sqrt{A \cdot D}} \right)^2 \cdot e^{-\pi \cdot D \cdot \frac{D}{A}} $$ **Where:** - $W_a$ = Wafer area - $S$ = Chip size - $A$ = Chip area - $D$ = Defect density ## 3. Mechanical Blade Dicing Models ### 3.1 Street Width Capability $$ W_{street} = P + K + n \cdot \sigma $$ **Where:** - $W_{street}$ = Minimum street width required - $P$ = Cut placement accuracy factor - $K$ = Kerf width factor - $n$ = Process capability factor (typically 3 for 3σ) - $\sigma$ = Standard deviation of cut position ### 3.2 Kerf Width Calculation $$ K_{eff} = K_{blade} + 2 \cdot \delta_{wear} + \delta_{runout} $$ **Where:** - $K_{eff}$ = Effective kerf width - $K_{blade}$ = Nominal blade thickness - $\delta_{wear}$ = Blade wear contribution - $\delta_{runout}$ = Spindle runout contribution ### 3.3 Maximum Undeformed Chip Thickness For grinding/dicing processes: $$ h_{max} = \left( \frac{4 \cdot v_f}{v_s \cdot C \cdot r} \right)^{1/2} \cdot \left( \frac{a_p}{d_s} \right)^{1/4} $$ **Where:** - $h_{max}$ = Maximum undeformed chip thickness - $v_f$ = Feed rate (mm/s) - $v_s$ = Blade surface velocity (m/s) - $C$ = Active cutting point density - $r$ = Chip width-to-thickness ratio - $a_p$ = Depth of cut - $d_s$ = Blade diameter **Optimal range for damage-free cutting:** $$ 10 \text{ nm} \leq h_{max} \leq 30 \text{ nm} $$ ### 3.4 Cutting Force Model #### 3.4.1 Total Cutting Force $$ \vec{F}_{total} = \vec{F}_t + \vec{F}_n + \vec{F}_f $$ **Where:** - $\vec{F}_t$ = Tangential cutting force - $\vec{F}_n$ = Normal force - $\vec{F}_f$ = Friction force #### 3.4.2 Tangential Force Component $$ F_t = K_s \cdot A_c = K_s \cdot a_p \cdot f $$ **Where:** - $K_s$ = Specific cutting force (N/mm²) - $A_c$ = Chip cross-sectional area - $a_p$ = Depth of cut - $f$ = Feed per revolution #### 3.4.3 Spindle Torque $$ T = F_t \cdot \frac{d_s}{2} $$ ### 3.5 Material Removal Rate (MRR) $$ \text{MRR} = v_f \cdot a_p \cdot K_{eff} $$ **Units:** mm³/s ### 3.6 Blade Wear Model $$ \frac{dV_w}{dL} = K_w \cdot \frac{F_n}{H_b} $$ **Where:** - $V_w$ = Blade wear volume - $L$ = Cutting length - $K_w$ = Wear coefficient - $F_n$ = Normal force - $H_b$ = Blade hardness ## 4. Chipping and Fracture Mechanics Models ### 4.1 Linear Elastic Fracture Mechanics (LEFM) #### 4.1.1 Stress Intensity Factor (Mode I) $$ K_I = \sigma \sqrt{\pi a} \cdot f\left(\frac{a}{W}\right) $$ **Where:** - $K_I$ = Mode I stress intensity factor (MPa$\cdot$√m) - $\sigma$ = Applied stress - $a$ = Crack length - $W$ = Specimen width - $f(a/W)$ = Geometry correction factor #### 4.1.2 Fracture Criterion Fracture occurs when: $$ K_I \geq K_{IC} $$ **Silicon fracture toughness values:** | Plane | $K_{IC}$ (MPa$\cdot$√m) | |-------|-------------------| | {111} | 0.82 | | {110} | 0.90 | | {100} | 0.95 | ### 4.2 Critical Crack Length for Backside Chipping Based on single edge cracked pure bending: $$ K_C = \frac{6M}{bh^2} \sqrt{\pi a_c} \cdot Y(\alpha) $$ **Where:** - $\alpha = \frac{a_c}{h}$ - For $a_c \ll h$: $Y(\alpha) \approx 1.122$ **Solving for critical crack length:** $$ a_c = \frac{1}{\pi} \left( \frac{K_C \cdot b}{6M_x \cdot Y} \right)^2 \cdot h^4 $$ **With uniform loading:** $$ M_x = \frac{q \cdot l^2}{2} $$ ### 4.3 Chipping Mode Classification Based on crystallographic orientation of silicon: | Mode | Angle | Mechanism | |------|-------|-----------| | Type 1 | 30° | Radial crack along ⟨110⟩ meets {111} cleavage | | Type 2 | 60° | Radial crack intersection | | Type 3 | 90° | Perpendicular radial crack | | Type 4 | Irregular | Multiple crack interactions | ### 4.4 Chipping Size Prediction $$ C_{size} = k \cdot \left( \frac{F_n}{H} \right)^{2/3} \cdot \left( \frac{E}{K_{IC}} \right)^{1/2} $$ **Where:** - $C_{size}$ = Chipping size - $F_n$ = Normal force - $H$ = Material hardness - $E$ = Young's modulus - $K_{IC}$ = Fracture toughness ### 4.5 Die Strength (Three-Point Bending) $$ \sigma_{break} = \frac{3 \cdot F \cdot L}{2 \cdot b \cdot t^2} $$ **Where:** - $\sigma_{break}$ = Breaking strength (MPa) - $F$ = Applied force at fracture - $L$ = Support span - $b$ = Die width - $t$ = Die thickness ### 4.6 Weibull Distribution for Die Strength $$ P_f(\sigma) = 1 - \exp\left[ -\left( \frac{\sigma}{\sigma_0} \right)^m \right] $$ **Where:** - $P_f$ = Probability of failure - $\sigma_0$ = Characteristic strength - $m$ = Weibull modulus (shape parameter) ## 5. Laser Dicing Mathematical Models ### 5.1 Laser Beam Fundamentals #### 5.1.1 Gaussian Beam Intensity Profile $$ I(r) = I_0 \cdot \exp\left( -\frac{2r^2}{w^2} \right) $$ **Where:** - $I_0$ = Peak intensity at beam center - $r$ = Radial distance from center - $w$ = Beam radius (1/e² point) #### 5.1.2 Peak Fluence $$ F_0 = \frac{2 \cdot E_p}{\pi \cdot w_0^2} $$ **Where:** - $F_0$ = Peak fluence (J/cm²) - $E_p$ = Pulse energy (J) - $w_0$ = Beam waist radius ### 5.2 Thermal Diffusion Model #### 5.2.1 Thermal Diffusion Length $$ L_{th} = \sqrt{D_{th} \cdot t_p} $$ **Where:** - $L_{th}$ = Thermal diffusion length - $D_{th}$ = Thermal diffusivity - Silicon: $D_{th} \approx 0.8$ cm²/s - $t_p$ = Pulse duration #### 5.2.2 Heat Penetration Depth $$ \delta_{th} = 2\sqrt{D_{th} \cdot t_p} $$ ### 5.3 Ablation Threshold #### 5.3.1 Single-Shot Ablation Threshold $$ F_{th} = \frac{\rho \cdot C_p \cdot (T_m - T_0) + \rho \cdot L_m}{\alpha \cdot (1-R)} $$ **Where:** - $\rho$ = Density - $C_p$ = Specific heat capacity - $T_m$ = Melting temperature - $T_0$ = Initial temperature - $L_m$ = Latent heat of melting - $\alpha$ = Absorption coefficient - $R$ = Reflectivity #### 5.3.2 Multi-Pulse Incubation Effect $$ F_{th}(N) = F_{th}(1) \cdot N^{S-1} $$ **Where:** - $N$ = Number of pulses - $S$ = Incubation coefficient ($S < 1$ for incubation) ### 5.4 Ablation Depth Model #### 5.4.1 Logarithmic Ablation Law $$ d = \frac{1}{\alpha_{eff}} \ln\left( \frac{F_0}{F_{th}} \right) $$ **Where:** - $d$ = Ablation depth per pulse - $\alpha_{eff}$ = Effective absorption coefficient #### 5.4.2 Total Ablation Depth $$ D_{total} = N_{eff} \cdot d = \frac{N_{eff}}{\alpha_{eff}} \ln\left( \frac{F_0}{F_{th}} \right) $$ ### 5.5 Heat-Affected Zone (HAZ) Model #### 5.5.1 HAZ Definition Criterion $$ \text{HAZ}: \quad T_{max}(r,z) \geq T_{melt} $$ For silicon: $T_{melt} = 1685$ K #### 5.5.2 HAZ Width vs Pulse Duration Empirical relationship: $$ W_{HAZ} \approx k \cdot t_p^{0.4} $$ | Pulse Duration | HAZ Width | |----------------|-----------| | 1 ns | ~0.3 $\mu$m | | 50 ns | ~1.0 $\mu$m | | 200 ns | ~1.7 $\mu$m | ### 5.6 Net Fluence (Accumulated Energy) $$ F_{net} = \frac{P_{avg}}{v_{scan} \cdot w_0} \cdot N_{passes} $$ **Where:** - $P_{avg}$ = Average laser power - $v_{scan}$ = Scanning velocity - $w_0$ = Beam waist - $N_{passes}$ = Number of scan passes ### 5.7 Pulse Overlap #### 5.7.1 Spatial Overlap $$ O_s = 1 - \frac{v_{scan}}{f_{rep} \cdot 2w_0} $$ **Where:** - $O_s$ = Spatial overlap (0 to 1) - $f_{rep}$ = Pulse repetition rate #### 5.7.2 Effective Pulse Number $$ N_{eff} = \frac{f_{rep} \cdot 2w_0}{v_{scan}} $$ ### 5.8 Die Breaking Strength Correlation $$ \sigma_{break} = A - B \cdot \ln(F_{net}) $$ **Where:** - $A, B$ = Material-dependent constants - Higher $F_{net}$ → Lower breaking strength due to HAZ ## 6. Stealth Dicing Mathematical Framework ### 6.1 Principle Overview Stealth dicing uses IR laser focused inside silicon to create internal modified layers without surface damage. **Key physics:** - Two-photon absorption (TPA) - Temperature-dependent absorption coefficient - Internal stress generation - Controlled crack propagation ### 6.2 Nonlinear Absorption Model #### 6.2.1 Two-Photon Absorption $$ \frac{dI}{dz} = -\alpha I - \beta I^2 $$ **Where:** - $\alpha$ = Linear absorption coefficient - $\beta$ = Two-photon absorption coefficient - $I$ = Intensity #### 6.2.2 Temperature-Dependent Absorption $$ \alpha(T) = \alpha_0 \cdot \exp\left( \frac{T - T_0}{T_c} \right) $$ **Where:** - $\alpha_0$ = Absorption coefficient at reference temperature - $T_c$ = Characteristic temperature **Critical relationship:** As temperature increases, silicon's bandgap shrinks: $$ E_g(T) = E_g(0) - \frac{\alpha_{SV} \cdot T^2}{T + \beta_{SV}} $$ ### 6.3 Modified Layer Formation #### 6.3.1 Energy Density in Modified Region $$ E_v = \frac{E_p \cdot \eta}{V_{mod}} $$ **Where:** - $E_v$ = Volumetric energy density (kJ/cm³) - $E_p$ = Pulse energy - $\eta$ = Absorption efficiency - $V_{mod}$ = Modified volume **Typical range:** $0.5 \leq E_v \leq 4.6$ kJ/cm³ #### 6.3.2 Modified Layer Dimensions $$ L_{mod} = f(E_p, z_f, NA) $$ **Where:** - $L_{mod}$ = Modified layer length - $E_p$ = Pulse energy - $z_f$ = Focal depth - $NA$ = Numerical aperture ### 6.4 Pulse Number per Spot $$ N = \frac{f_{rep} \cdot S}{v_{scan}} $$ **Where:** - $N$ = Number of pulses per spot - $f_{rep}$ = Repetition rate - $S$ = Focal spot size - $v_{scan}$ = Scanning speed ### 6.5 Heat Conduction Model #### 6.5.1 3D Heat Equation $$ \rho C_p \frac{\partial T}{\partial t} = k \nabla^2 T + Q(x,y,z,t) $$ **Where:** - $Q$ = Volumetric heat source from laser absorption #### 6.5.2 Gaussian Heat Source $$ Q(r,z,t) = \frac{\alpha P}{\pi w^2} \exp\left( -\frac{2r^2}{w^2} \right) \exp(-\alpha z) \cdot g(t) $$ ### 6.6 Thermal Stress Model #### 6.6.1 Thermoelastic Stress $$ \sigma_{th} = \frac{E \alpha_{th} \Delta T}{1 - \nu} $$ **Where:** - $E$ = Young's modulus - $\alpha_{th}$ = Thermal expansion coefficient - $\nu$ = Poisson's ratio - $\Delta T$ = Temperature change #### 6.6.2 Stress Intensity Factor from Modified Layer $$ K_I = \sigma_{th} \sqrt{\pi a_{mod}} \cdot f(geometry) $$ ### 6.7 Crack Propagation Criterion Crack propagates when: $$ K_I \geq K_{IC} $$ **Crack propagation direction** follows maximum tangential stress criterion: $$ \theta_c = 2 \arctan\left( \frac{1}{4} \left[ \frac{K_I}{K_{II}} - \text{sign}(K_{II})\sqrt{\left(\frac{K_I}{K_{II}}\right)^2 + 8} \right] \right) $$ ## 7. Plasma Dicing (DRIE/Bosch Process) Models ### 7.1 Bosch Process Overview Three-step cyclic process: 1. **Passivation deposition** (C₄F₈) 2. **Bottom film etching** (SF₆ + ion bombardment) 3. **Silicon etching** (SF₆ radicals) ### 7.2 Etch Rate Fundamentals #### 7.2.1 Chemical Etch Rate $$ R_{chem} = k_0 \cdot [F^*] \cdot \exp\left( -\frac{E_a}{k_B T} \right) $$ **Where:** - $k_0$ = Pre-exponential factor - $[F^*]$ = Fluorine radical concentration - $E_a$ = Activation energy - $k_B$ = Boltzmann constant - $T$ = Surface temperature #### 7.2.2 Ion-Enhanced Etch Rate $$ R_{total} = R_{chem} + R_{ion} = R_{chem} \cdot (1 + \gamma \cdot \Gamma_i \cdot \sqrt{E_i}) $$ **Where:** - $\gamma$ = Ion enhancement factor - $\Gamma_i$ = Ion flux - $E_i$ = Ion energy ### 7.3 Loading Effect Model $$ R(A_{open}) = R_0 \cdot \left( 1 + \frac{A_{open}}{A_{ref}} \right)^{-n} $$ **Where:** - $R_0$ = Baseline etch rate - $A_{open}$ = Open area (exposed silicon) - $A_{ref}$ = Reference area - $n$ = Loading effect exponent ### 7.4 Aspect Ratio Dependent Etching (ARDE) #### 7.4.1 Knudsen Transport Model $$ \frac{R(AR)}{R_0} = \frac{1}{1 + \frac{AR}{AR_{crit}}} $$ **Where:** - $AR$ = Aspect ratio = depth/width - $AR_{crit}$ = Critical aspect ratio #### 7.4.2 RIE Lag $$ \Delta d = d_0 \cdot \left( 1 - \frac{w_{narrow}}{w_{wide}} \right)^{\beta} $$ **Where:** - $\Delta d$ = Depth difference - $w_{narrow}, w_{wide}$ = Feature widths - $\beta$ = Lag exponent ### 7.5 Selectivity $$ S = \frac{R_{Si}}{R_{mask}} $$ **Typical values:** | Mask Material | Selectivity | |---------------|-------------| | Photoresist | 50-100:1 | | SiO₂ | 100-200:1 | | Al | >200:1 | ### 7.6 Scallop Formation Model Scallop depth per cycle: $$ \delta_{scallop} = R_{Si} \cdot t_{etch} - R_{pass} \cdot t_{pass} $$ **Where:** - $t_{etch}$ = Etch step duration - $t_{pass}$ = Passivation step duration **Scallop pitch:** $$ p_{scallop} = R_{vertical} \cdot (t_{etch} + t_{pass}) $$ ### 7.7 Sidewall Angle $$ \theta = 90° - \arctan\left( \frac{R_{lateral}}{R_{vertical}} \right) $$ **Process control:** - Increase passivation → More vertical (closer to 90°) - Increase etch time → More tapered ### 7.8 Etch Depth Calculation $$ d_{total} = N_{cycles} \cdot d_{cycle} $$ **Where:** - $N_{cycles}$ = Number of Bosch cycles - $d_{cycle}$ = Depth per cycle (typically 0.5-5 $\mu$m) ## 8. Process Comparison ### 8.1 Key Performance Metrics | Parameter | Blade Dicing | Laser Ablation | Stealth Dicing | Plasma Dicing | |-----------|--------------|----------------|----------------|---------------| | **Kerf Width** | 20-100 $\mu$m | 10-30 $\mu$m | ~0 $\mu$m | <10 $\mu$m | | **Die Strength** | Lower | Medium | High | Highest | | **Throughput** | Medium | High | High | Variable | | **HAZ** | Minimal | Present | Internal only | None | | **Debris** | Yes (wet) | Yes | No | No | | **Process** | Wet | Dry | Dry | Dry | ### 8.2 Die Count Advantage (Plasma vs Blade) $$ \Delta N = \frac{\pi D^2}{4} \cdot \left( \frac{1}{(s+k_p)^2} - \frac{1}{(s+k_b)^2} \right) $$ **Where:** - $k_p$ = Plasma kerf (~5 $\mu$m) - $k_b$ = Blade kerf (~50 $\mu$m) - $s$ = Die size ### 8.3 Cost Model $$ C_{total} = C_{equipment} + C_{consumables} + C_{yield\_loss} $$ **Where:** $$ C_{yield\_loss} = (1 - Y) \cdot N_{dies} \cdot C_{die} $$ ## 9. Advanced Modeling Considerations ### 9.1 AI-Based Process Optimization Machine learning models for process optimization: $$ \vec{y} = f_{ML}(\vec{x}; \vec{\theta}) $$ **Where:** - $\vec{x}$ = Input parameters (speed, power, pressure, etc.) - $\vec{y}$ = Output metrics (kerf, chipping, die strength) - $\vec{\theta}$ = Model parameters ### 9.2 Finite Element Analysis (FEA) #### 9.2.1 Governing Equations **Mechanical equilibrium:** $$ \nabla \cdot \boldsymbol{\sigma} + \vec{f} = \rho \frac{\partial^2 \vec{u}}{\partial t^2} $$ **Thermal diffusion:** $$ \rho C_p \frac{\partial T}{\partial t} = \nabla \cdot (k \nabla T) + Q $$ #### 9.2.2 Coupled Thermo-Mechanical $$ \boldsymbol{\sigma} = \mathbf{C} : (\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}_{th}) $$ **Where:** $$ \boldsymbol{\varepsilon}_{th} = \alpha_{th} (T - T_{ref}) \mathbf{I} $$ ### 9.3 Extended Finite Element Method (XFEM) For crack propagation modeling: $$ u^h(\vec{x}) = \sum_i N_i(\vec{x}) u_i + \sum_j N_j(\vec{x}) H(\vec{x}) a_j + \sum_k N_k(\vec{x}) \sum_{l=1}^{4} F_l(\vec{x}) b_k^l $$ **Where:** - $H(\vec{x})$ = Heaviside enrichment function - $F_l(\vec{x})$ = Crack tip enrichment functions ### 9.4 Hybrid Dicing Optimization Objective function: $$ \min_{\vec{p}} \left[ w_1 \cdot C_{chipping} + w_2 \cdot C_{HAZ} + w_3 \cdot C_{time} \right] $$ Subject to constraints: - $K_{eff} \leq K_{max}$ - $\sigma_{break} \geq \sigma_{min}$ - $C_{chip} \leq C_{max}$ ## 10. Key Equations ### 10.1 Yield and Geometry | Model | Equation | |-------|----------| | Dies per Wafer | $\text{DPW} = \frac{\pi (D/2)^2}{A_{die} + A_{scribe}} \times U$ | | Poisson Yield | $Y = e^{-D_0 \cdot A}$ | | Murphy Yield | $Y = \left( \frac{1 - e^{-DA}}{DA} \right)^2$ | ### 10.2 Mechanical Dicing | Model | Equation | |-------|----------| | Kerf Width | $K_{eff} = K_{blade} + 2\delta_{wear} + \delta_{runout}$ | | Cutting Force | $F_t = K_s \cdot a_p \cdot f$ | | MRR | $\text{MRR} = v_f \cdot a_p \cdot K_{eff}$ | ### 10.3 Fracture Mechanics | Model | Equation | |-------|----------| | Stress Intensity Factor | $K_I = \sigma \sqrt{\pi a} \cdot f(a/W)$ | | Critical Crack Length | $a_c = \frac{1}{\pi} \left( \frac{K_C b}{6MY} \right)^2 h^4$ | | Die Strength | $\sigma_{break} = \frac{3FL}{2bt^2}$ | ### 10.4 Laser Processing | Model | Equation | |-------|----------| | Peak Fluence | $F_0 = \frac{2E_p}{\pi w_0^2}$ | | Thermal Diffusion Length | $L_{th} = \sqrt{D_{th} \cdot t_p}$ | | Ablation Depth | $d = \frac{1}{\alpha_{eff}} \ln\left(\frac{F_0}{F_{th}}\right)$ | | Net Fluence | $F_{net} = \frac{P_{avg}}{v_{scan} \cdot w_0} \cdot N_{passes}$ | ### 10.5 Plasma Dicing | Model | Equation | |-------|----------| | Etch Rate | $R = k_0 [F^*] \exp(-E_a/k_BT)$ | | ARDE | $R(AR)/R_0 = 1/(1 + AR/AR_{crit})$ | | Selectivity | $S = R_{Si}/R_{mask}$ | ### Material Properties (Silicon) | Property | Value | Units | |----------|-------|-------| | Young's Modulus (⟨100⟩) | 130 | GPa | | Young's Modulus (⟨111⟩) | 188 | GPa | | Fracture Toughness {111} | 0.82 | MPa$\cdot$√m | | Thermal Conductivity | 149 | W/(m$\cdot$K) | | Thermal Diffusivity | 0.8 | cm²/s | | Melting Point | 1685 | K | | Density | 2330 | kg/m³ | ## Notation Summary | Symbol | Description | Units | |--------|-------------|-------| | $D$ | Wafer diameter | mm | | $A$ | Area | mm² | | $K$ | Kerf width | $\mu$m | | $v_f$ | Feed rate | mm/s | | $v_s$ | Blade velocity | m/s | | $F$ | Force | N | | $\sigma$ | Stress | MPa | | $K_I$ | Stress intensity factor | MPa$\cdot$√m | | $K_{IC}$ | Fracture toughness | MPa$\cdot$√m | | $E_p$ | Pulse energy | J | | $F_0$ | Peak fluence | J/cm² | | $t_p$ | Pulse duration | ns, ps, fs | | $R$ | Etch rate | $\mu$m/min | | $Y$ | Yield | fraction |
Map electric fields in materials.
Measure overlay using diffraction gratings.
# Mathematical Modeling of Diffusion and Ion Implantation in Semiconductor Manufacturing Part I: Diffusion Modeling Fundamental Equations Dopant redistribution in silicon at elevated temperatures is governed by Fick's Laws . Fick's First Law Relates flux to concentration gradient: $$ J = -D \frac{\partial C}{\partial x} $$ Where: - $J$ — Atomic flux (atoms/cm²·s) - $D$ — Diffusion coefficient (cm²/s) - $C$ — Concentration (atoms/cm³) - $x$ — Position (cm) Fick's Second Law The diffusion equation follows from continuity: $$ \frac{\partial C}{\partial t} = D \frac{\partial^2 C}{\partial x^2} $$ This parabolic PDE admits analytical solutions for idealized boundary conditions. Temperature Dependence The diffusion coefficient follows an Arrhenius relationship : $$ D(T) = D_0 \exp\left(-\frac{E_a}{kT}\right) $$ Parameters: - $D_0$ — Pre-exponential factor (cm²/s) - $E_a$ — Activation energy (eV) - $k$ — Boltzmann's constant ($8.617 \times 10^{-5}$ eV/K) - $T$ — Absolute temperature (K) Typical Values for Phosphorus in Silicon: | Parameter | Value | |-----------|-------| | $D_0$ | $3.85$ cm²/s | | $E_a$ | $3.66$ eV | Diffusion approximately doubles every 10–15°C near typical process temperatures (900–1100°C). Classical Analytical Solutions Case 1: Constant Surface Concentration (Predeposition) Boundary Conditions: - $C(0, t) = C_s$ (constant surface concentration) - $C(\infty, t) = 0$ (zero at infinite depth) - $C(x, 0) = 0$ (initially undoped) Solution: $$ C(x,t) = C_s \cdot \text{erfc}\left(\frac{x}{2\sqrt{Dt}}\right) $$ Complementary Error Function: $$ \text{erfc}(z) = 1 - \text{erf}(z) = \frac{2}{\sqrt{\pi}} \int_z^{\infty} e^{-u^2} \, du $$ Total Incorporated Dose: $$ Q(t) = \frac{2 C_s \sqrt{Dt}}{\sqrt{\pi}} $$ Case 2: Fixed Dose (Drive-in Diffusion) Boundary Conditions: - $\displaystyle\int_0^{\infty} C \, dx = Q$ (constant total dose) - $\displaystyle\frac{\partial C}{\partial x}\bigg|_{x=0} = 0$ (no flux at surface) Solution (Gaussian Profile): $$ C(x,t) = \frac{Q}{\sqrt{\pi Dt}} \exp\left(-\frac{x^2}{4Dt}\right) $$ Peak Surface Concentration: $$ C(0,t) = \frac{Q}{\sqrt{\pi Dt}} $$ Junction Depth Calculation The metallurgical junction forms where dopant concentration equals background doping $C_B$. For erfc Profile: $$ x_j = 2\sqrt{Dt} \cdot \text{erfc}^{-1}\left(\frac{C_B}{C_s}\right) $$ For Gaussian Profile: $$ x_j = 2\sqrt{Dt \cdot \ln\left(\frac{Q}{C_B \sqrt{\pi Dt}}\right)} $$ Concentration-Dependent Diffusion At high doping concentrations (approaching or exceeding intrinsic carrier concentration $n_i$), diffusivity becomes concentration-dependent. Generalized Model: $$ D = D^0 + D^{-}\frac{n}{n_i} + D^{+}\frac{p}{n_i} + D^{=}\left(\frac{n}{n_i}\right)^2 $$ Physical Interpretation: | Term | Mechanism | |------|-----------| | $D^0$ | Neutral vacancy diffusion | | $D^{-}$ | Singly negative vacancy diffusion | | $D^{+}$ | Positive vacancy diffusion | | $D^{=}$ | Doubly negative vacancy diffusion | Resulting Nonlinear PDE: $$ \frac{\partial C}{\partial t} = \frac{\partial}{\partial x}\left(D(C) \frac{\partial C}{\partial x}\right) $$ This requires numerical solution methods. Point Defect Mediated Diffusion Modern process modeling couples dopant diffusion to point defect dynamics. Governing System of PDEs: $$ \frac{\partial C_I}{\partial t} = \nabla \cdot (D_I \nabla C_I) - k_{IV} C_I C_V + G_I - R_I $$ $$ \frac{\partial C_V}{\partial t} = \nabla \cdot (D_V \nabla C_V) - k_{IV} C_I C_V + G_V - R_V $$ $$ \frac{\partial C_A}{\partial t} = \nabla \cdot (D_{AI} C_I \nabla C_A) + \text{(clustering terms)} $$ Variable Definitions: - $C_I$ — Interstitial concentration - $C_V$ — Vacancy concentration - $C_A$ — Dopant atom concentration - $k_{IV}$ — Interstitial-vacancy recombination rate - $G$ — Generation rate - $R$ — Surface recombination rate Part II: Ion Implantation Modeling Energy Loss Mechanisms Implanted ions lose energy through two mechanisms: Total Stopping Power: $$ S(E) = -\frac{dE}{dx} = S_n(E) + S_e(E) $$ Nuclear Stopping (Elastic Collisions) Dominates at low energies : $$ S_n(E) = \frac{\pi a^2 \gamma E \cdot s_n(\varepsilon)}{1 + M_2/M_1} $$ Where: - $\gamma = \displaystyle\frac{4 M_1 M_2}{(M_1 + M_2)^2}$ — Energy transfer factor - $a$ — Screening length - $s_n(\varepsilon)$ — Reduced nuclear stopping Electronic Stopping (Inelastic Interactions) Dominates at high energies : $$ S_e(E) \propto \sqrt{E} $$ (at intermediate energies) LSS Theory Lindhard, Scharff, and Schiøtt developed universal scaling using reduced units. Reduced Energy: $$ \varepsilon = \frac{a M_2 E}{Z_1 Z_2 e^2 (M_1 + M_2)} $$ Reduced Path Length: $$ \rho = 4\pi a^2 N \frac{M_1 M_2}{(M_1 + M_2)^2} \cdot x $$ This allows tabulation of universal range curves applicable across ion-target combinations. Gaussian Profile Approximation First-Order Implant Profile: $$ C(x) = \frac{\Phi}{\sqrt{2\pi} \, \Delta R_p} \exp\left(-\frac{(x - R_p)^2}{2 \Delta R_p^2}\right) $$ Parameters: | Symbol | Name | Units | |--------|------|-------| | $\Phi$ | Dose | ions/cm² | | $R_p$ | Projected range (mean stopping depth) | cm | | $\Delta R_p$ | Range straggle (standard deviation) | cm | Peak Concentration: $$ C_{\text{peak}} = \frac{\Phi}{\sqrt{2\pi} \, \Delta R_p} \approx \frac{0.4 \, \Phi}{\Delta R_p} $$ Higher-Order Moment Distributions The Gaussian approximation fails for many practical cases. The Pearson IV distribution uses four statistical moments: | Moment | Symbol | Physical Meaning | |--------|--------|------------------| | 1st | $R_p$ | Projected range | | 2nd | $\Delta R_p$ | Range straggle | | 3rd | $\gamma$ | Skewness | | 4th | $\beta$ | Kurtosis | Pearson IV Form: $$ C(x) = \frac{K}{\left[(x-a)^2 + b^2\right]^m} \exp\left(-\nu \arctan\frac{x-a}{b}\right) $$ Parameters $(a, b, m, \nu, K)$ are derived from the four moments through algebraic relations. Skewness Behavior: - Light ions (B) in heavy substrates → Negative skewness (tail toward surface) - Heavy ions (As, Sb) in silicon → Positive skewness (tail toward bulk) Dual Pearson Model For channeling tails or complex profiles: $$ C(x) = f \cdot C_1(x) + (1-f) \cdot C_2(x) $$ Where: - $C_1(x)$, $C_2(x)$ — Two Pearson distributions with different parameters - $f$ — Weight fraction Lateral Distribution Ions scatter laterally as well: $$ C(x, r) = C(x) \cdot \frac{1}{2\pi \Delta R_{\perp}^2} \exp\left(-\frac{r^2}{2 \Delta R_{\perp}^2}\right) $$ For Amorphous Targets: $$ \Delta R_{\perp} \approx \frac{\Delta R_p}{\sqrt{3}} $$ Lateral straggle is critical for device scaling—it limits minimum feature sizes. Monte Carlo Simulation (TRIM/SRIM) For accurate profiles, especially in multilayer or crystalline structures, Monte Carlo methods track individual ion trajectories. Algorithm: 1. Initialize ion position, direction, energy 2. Select free flight path: $\lambda = 1/(N\pi a^2)$ 3. Calculate impact parameter and scattering angle via screened Coulomb potential 4. Energy transfer to recoil: $$T = T_m \sin^2\left(\frac{\theta}{2}\right)$$ where $T_m = \gamma E$ 5. Apply electronic energy loss over path segment 6. Update ion position/direction; cascade recoils if $T > E_d$ (displacement energy) 7. Repeat until $E < E_{\text{cutoff}}$ 8. Accumulate statistics over $10^4 - 10^6$ ion histories ZBL Interatomic Potential: $$ V(r) = \frac{Z_1 Z_2 e^2}{r} \, \phi(r/a) $$ Where $\phi$ is the screening function tabulated from quantum mechanical calculations. Channeling In crystalline silicon, ions aligned with crystal axes experience reduced stopping. Critical Angle for Channeling: $$ \psi_c \approx \sqrt{\frac{2 Z_1 Z_2 e^2}{E \, d}} $$ Where: - $d$ — Atomic spacing along the channel - $E$ — Ion energy Effects: - Channeled ions penetrate 2–10× deeper - Creates extended tails in profiles - Modern implants use 7° tilt or random-equivalent conditions to minimize Damage Accumulation Implant damage is quantified by: $$ D(x) = \Phi \int_0^{\infty} \nu(E) \cdot F(x, E) \, dE $$ Where: - $\nu(E)$ — Kinchin-Pease damage function (displaced atoms per ion) - $F(x, E)$ — Energy deposition profile Amorphization Threshold for Silicon: $$ \sim 10^{22} \text{ displacements/cm}^3 $$ (approximately 10–15% of atoms displaced) Part III: Post-Implant Diffusion and Transient Enhanced Diffusion Transient Enhanced Diffusion (TED) After implantation, excess interstitials dramatically enhance diffusion until they anneal: $$ D_{\text{eff}} = D^* \left(1 + \frac{C_I}{C_I^*}\right) $$ Where: - $C_I^*$ — Equilibrium interstitial concentration "+1" Model for Boron: $$ \frac{\partial C_B}{\partial t} = \frac{\partial}{\partial x}\left[D_B \left(1 + \frac{C_I}{C_I^*}\right) \frac{\partial C_B}{\partial x}\right] $$ Impact: TED can cause junction depths 2–5× deeper than equilibrium diffusion would predict—critical for modern shallow junctions. {311} Defect Dissolution Kinetics Interstitials cluster into rod-like {311} defects that slowly dissolve: $$ \frac{dN_{311}}{dt} = -\nu_0 \exp\left(-\frac{E_a}{kT}\right) N_{311} $$ The released interstitials sustain TED, explaining why TED persists for times much longer than point defect diffusion would suggest. Part IV: Numerical Methods Finite Difference Discretization For the diffusion equation on uniform grid $(x_i, t_n)$: Explicit (Forward Euler) $$ \frac{C_i^{n+1} - C_i^n}{\Delta t} = D \frac{C_{i+1}^n - 2C_i^n + C_{i-1}^n}{\Delta x^2} $$ Stability Requirement (CFL Condition): $$ \Delta t < \frac{\Delta x^2}{2D} $$ Implicit (Backward Euler) $$ \frac{C_i^{n+1} - C_i^n}{\Delta t} = D \frac{C_{i+1}^{n+1} - 2C_i^{n+1} + C_{i-1}^{n+1}}{\Delta x^2} $$ - Unconditionally stable - Requires solving tridiagonal system each timestep Crank-Nicolson Method - Average of explicit and implicit schemes - Second-order accurate in time - Results in tridiagonal system Adaptive Meshing Concentration gradients vary by orders of magnitude. Adaptive grids refine near: - Junctions - Surface - Implant peaks - Moving interfaces Grid Spacing Scaling: $$ \Delta x \propto \frac{C}{|\nabla C|} $$ Process Simulation Flow (TCAD) Modern simulators (Sentaurus Process, ATHENA, FLOOPS) integrate: 1. Implantation → Monte Carlo or analytical tables 2. Damage model → Amorphization, defect clustering 3. Annealing → Coupled dopant-defect PDEs 4. Oxidation → Deal-Grove kinetics, stress effects, OED 5. Silicidation, epitaxy, etc. → Specialized models Output feeds device simulation (drift-diffusion, Monte Carlo transport). Part V: Key Process Design Equations Thermal Budget The characteristic diffusion length after multiple thermal steps: $$ \sqrt{Dt}_{\text{total}} = \sqrt{\sum_i D_i t_i} $$ For Varying Temperature $T(t)$: $$ Dt = \int_0^{t_f} D_0 \exp\left(-\frac{E_a}{kT(t')}\right) dt' $$ Sheet Resistance $$ R_s = \frac{1}{q \displaystyle\int_0^{x_j} \mu(C) \cdot C(x) \, dx} $$ For Uniform Mobility Approximation: $$ R_s \approx \frac{1}{q \mu Q} $$ Electrical measurements to profile parameters. Implant Dose-Energy Selection Target Peak Concentration: $$ C_{\text{peak}} = \frac{0.4 \, \Phi}{\Delta R_p(E)} $$ Target Depth (Empirical): $$ R_p(E) \approx A \cdot E^n $$ Where: - $n \approx 0.6 - 0.8$ (depending on energy regime) - $A$ — Ion-target dependent constant Key Mathematical Tools: | Process | Core Equation | Solution Method | |---------|---------------|-----------------| | Thermal diffusion | $\displaystyle\frac{\partial C}{\partial t} = \nabla \cdot (D \nabla C)$ | Analytical (erfc, Gaussian) or FEM/FDM | | Implant profile | 4-moment Pearson distribution | Lookup tables or Monte Carlo | | Damage evolution | Coupled defect-dopant kinetics | Stiff ODE solvers | | TED | $D_{\text{eff}} = D^*(1 + C_I/C_I^*)$ | Coupled PDEs | | 2D/3D profiles | $\nabla \cdot (D \nabla C)$ in 2D/3D | Finite element methods | Common Dopant Properties in Silicon: | Dopant | Type | $D_0$ (cm²/s) | $E_a$ (eV) | Typical Use | |--------|------|---------------|------------|-------------| | Boron (B) | p-type | 0.76 | 3.46 | Source/drain, channel doping | | Phosphorus (P) | n-type | 3.85 | 3.66 | Source/drain, n-well | | Arsenic (As) | n-type | 0.32 | 3.56 | Shallow junctions | | Antimony (Sb) | n-type | 0.214 | 3.65 | Buried layers |
# Mathematical Modeling of Diffusion 1. Fundamental Governing Equations 1.1 Fick's Laws of Diffusion The foundation of diffusion modeling in semiconductor manufacturing rests on Fick's laws : Fick's First Law The flux is proportional to the concentration gradient: $$ J = -D \frac{\partial C}{\partial x} $$ Where: - $J$ = flux (atoms/cm²·s) - $D$ = diffusion coefficient (cm²/s) - $C$ = concentration (atoms/cm³) - $x$ = position (cm) Note: The negative sign indicates diffusion occurs from high to low concentration regions. Fick's Second Law Derived from the continuity equation combined with Fick's first law: $$ \frac{\partial C}{\partial t} = D \frac{\partial^2 C}{\partial x^2} $$ Key characteristics: - This is a parabolic partial differential equation - Mathematically identical to the heat equation - Assumes constant diffusion coefficient $D$ 1.2 Temperature Dependence (Arrhenius Relationship) The diffusion coefficient follows the Arrhenius relationship: $$ D(T) = D_0 \exp\left(-\frac{E_a}{kT}\right) $$ Where: - $D_0$ = pre-exponential factor (cm²/s) - $E_a$ = activation energy (eV) - $k$ = Boltzmann constant ($8.617 \times 10^{-5}$ eV/K) - $T$ = absolute temperature (K) 1.3 Typical Dopant Parameters in Silicon | Dopant | $D_0$ (cm²/s) | $E_a$ (eV) | $D$ at 1100°C (cm²/s) | |--------|---------------|------------|------------------------| | Boron (B) | ~10.5 | ~3.69 | ~$10^{-13}$ | | Phosphorus (P) | ~10.5 | ~3.69 | ~$10^{-13}$ | | Arsenic (As) | ~0.32 | ~3.56 | ~$10^{-14}$ | | Antimony (Sb) | ~5.6 | ~3.95 | ~$10^{-14}$ | 2. Analytical Solutions for Standard Boundary Conditions 2.1 Constant Surface Concentration (Predeposition) Boundary and Initial Conditions - $C(0,t) = C_s$ — surface held at solid solubility - $C(x,0) = 0$ — initially undoped wafer - $C(\infty,t) = 0$ — semi-infinite substrate Solution: Complementary Error Function Profile $$ C(x,t) = C_s \cdot \text{erfc}\left(\frac{x}{2\sqrt{Dt}}\right) $$ Where the complementary error function is defined as: $$ \text{erfc}(\eta) = 1 - \text{erf}(\eta) = 1 - \frac{2}{\sqrt{\pi}}\int_0^\eta e^{-u^2} \, du $$ Total Dose Introduced $$ Q = \int_0^\infty C(x,t) \, dx = \frac{2 C_s \sqrt{Dt}}{\sqrt{\pi}} \approx 1.13 \, C_s \sqrt{Dt} $$ Key Properties - Surface concentration remains constant at $C_s$ - Profile penetrates deeper with increasing $\sqrt{Dt}$ - Characteristic diffusion length: $L_D = 2\sqrt{Dt}$ 2.2 Fixed Dose / Gaussian Drive-in Boundary and Initial Conditions - Total dose $Q$ is conserved (no dopant enters or leaves) - Zero flux at surface: $\left.\frac{\partial C}{\partial x}\right|_{x=0} = 0$ - Delta-function or thin layer initial condition Solution: Gaussian Profile $$ C(x,t) = \frac{Q}{\sqrt{\pi Dt}} \exp\left(-\frac{x^2}{4Dt}\right) $$ Time-Dependent Surface Concentration $$ C_s(t) = C(0,t) = \frac{Q}{\sqrt{\pi Dt}} $$ Key characteristics: - Surface concentration decreases with time as $t^{-1/2}$ - Profile broadens while maintaining total dose - Peak always at surface ($x = 0$) 2.3 Junction Depth Calculation The junction depth $x_j$ is the position where dopant concentration equals background concentration $C_B$: For erfc Profile $$ x_j = 2\sqrt{Dt} \cdot \text{erfc}^{-1}\left(\frac{C_B}{C_s}\right) $$ For Gaussian Profile $$ x_j = 2\sqrt{Dt \cdot \ln\left(\frac{Q}{C_B \sqrt{\pi Dt}}\right)} $$ 3. Green's Function Method 3.1 General Solution for Arbitrary Initial Conditions For an arbitrary initial profile $C_0(x')$, the solution is a convolution with the Gaussian kernel (Green's function): $$ C(x,t) = \int_{-\infty}^{\infty} C_0(x') \cdot \frac{1}{2\sqrt{\pi Dt}} \exp\left(-\frac{(x-x')^2}{4Dt}\right) dx' $$ Physical interpretation: - Each point in the initial distribution spreads as a Gaussian - The final profile is the superposition of all spreading contributions 3.2 Application: Ion-Implanted Gaussian Profile Initial Implant Profile $$ C_0(x) = \frac{Q}{\sqrt{2\pi} \, \Delta R_p} \exp\left(-\frac{(x - R_p)^2}{2 \Delta R_p^2}\right) $$ Where: - $Q$ = implanted dose (atoms/cm²) - $R_p$ = projected range (mean depth) - $\Delta R_p$ = straggle (standard deviation) Profile After Diffusion $$ C(x,t) = \frac{Q}{\sqrt{2\pi \, \sigma_{eff}^2}} \exp\left(-\frac{(x - R_p)^2}{2 \sigma_{eff}^2}\right) $$ Effective Straggle $$ \sigma_{eff} = \sqrt{\Delta R_p^2 + 2Dt} $$ Key observations: - Peak remains at $R_p$ (no shift in position) - Peak concentration decreases - Profile broadens symmetrically 4. Concentration-Dependent Diffusion 4.1 Nonlinear Diffusion Equation At high dopant concentrations (above intrinsic carrier concentration $n_i$), diffusion becomes concentration-dependent : $$ \frac{\partial C}{\partial t} = \frac{\partial}{\partial x}\left(D(C) \frac{\partial C}{\partial x}\right) $$ 4.2 Concentration-Dependent Diffusivity Models Simple Power Law Model $$ D(C) = D^i \left(1 + \left(\frac{C}{n_i}\right)^r\right) $$ Charged Defect Model (Fair's Equation) $$ D = D^0 + D^- \frac{n}{n_i} + D^{=} \left(\frac{n}{n_i}\right)^2 + D^+ \frac{p}{n_i} $$ Where: - $D^0$ = neutral defect contribution - $D^-$ = singly negative defect contribution - $D^{=}$ = doubly negative defect contribution - $D^+$ = positive defect contribution - $n, p$ = electron and hole concentrations 4.3 Electric Field Enhancement High concentration gradients create internal electric fields that enhance diffusion: $$ J = -D \frac{\partial C}{\partial x} - \mu C \mathcal{E} $$ For extrinsic conditions with a single dopant species: $$ J = -hD \frac{\partial C}{\partial x} $$ Field enhancement factor: $$ h = 1 + \frac{C}{n + p} $$ - For fully ionized n-type dopant at high concentration: $h \approx 2$ - Results in approximately 2× faster effective diffusion 4.4 Resulting Profile Shapes - Phosphorus: "Kink-and-tail" profile at high concentrations - Arsenic: Box-like profiles due to clustering - Boron: Enhanced tail diffusion in oxidizing ambient 5. Point Defect-Mediated Diffusion 5.1 Diffusion Mechanisms Dopants don't diffuse as isolated atoms—they move via defect complexes : Vacancy Mechanism $$ A + V \rightleftharpoons AV \quad \text{(dopant-vacancy pair forms, diffuses, dissociates)} $$ Interstitial Mechanism $$ A + I \rightleftharpoons AI \quad \text{(dopant-interstitial pair)} $$ Kick-out Mechanism $$ A_s + I \rightleftharpoons A_i \quad \text{(substitutional ↔ interstitial)} $$ 5.2 Effective Diffusivity $$ D_{eff} = D_V \frac{C_V}{C_V^*} + D_I \frac{C_I}{C_I^*} $$ Where: - $D_V, D_I$ = diffusivity via vacancy/interstitial mechanism - $C_V, C_I$ = actual vacancy/interstitial concentrations - $C_V^*, C_I^*$ = equilibrium concentrations Fractional interstitialcy: $$ f_I = \frac{D_I}{D_V + D_I} $$ | Dopant | $f_I$ | Dominant Mechanism | |--------|-------|-------------------| | Boron | ~1.0 | Interstitial | | Phosphorus | ~0.9 | Interstitial | | Arsenic | ~0.4 | Mixed | | Antimony | ~0.02 | Vacancy | 5.3 Coupled Reaction-Diffusion System The full model requires solving coupled PDEs : Dopant Equation $$ \frac{\partial C_A}{\partial t} = \nabla \cdot \left(D_A \frac{C_I}{C_I^*} \nabla C_A\right) $$ Interstitial Balance $$ \frac{\partial C_I}{\partial t} = D_I \nabla^2 C_I + G - k_{IV}\left(C_I C_V - C_I^* C_V^*\right) $$ Vacancy Balance $$ \frac{\partial C_V}{\partial t} = D_V \nabla^2 C_V + G - k_{IV}\left(C_I C_V - C_I^* C_V^*\right) $$ Where: - $G$ = defect generation rate - $k_{IV}$ = bulk recombination rate constant 5.4 Transient Enhanced Diffusion (TED) After ion implantation, excess interstitials cause anomalously rapid diffusion : The "+1" Model: $$ \int_0^\infty (C_I - C_I^*) \, dx \approx \Phi \quad \text{(implant dose)} $$ Enhancement factor: $$ \frac{D_{eff}}{D^*} = \frac{C_I}{C_I^*} \gg 1 \quad \text{(transient)} $$ Key characteristics: - Enhancement decays as interstitials recombine - Time constant: typically 10-100 seconds at 1000°C - Critical for shallow junction formation 6. Oxidation Effects 6.1 Oxidation-Enhanced Diffusion (OED) During thermal oxidation, silicon interstitials are injected into the substrate: $$ \frac{C_I}{C_I^*} = 1 + A \left(\frac{dx_{ox}}{dt}\right)^n $$ Effective diffusivity: $$ D_{eff} = D^* \left[1 + f_I \left(\frac{C_I}{C_I^*} - 1\right)\right] $$ Dopants enhanced by oxidation: - Boron (high $f_I$) - Phosphorus (high $f_I$) 6.2 Oxidation-Retarded Diffusion (ORD) Growing oxide absorbs vacancies , reducing vacancy concentration: $$ \frac{C_V}{C_V^*} < 1 $$ Dopants retarded by oxidation: - Antimony (low $f_I$, primarily vacancy-mediated) 6.3 Segregation at SiO₂/Si Interface Dopants redistribute at the interface according to the segregation coefficient : $$ m = \frac{C_{Si}}{C_{SiO_2}}\bigg|_{\text{interface}} $$ | Dopant | Segregation Coefficient $m$ | Behavior | |--------|----------------------------|----------| | Boron | ~0.3 | Pile-down (into oxide) | | Phosphorus | ~10 | Pile-up (into silicon) | | Arsenic | ~10 | Pile-up | 7. Numerical Methods 7.1 Finite Difference Method Discretize space and time on grid $(x_i, t^n)$: Explicit Scheme (FTCS) $$ \frac{C_i^{n+1} - C_i^n}{\Delta t} = D \frac{C_{i+1}^n - 2C_i^n + C_{i-1}^n}{(\Delta x)^2} $$ Rearranged: $$ C_i^{n+1} = C_i^n + \alpha \left(C_{i+1}^n - 2C_i^n + C_{i-1}^n\right) $$ Where Fourier number: $$ \alpha = \frac{D \Delta t}{(\Delta x)^2} $$ Stability requirement (von Neumann analysis): $$ \alpha \leq \frac{1}{2} $$ Implicit Scheme (BTCS) $$ \frac{C_i^{n+1} - C_i^n}{\Delta t} = D \frac{C_{i+1}^{n+1} - 2C_i^{n+1} + C_{i-1}^{n+1}}{(\Delta x)^2} $$ - Unconditionally stable (no restriction on $\alpha$) - Requires solving tridiagonal system at each time step Crank-Nicolson Scheme (Second-Order Accurate) $$ C_i^{n+1} - C_i^n = \frac{\alpha}{2}\left[(C_{i+1}^{n+1} - 2C_i^{n+1} + C_{i-1}^{n+1}) + (C_{i+1}^n - 2C_i^n + C_{i-1}^n)\right] $$ Properties: - Unconditionally stable - Second-order accurate in both space and time - Results in tridiagonal system: solved by Thomas algorithm 7.2 Handling Concentration-Dependent Diffusion Use iterative methods: 1. Estimate $D^{(k)}$ from current concentration $C^{(k)}$ 2. Solve linear diffusion equation for $C^{(k+1)}$ 3. Update diffusivity: $D^{(k+1)} = D(C^{(k+1)})$ 4. Iterate until $\|C^{(k+1)} - C^{(k)}\| < \epsilon$ 7.3 Moving Boundary Problems For oxidation with moving Si/SiO₂ interface: Approaches: - Coordinate transformation: Map to fixed domain via $\xi = x/s(t)$ - Front-tracking methods: Explicitly track interface position - Level-set methods: Implicit interface representation - Phase-field methods: Diffuse interface approximation 8. Thermal Budget Concept 8.1 The Dt Product Diffusion profiles scale with $\sqrt{Dt}$. The thermal budget quantifies total diffusion: $$ (Dt)_{total} = \sum_i D(T_i) \cdot t_i $$ 8.2 Continuous Temperature Profile For time-varying temperature: $$ (Dt)_{eff} = \int_0^{t_{total}} D(T(\tau)) \, d\tau $$ 8.3 Equivalent Time at Reference Temperature $$ t_{eq} = \sum_i t_i \exp\left(\frac{E_a}{k}\left(\frac{1}{T_{ref}} - \frac{1}{T_i}\right)\right) $$ 8.4 Combining Multiple Diffusion Steps For sequential Gaussian redistributions: $$ \sigma_{final} = \sqrt{\sum_i 2D_i t_i} $$ For erfc profiles, use effective $(Dt)_{total}$: $$ C(x) = C_s \cdot \text{erfc}\left(\frac{x}{2\sqrt{(Dt)_{total}}}\right) $$ 9. Key Dimensionless Parameters | Parameter | Definition | Physical Meaning | |-----------|------------|------------------| | Fourier Number | $Fo = \dfrac{Dt}{L^2}$ | Diffusion time vs. characteristic length | | Damköhler Number | $Da = \dfrac{kL^2}{D}$ | Reaction rate vs. diffusion rate | | Péclet Number | $Pe = \dfrac{vL}{D}$ | Advection (drift) vs. diffusion | | Biot Number | $Bi = \dfrac{hL}{D}$ | Surface transfer vs. bulk diffusion | 10. Process Simulation Software 10.1 Commercial and Research Tools | Simulator | Developer | Key Capabilities | |-----------|-----------|------------------| | Sentaurus Process | Synopsys | Full 3D, atomistic KMC, advanced models | | Athena | Silvaco | Integrated with device simulation (Atlas) | | SUPREM-IV | Stanford | Classic 1D/2D, widely validated | | FLOOPS | U. Florida | Research-oriented, extensible | | Victory Process | Silvaco | Modern 3D process simulation | 10.2 Physical Models Incorporated - Multiple coupled dopant species - Full point-defect dynamics (I, V, clusters) - Stress-dependent diffusion - Cluster nucleation and dissolution - Atomistic kinetic Monte Carlo (KMC) options - Quantum corrections for ultra-shallow junctions Mathematical Modeling Hierarchy: Level 1: Simple Analytical Models $$ \frac{\partial C}{\partial t} = D \frac{\partial^2 C}{\partial x^2} $$ - Constant $D$ - erfc and Gaussian solutions - Junction depth calculations Level 2: Intermediate Complexity $$ \frac{\partial C}{\partial t} = \frac{\partial}{\partial x}\left(D(C) \frac{\partial C}{\partial x}\right) $$ - Concentration-dependent $D$ - Electric field effects - Nonlinear PDEs requiring numerical methods Level 3: Advanced Coupled Models $$ \begin{aligned} \frac{\partial C_A}{\partial t} &= \nabla \cdot \left(D_A \frac{C_I}{C_I^*} \nabla C_A\right) \\[6pt] \frac{\partial C_I}{\partial t} &= D_I \nabla^2 C_I + G - k_{IV}(C_I C_V - C_I^* C_V^*) \end{aligned} $$ - Coupled dopant-defect systems - TED, OED/ORD effects - Process simulators required Level 4: State-of-the-Art - Atomistic kinetic Monte Carlo - Molecular dynamics for interface phenomena - Ab initio calculations for defect properties - Essential for sub-10nm technology nodes Key Insight The fundamental scaling of semiconductor diffusion is governed by $\sqrt{Dt}$, but the effective diffusion coefficient $D$ depends on: - Temperature (Arrhenius) - Concentration (charged defects) - Point defect supersaturation (TED) - Processing ambient (oxidation) - Mechanical stress This complexity requires sophisticated physical models for modern nanometer-scale devices.
How far acid or other species diffuse in resist.
# Mathematical Modeling of Diffusion in Semiconductor Manufacturing ## 1. Fundamental Governing Equations ### 1.1 Fick's Laws of Diffusion The foundation of diffusion modeling in semiconductor manufacturing rests on **Fick's laws**: #### Fick's First Law The flux is proportional to the concentration gradient: $$ J = -D \frac{\partial C}{\partial x} $$ **Where:** - $J$ = flux (atoms/cm²·s) - $D$ = diffusion coefficient (cm²/s) - $C$ = concentration (atoms/cm³) - $x$ = position (cm) > **Note:** The negative sign indicates diffusion occurs from high to low concentration regions. #### Fick's Second Law Derived from the continuity equation combined with Fick's first law: $$ \frac{\partial C}{\partial t} = D \frac{\partial^2 C}{\partial x^2} $$ **Key characteristics:** - This is a **parabolic partial differential equation** - Mathematically identical to the heat equation - Assumes constant diffusion coefficient $D$ ### 1.2 Temperature Dependence (Arrhenius Relationship) The diffusion coefficient follows the Arrhenius relationship: $$ D(T) = D_0 \exp\left(-\frac{E_a}{kT}\right) $$ **Where:** - $D_0$ = pre-exponential factor (cm²/s) - $E_a$ = activation energy (eV) - $k$ = Boltzmann constant ($8.617 \times 10^{-5}$ eV/K) - $T$ = absolute temperature (K) ### 1.3 Typical Dopant Parameters in Silicon | Dopant | $D_0$ (cm²/s) | $E_a$ (eV) | $D$ at 1100°C (cm²/s) | |--------|---------------|------------|------------------------| | Boron (B) | ~10.5 | ~3.69 | ~$10^{-13}$ | | Phosphorus (P) | ~10.5 | ~3.69 | ~$10^{-13}$ | | Arsenic (As) | ~0.32 | ~3.56 | ~$10^{-14}$ | | Antimony (Sb) | ~5.6 | ~3.95 | ~$10^{-14}$ | ## 2. Analytical Solutions for Standard Boundary Conditions ### 2.1 Constant Surface Concentration (Predeposition) #### Boundary and Initial Conditions - $C(0,t) = C_s$ — surface held at solid solubility - $C(x,0) = 0$ — initially undoped wafer - $C(\infty,t) = 0$ — semi-infinite substrate #### Solution: Complementary Error Function Profile $$ C(x,t) = C_s \cdot \text{erfc}\left(\frac{x}{2\sqrt{Dt}}\right) $$ **Where the complementary error function is defined as:** $$ \text{erfc}(\eta) = 1 - \text{erf}(\eta) = 1 - \frac{2}{\sqrt{\pi}}\int_0^\eta e^{-u^2} \, du $$ #### Total Dose Introduced $$ Q = \int_0^\infty C(x,t) \, dx = \frac{2 C_s \sqrt{Dt}}{\sqrt{\pi}} \approx 1.13 \, C_s \sqrt{Dt} $$ #### Key Properties - Surface concentration remains constant at $C_s$ - Profile penetrates deeper with increasing $\sqrt{Dt}$ - Characteristic diffusion length: $L_D = 2\sqrt{Dt}$ ### 2.2 Fixed Dose / Gaussian Drive-in #### Boundary and Initial Conditions - Total dose $Q$ is conserved (no dopant enters or leaves) - Zero flux at surface: $\left.\frac{\partial C}{\partial x}\right|_{x=0} = 0$ - Delta-function or thin layer initial condition #### Solution: Gaussian Profile $$ C(x,t) = \frac{Q}{\sqrt{\pi Dt}} \exp\left(-\frac{x^2}{4Dt}\right) $$ #### Time-Dependent Surface Concentration $$ C_s(t) = C(0,t) = \frac{Q}{\sqrt{\pi Dt}} $$ **Key characteristics:** - Surface concentration **decreases** with time as $t^{-1/2}$ - Profile broadens while maintaining total dose - Peak always at surface ($x = 0$) ### 2.3 Junction Depth Calculation The **junction depth** $x_j$ is the position where dopant concentration equals background concentration $C_B$: #### For erfc Profile $$ x_j = 2\sqrt{Dt} \cdot \text{erfc}^{-1}\left(\frac{C_B}{C_s}\right) $$ #### For Gaussian Profile $$ x_j = 2\sqrt{Dt \cdot \ln\left(\frac{Q}{C_B \sqrt{\pi Dt}}\right)} $$ ## 3. Green's Function Method ### 3.1 General Solution for Arbitrary Initial Conditions For an arbitrary initial profile $C_0(x')$, the solution is a **convolution** with the Gaussian kernel (Green's function): $$ C(x,t) = \int_{-\infty}^{\infty} C_0(x') \cdot \frac{1}{2\sqrt{\pi Dt}} \exp\left(-\frac{(x-x')^2}{4Dt}\right) dx' $$ **Physical interpretation:** - Each point in the initial distribution spreads as a Gaussian - The final profile is the superposition of all spreading contributions ### 3.2 Application: Ion-Implanted Gaussian Profile #### Initial Implant Profile $$ C_0(x) = \frac{Q}{\sqrt{2\pi} \, \Delta R_p} \exp\left(-\frac{(x - R_p)^2}{2 \Delta R_p^2}\right) $$ **Where:** - $Q$ = implanted dose (atoms/cm²) - $R_p$ = projected range (mean depth) - $\Delta R_p$ = straggle (standard deviation) #### Profile After Diffusion $$ C(x,t) = \frac{Q}{\sqrt{2\pi \, \sigma_{eff}^2}} \exp\left(-\frac{(x - R_p)^2}{2 \sigma_{eff}^2}\right) $$ #### Effective Straggle $$ \sigma_{eff} = \sqrt{\Delta R_p^2 + 2Dt} $$ **Key observations:** - Peak remains at $R_p$ (no shift in position) - Peak concentration decreases - Profile broadens symmetrically ## 4. Concentration-Dependent Diffusion ### 4.1 Nonlinear Diffusion Equation At high dopant concentrations (above intrinsic carrier concentration $n_i$), diffusion becomes **concentration-dependent**: $$ \frac{\partial C}{\partial t} = \frac{\partial}{\partial x}\left(D(C) \frac{\partial C}{\partial x}\right) $$ ### 4.2 Concentration-Dependent Diffusivity Models #### Simple Power Law Model $$ D(C) = D^i \left(1 + \left(\frac{C}{n_i}\right)^r\right) $$ #### Charged Defect Model (Fair's Equation) $$ D = D^0 + D^- \frac{n}{n_i} + D^{=} \left(\frac{n}{n_i}\right)^2 + D^+ \frac{p}{n_i} $$ **Where:** - $D^0$ = neutral defect contribution - $D^-$ = singly negative defect contribution - $D^{=}$ = doubly negative defect contribution - $D^+$ = positive defect contribution - $n, p$ = electron and hole concentrations ### 4.3 Electric Field Enhancement High concentration gradients create internal electric fields that enhance diffusion: $$ J = -D \frac{\partial C}{\partial x} - \mu C \mathcal{E} $$ For extrinsic conditions with a single dopant species: $$ J = -hD \frac{\partial C}{\partial x} $$ **Field enhancement factor:** $$ h = 1 + \frac{C}{n + p} $$ - For fully ionized n-type dopant at high concentration: $h \approx 2$ - Results in approximately 2× faster effective diffusion ### 4.4 Resulting Profile Shapes - **Phosphorus:** "Kink-and-tail" profile at high concentrations - **Arsenic:** Box-like profiles due to clustering - **Boron:** Enhanced tail diffusion in oxidizing ambient ## 5. Point Defect-Mediated Diffusion ### 5.1 Diffusion Mechanisms Dopants don't diffuse as isolated atoms—they move via **defect complexes**: #### Vacancy Mechanism $$ A + V \rightleftharpoons AV \quad \text{(dopant-vacancy pair forms, diffuses, dissociates)} $$ #### Interstitial Mechanism $$ A + I \rightleftharpoons AI \quad \text{(dopant-interstitial pair)} $$ #### Kick-out Mechanism $$ A_s + I \rightleftharpoons A_i \quad \text{(substitutional ↔ interstitial)} $$ ### 5.2 Effective Diffusivity $$ D_{eff} = D_V \frac{C_V}{C_V^*} + D_I \frac{C_I}{C_I^*} $$ **Where:** - $D_V, D_I$ = diffusivity via vacancy/interstitial mechanism - $C_V, C_I$ = actual vacancy/interstitial concentrations - $C_V^*, C_I^*$ = equilibrium concentrations **Fractional interstitialcy:** $$ f_I = \frac{D_I}{D_V + D_I} $$ | Dopant | $f_I$ | Dominant Mechanism | |--------|-------|-------------------| | Boron | ~1.0 | Interstitial | | Phosphorus | ~0.9 | Interstitial | | Arsenic | ~0.4 | Mixed | | Antimony | ~0.02 | Vacancy | ### 5.3 Coupled Reaction-Diffusion System The full model requires solving **coupled PDEs**: #### Dopant Equation $$ \frac{\partial C_A}{\partial t} = \nabla \cdot \left(D_A \frac{C_I}{C_I^*} \nabla C_A\right) $$ #### Interstitial Balance $$ \frac{\partial C_I}{\partial t} = D_I \nabla^2 C_I + G - k_{IV}\left(C_I C_V - C_I^* C_V^*\right) $$ #### Vacancy Balance $$ \frac{\partial C_V}{\partial t} = D_V \nabla^2 C_V + G - k_{IV}\left(C_I C_V - C_I^* C_V^*\right) $$ **Where:** - $G$ = defect generation rate - $k_{IV}$ = bulk recombination rate constant ### 5.4 Transient Enhanced Diffusion (TED) After ion implantation, excess interstitials cause **anomalously rapid diffusion**: **The "+1" Model:** $$ \int_0^\infty (C_I - C_I^*) \, dx \approx \Phi \quad \text{(implant dose)} $$ **Enhancement factor:** $$ \frac{D_{eff}}{D^*} = \frac{C_I}{C_I^*} \gg 1 \quad \text{(transient)} $$ **Key characteristics:** - Enhancement decays as interstitials recombine - Time constant: typically 10-100 seconds at 1000°C - Critical for shallow junction formation ## 6. Oxidation Effects ### 6.1 Oxidation-Enhanced Diffusion (OED) During thermal oxidation, silicon interstitials are **injected** into the substrate: $$ \frac{C_I}{C_I^*} = 1 + A \left(\frac{dx_{ox}}{dt}\right)^n $$ **Effective diffusivity:** $$ D_{eff} = D^* \left[1 + f_I \left(\frac{C_I}{C_I^*} - 1\right)\right] $$ **Dopants enhanced by oxidation:** - Boron (high $f_I$) - Phosphorus (high $f_I$) ### 6.2 Oxidation-Retarded Diffusion (ORD) Growing oxide **absorbs vacancies**, reducing vacancy concentration: $$ \frac{C_V}{C_V^*} < 1 $$ **Dopants retarded by oxidation:** - Antimony (low $f_I$, primarily vacancy-mediated) ### 6.3 Segregation at SiO₂/Si Interface Dopants redistribute at the interface according to the **segregation coefficient**: $$ m = \frac{C_{Si}}{C_{SiO_2}}\bigg|_{\text{interface}} $$ | Dopant | Segregation Coefficient $m$ | Behavior | |--------|----------------------------|----------| | Boron | ~0.3 | Pile-down (into oxide) | | Phosphorus | ~10 | Pile-up (into silicon) | | Arsenic | ~10 | Pile-up | ## 7. Numerical Methods ### 7.1 Finite Difference Method Discretize space and time on grid $(x_i, t^n)$: #### Explicit Scheme (FTCS) $$ \frac{C_i^{n+1} - C_i^n}{\Delta t} = D \frac{C_{i+1}^n - 2C_i^n + C_{i-1}^n}{(\Delta x)^2} $$ **Rearranged:** $$ C_i^{n+1} = C_i^n + \alpha \left(C_{i+1}^n - 2C_i^n + C_{i-1}^n\right) $$ **Where Fourier number:** $$ \alpha = \frac{D \Delta t}{(\Delta x)^2} $$ **Stability requirement (von Neumann analysis):** $$ \alpha \leq \frac{1}{2} $$ #### Implicit Scheme (BTCS) $$ \frac{C_i^{n+1} - C_i^n}{\Delta t} = D \frac{C_{i+1}^{n+1} - 2C_i^{n+1} + C_{i-1}^{n+1}}{(\Delta x)^2} $$ - **Unconditionally stable** (no restriction on $\alpha$) - Requires solving tridiagonal system at each time step #### Crank-Nicolson Scheme (Second-Order Accurate) $$ C_i^{n+1} - C_i^n = \frac{\alpha}{2}\left[(C_{i+1}^{n+1} - 2C_i^{n+1} + C_{i-1}^{n+1}) + (C_{i+1}^n - 2C_i^n + C_{i-1}^n)\right] $$ **Properties:** - Unconditionally stable - Second-order accurate in both space and time - Results in tridiagonal system: solved by **Thomas algorithm** ### 7.2 Handling Concentration-Dependent Diffusion Use iterative methods: 1. Estimate $D^{(k)}$ from current concentration $C^{(k)}$ 2. Solve linear diffusion equation for $C^{(k+1)}$ 3. Update diffusivity: $D^{(k+1)} = D(C^{(k+1)})$ 4. Iterate until $\|C^{(k+1)} - C^{(k)}\| < \epsilon$ ### 7.3 Moving Boundary Problems For oxidation with moving Si/SiO₂ interface: **Approaches:** - **Coordinate transformation:** Map to fixed domain via $\xi = x/s(t)$ - **Front-tracking methods:** Explicitly track interface position - **Level-set methods:** Implicit interface representation - **Phase-field methods:** Diffuse interface approximation ## 8. Thermal Budget Concept ### 8.1 The Dt Product Diffusion profiles scale with $\sqrt{Dt}$. The **thermal budget** quantifies total diffusion: $$ (Dt)_{total} = \sum_i D(T_i) \cdot t_i $$ ### 8.2 Continuous Temperature Profile For time-varying temperature: $$ (Dt)_{eff} = \int_0^{t_{total}} D(T(\tau)) \, d\tau $$ ### 8.3 Equivalent Time at Reference Temperature $$ t_{eq} = \sum_i t_i \exp\left(\frac{E_a}{k}\left(\frac{1}{T_{ref}} - \frac{1}{T_i}\right)\right) $$ ### 8.4 Combining Multiple Diffusion Steps For sequential Gaussian redistributions: $$ \sigma_{final} = \sqrt{\sum_i 2D_i t_i} $$ For erfc profiles, use effective $(Dt)_{total}$: $$ C(x) = C_s \cdot \text{erfc}\left(\frac{x}{2\sqrt{(Dt)_{total}}}\right) $$ ## 9. Key Dimensionless Parameters | Parameter | Definition | Physical Meaning | |-----------|------------|------------------| | **Fourier Number** | $Fo = \dfrac{Dt}{L^2}$ | Diffusion time vs. characteristic length | | **Damköhler Number** | $Da = \dfrac{kL^2}{D}$ | Reaction rate vs. diffusion rate | | **Péclet Number** | $Pe = \dfrac{vL}{D}$ | Advection (drift) vs. diffusion | | **Biot Number** | $Bi = \dfrac{hL}{D}$ | Surface transfer vs. bulk diffusion | ## 10. Process Simulation Software ### 10.1 Commercial and Research Tools | Simulator | Developer | Key Capabilities | |-----------|-----------|------------------| | **Sentaurus Process** | Synopsys | Full 3D, atomistic KMC, advanced models | | **Athena** | Silvaco | Integrated with device simulation (Atlas) | | **SUPREM-IV** | Stanford | Classic 1D/2D, widely validated | | **FLOOPS** | U. Florida | Research-oriented, extensible | | **Victory Process** | Silvaco | Modern 3D process simulation | ### 10.2 Physical Models Incorporated - Multiple coupled dopant species - Full point-defect dynamics (I, V, clusters) - Stress-dependent diffusion - Cluster nucleation and dissolution - Atomistic kinetic Monte Carlo (KMC) options - Quantum corrections for ultra-shallow junctions ## Mathematical Modeling Hierarchy ### Level 1: Simple Analytical Models $$ \frac{\partial C}{\partial t} = D \frac{\partial^2 C}{\partial x^2} $$ - Constant $D$ - erfc and Gaussian solutions - Junction depth calculations ### Level 2: Intermediate Complexity $$ \frac{\partial C}{\partial t} = \frac{\partial}{\partial x}\left(D(C) \frac{\partial C}{\partial x}\right) $$ - Concentration-dependent $D$ - Electric field effects - Nonlinear PDEs requiring numerical methods ### Level 3: Advanced Coupled Models $$ \begin{aligned} \frac{\partial C_A}{\partial t} &= \nabla \cdot \left(D_A \frac{C_I}{C_I^*} \nabla C_A\right) \\[6pt] \frac{\partial C_I}{\partial t} &= D_I \nabla^2 C_I + G - k_{IV}(C_I C_V - C_I^* C_V^*) \end{aligned} $$ - Coupled dopant-defect systems - TED, OED/ORD effects - Process simulators required ### Level 4: State-of-the-Art - Atomistic kinetic Monte Carlo - Molecular dynamics for interface phenomena - Ab initio calculations for defect properties - Essential for sub-10nm technology nodes ## Key Insight The fundamental scaling of semiconductor diffusion is governed by $\sqrt{Dt}$, but the effective diffusion coefficient $D$ depends on: - Temperature (Arrhenius) - Concentration (charged defects) - Point defect supersaturation (TED) - Processing ambient (oxidation) - Mechanical stress This complexity requires sophisticated physical models for modern nanometer-scale devices.
Virtual replica for simulation.
Allowed variation in dimensions.
Bond wafers without intermediate layer.
Ability to detect small differences.
# Design of Experiments (DOE) in Semiconductor Manufacturing DOE is a statistical methodology for systematically investigating relationships between process parameters and responses (yield, thickness, defects, etc.). 1. Fundamental Mathematical Model First-order linear model: y = β₀ + Σᵢβᵢxᵢ + ε Second-order model (with curvature and interactions): y = β₀ + Σᵢβᵢxᵢ + Σᵢβᵢᵢxᵢ² + Σᵢ<ⱼβᵢⱼxᵢxⱼ + ε Where: • y = response (oxide thickness, threshold voltage) • xᵢ = coded factor levels (scaled to [-1, +1]) • β = model coefficients • ε = random error ~ N(0, σ²) 2. Matrix Formulation Model in matrix form: Y = Xβ + ε Least squares estimation: β̂ = (X'X)⁻¹X'Y Variance-covariance of estimates: Var(β̂) = σ²(X'X)⁻¹ 3. Factorial Designs Full Factorial (2ᵏ) For k factors at 2 levels: requires 2ᵏ runs. Orthogonality property: X'X = nI All effects estimated independently with equal precision. Fractional Factorial (2ᵏ⁻ᵖ) Resolution determines confounding: • Resolution III: Main effects aliased with 2FIs • Resolution IV: Main effects clear; 2FIs aliased with each other • Resolution V: Main effects and 2FIs all estimable For 2⁵⁻² design with generators D = AB, E = AC: • Defining relation: I = ABD = ACE = BCDE • Find aliases by multiplying effect by defining relation 4. Response Surface Methodology (RSM) Central Composite Design (CCD) Combines: • 2ᵏ or 2ᵏ⁻ᵖ factorial points • 2k axial points at ±α from center • n₀ center points Rotatability condition: α = (2ᵏ)¹/⁴ = F¹/⁴ • For k=2: α = √2 ≈ 1.414 • For k=3: α = 2³/⁴ ≈ 1.682 Box-Behnken Design • 3 levels per factor • No corner points (useful when extremes are dangerous) • More economical than CCD for 3+ factors 5. Optimal Design Theory D-optimal: Maximize |X'X| • Minimizes volume of joint confidence region A-optimal: Minimize trace[(X'X)⁻¹] • Minimizes average variance of estimates I-optimal: Minimize integrated prediction variance: ∫ Var[ŷ(x)] dx G-optimal: Minimize maximum prediction variance 6. Analysis of Variance (ANOVA) Sum of squares decomposition: SSₜₒₜₐₗ = SSₘₒdₑₗ + SSᵣₑₛᵢdᵤₐₗ SSₘₒdₑₗ = Σᵢ(ŷᵢ - ȳ)² SSᵣₑₛᵢdᵤₐₗ = Σᵢ(yᵢ - ŷᵢ)² F-test for significance: F = MSₑffₑcₜ / MSₑᵣᵣₒᵣ = (SSₑffₑcₜ/dfₑffₑcₜ) / (SSₑᵣᵣₒᵣ/dfₑᵣᵣₒᵣ) Effect estimation: Effectₐ = ȳₐ₊ - ȳₐ₋ β̂ₐ = Effectₐ / 2 7. Semiconductor-Specific Designs Split-Plot Designs For hard-to-change factors (temperature, pressure) vs easy-to-change (gas flow): yᵢⱼₖ = μ + αᵢ + δᵢⱼ + βₖ + (αβ)ᵢₖ + εᵢⱼₖ Where: • αᵢ = whole-plot factor (hard to change) • δᵢⱼ = whole-plot error • βₖ = subplot factor (easy to change) • εᵢⱼₖ = subplot error Variance Components (Nested Designs) For Lots → Wafers → Dies → Measurements: σ²ₜₒₜₐₗ = σ²ₗₒₜ + σ²wₐfₑᵣ + σ²dᵢₑ + σ²ₘₑₐₛ Mixture Designs For etch gas chemistry where components sum to 1: Σᵢxᵢ = 1 Uses simplex-lattice designs and Scheffé models. 8. Robust Parameter Design (Taguchi) Signal-to-Noise ratios: Nominal-is-best: S/N = 10·log₁₀(ȳ²/s²) Smaller-is-better: S/N = -10·log₁₀[(1/n)·Σyᵢ²] Larger-is-better: S/N = -10·log₁₀[(1/n)·Σ(1/yᵢ²)] 9. Sequential Optimization Steepest Ascent/Descent: ∇y = (β₁, β₂, ..., βₖ) Step sizes: Δxᵢ ∝ βᵢ × (range of xᵢ) 10. Model Diagnostics Coefficient of determination: R² = 1 - SSᵣₑₛᵢdᵤₐₗ/SSₜₒₜₐₗ Adjusted R²: R²ₐdⱼ = 1 - [SSᵣₑₛᵢdᵤₐₗ/(n-p)] / [SSₜₒₜₐₗ/(n-1)] PRESS statistic: PRESS = Σᵢ(yᵢ - ŷ₍ᵢ₎)² Prediction R²: R²ₚᵣₑd = 1 - PRESS/SSₜₒₜₐₗ Variance Inflation Factor: VIFⱼ = 1/(1 - R²ⱼ) VIF > 10 indicates problematic collinearity. 11. Power and Sample Size Minimum detectable effect: δ = σ × √[2(zₐ/₂ + zᵦ)²/n] Power calculation: Power = Φ(|δ|√n / (σ√2) - zₐ/₂) 12. Multivariate Optimization Desirability function for target T between L and U: d = [(y-L)/(T-L)]ˢ when L ≤ y ≤ T d = [(U-y)/(U-T)]ᵗ when T ≤ y ≤ U Overall desirability: D = (∏ᵢdᵢʷⁱ)^(1/Σwᵢ) 13. Process Capability Integration Cₚ = (USL - LSL) / 6σ Cₚₖ = min[(USL - μ)/3σ, (μ - LSL)/3σ] DOE improves Cₚₖ by centering and reducing variation. 14. Model Selection AIC: AIC = n·ln(SSE/n) + 2p BIC: BIC = n·ln(SSE/n) + p·ln(n) 15. Modern Advances Definitive Screening Designs (DSD) • Jones & Nachtsheim (2011) • Requires only 2k+1 runs for k factors • Estimates main effects, quadratic effects, and some 2FIs Bayesian DOE • Prior: p(β) • Posterior: p(β|Y) ∝ p(Y|β)p(β) • Expected Improvement for sequential selection Gaussian Process (Kriging) • Non-parametric, data-driven • Provides uncertainty quantification Summary DOE provides the rigorous framework for process optimization where: • Single experiments cost tens of thousands of dollars • Cycle times span weeks to months • Maximum information from minimum runs is essential
Tests within die area.
Packaging moisture controls.
Solid resist film instead of liquid.