die, dies, dicing, singulation, yield, wafer dicing, die separation, semiconductor dicing, wafer singulation
# 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 |
dielectric capping layer,beol
Protect low-k during processing.
dielectric loss, signal & power integrity
Dielectric loss dissipates signal energy in insulator materials increasing with frequency.
diff-gan graph, graph neural networks
Diffusion-GAN combines diffusion models with adversarial training for diverse graph generation.
differentiable architecture search, darts, neural architecture
Search architectures via gradient descent.
differentiable mpc, control theory
MPC with differentiable dynamics for end-to-end learning.
differentiable neural computer (dnc),differentiable neural computer,dnc,neural architecture
Advanced memory-augmented network.
differentiable physics engines, physics simulation
Physics simulators with gradients.
differentiable programming,programming
Programs where components are differentiable.
differentiable rasterization, 3d vision
Render Gaussians with gradients.
differentiable rendering, multimodal ai
Differentiable rendering enables gradient-based optimization of 3D representations from images.
differentiable rendering,computer vision
Rendering with gradients for optimization.
differential impedance, signal & power integrity
Differential impedance is the ratio of differential voltage to differential current in paired transmission lines.
differential phase contrast, dpc, metrology
Map electric fields in materials.
differential privacy in federated learning, federated learning
Add DP noise to protect privacy.
differential privacy rec, recommendation systems
Differential privacy in recommendations adds calibrated noise to protect individual user information during training.
differential privacy, training techniques
Differential privacy adds calibrated noise protecting individual data points during training.
differential privacy,ai safety
Add noise during training to prevent extracting individual training examples.
differential privacy,dp,noise
Differential privacy adds noise to protect individual data. Epsilon measures privacy level.
differential privacy,dp,noise
Differential privacy adds noise to prevent identifying individuals in training data. Formal privacy guarantee.
differential signaling, signal & power integrity
Differential signaling transmits signals as voltage difference between complementary lines rejecting common-mode noise.
differential signaling,design
Use signal pairs for noise immunity.
differential testing,software testing
Compare outputs across implementations.
diffpool, graph neural networks
Differentiable graph pooling learns hierarchical graph representations by clustering nodes through soft assignment matrices optimized end-to-end.
diffpool, graph neural networks
Differentiable graph pooling.
diffraction-based overlay, dbo, metrology
Measure overlay using diffraction gratings.
diffusers,huggingface,stable diffusion
Diffusers library for diffusion models. Stable Diffusion, ControlNet. Hugging Face.
diffusion and ion implantation,diffusion,ion implantation,dopant diffusion,fick law,implant profile,gaussian profile,pearson distribution,ted,transient enhanced diffusion,thermal budget,semiconductor doping
# 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 |
diffusion bonding, business & strategy
Diffusion bonding joins surfaces through atomic interdiffusion at interfaces.
diffusion coefficient,diffusion
Rate at which dopants diffuse depends on temp and material.
diffusion equations,fick laws,fick second law,semiconductor diffusion equations,dopant diffusion equations,arrhenius diffusion,junction depth calculation,transient enhanced diffusion,oxidation enhanced diffusion,numerical methods diffusion,thermal budget
# 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.
diffusion furnace,diffusion
Tube furnace for thermal processing of wafers.
diffusion language models, generative models
Apply diffusion process to discrete text generation.
diffusion length,lithography
How far acid or other species diffuse in resist.
diffusion model training, generative models
Train denoising models.
diffusion modeling, diffusion model, fick law modeling, dopant diffusion model, semiconductor diffusion model, thermal diffusion model, diffusion coefficient calculation, diffusion simulation, diffusion mathematics
# 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.
diffusion models for graphs, graph neural networks
Apply diffusion to graph generation.
diffusion models,generative models
Generative models that gradually denoise random noise into realistic images.
diffusion on graphs, graph neural networks
Model information diffusion on graph structure.
diffusion simulation, simulation
Model dopant movement at high temperature.
diffusion upscaler, multimodal ai
Diffusion-based upscalers generate high-resolution details through iterative refinement.
diffusion-lm, foundation model
Language model using continuous diffusion.
diffusion,denoising,generative
Diffusion models denoise to generate. Start from noise, iteratively refine. DALL-E, Stable Diffusion.
diffusion,stable diffusion,image gen
Diffusion models generate images from text by iteratively denoising. Stable Diffusion is open-source and customizable.
dify,llmops,platform
Dify is LLMOps platform. Prompt engineering, RAG, agents.
digital twin for robotics,robotics
Virtual replica for testing and training.
digital twin of semiconductor fab, digital manufacturing
Virtual replica for simulation.
digital twin,production
Virtual model of tool or process for simulation and optimization.
dilated attention in vision, computer vision
Attend to patches at intervals.
dilated attention,llm architecture
Attend to tokens at regular intervals for long-range dependencies.