rie
**Mathematical Modeling of Plasma Etching in Semiconductor Manufacturing**
**Introduction**
Plasma etching is a critical process in semiconductor manufacturing where reactive gases are ionized to create a plasma, which selectively removes material from a wafer surface. The mathematical modeling of this process spans multiple physics domains:
- **Electromagnetic theory** — RF power coupling and field distributions
- **Statistical mechanics** — Particle distributions and kinetic theory
- **Reaction kinetics** — Gas-phase and surface chemistry
- **Transport phenomena** — Species diffusion and convection
- **Surface science** — Etch mechanisms and selectivity
**Foundational Plasma Physics**
**Boltzmann Transport Equation**
The most fundamental description of plasma behavior is the **Boltzmann transport equation**, governing the evolution of the particle velocity distribution function $f(\mathbf{r}, \mathbf{v}, t)$:
$$
\frac{\partial f}{\partial t} + \mathbf{v} \cdot
abla f + \frac{\mathbf{F}}{m} \cdot
abla_v f = \left(\frac{\partial f}{\partial t}\right)_{\text{collision}}
$$
**Where:**
- $f(\mathbf{r}, \mathbf{v}, t)$ — Velocity distribution function
- $\mathbf{v}$ — Particle velocity
- $\mathbf{F}$ — External force (electromagnetic)
- $m$ — Particle mass
- RHS — Collision integral
**Fluid Moment Equations**
For computational tractability, velocity moments of the Boltzmann equation yield fluid equations:
**Continuity Equation (Mass Conservation)**
$$
\frac{\partial n}{\partial t} +
abla \cdot (n\mathbf{u}) = S - L
$$
**Where:**
- $n$ — Species number density $[\text{m}^{-3}]$
- $\mathbf{u}$ — Drift velocity $[\text{m/s}]$
- $S$ — Source term (generation rate)
- $L$ — Loss term (consumption rate)
**Momentum Conservation**
$$
\frac{\partial (nm\mathbf{u})}{\partial t} +
abla \cdot (nm\mathbf{u}\mathbf{u}) +
abla p = nq(\mathbf{E} + \mathbf{u} \times \mathbf{B}) - nm
u_m \mathbf{u}
$$
**Where:**
- $p = nk_BT$ — Pressure
- $q$ — Particle charge
- $\mathbf{E}$, $\mathbf{B}$ — Electric and magnetic fields
- $
u_m$ — Momentum transfer collision frequency $[\text{s}^{-1}]$
**Energy Conservation**
$$
\frac{\partial}{\partial t}\left(\frac{3}{2}nk_BT\right) +
abla \cdot \mathbf{q} + p
abla \cdot \mathbf{u} = Q_{\text{heating}} - Q_{\text{loss}}
$$
**Where:**
- $k_B = 1.38 \times 10^{-23}$ J/K — Boltzmann constant
- $\mathbf{q}$ — Heat flux vector
- $Q_{\text{heating}}$ — Power input (Joule heating, stochastic heating)
- $Q_{\text{loss}}$ — Energy losses (collisions, radiation)
**Electromagnetic Field Coupling**
**Maxwell's Equations**
For capacitively coupled plasma (CCP) and inductively coupled plasma (ICP) reactors:
$$
abla \times \mathbf{E} = -\frac{\partial \mathbf{B}}{\partial t}
$$
$$
abla \times \mathbf{H} = \mathbf{J} + \frac{\partial \mathbf{D}}{\partial t}
$$
$$
abla \cdot \mathbf{D} = \rho
$$
$$
abla \cdot \mathbf{B} = 0
$$
**Plasma Conductivity**
The plasma current density couples through the complex conductivity:
$$
\mathbf{J} = \sigma \mathbf{E}
$$
For RF plasmas, the **complex conductivity** is:
$$
\sigma = \frac{n_e e^2}{m_e(
u_m + i\omega)}
$$
**Where:**
- $n_e$ — Electron density
- $e = 1.6 \times 10^{-19}$ C — Elementary charge
- $m_e = 9.1 \times 10^{-31}$ kg — Electron mass
- $\omega$ — RF angular frequency
- $
u_m$ — Electron-neutral collision frequency
**Power Deposition**
Time-averaged power density deposited into the plasma:
$$
P = \frac{1}{2}\text{Re}(\mathbf{J} \cdot \mathbf{E}^*)
$$
**Typical values:**
- CCP: $0.1 - 1$ W/cm³
- ICP: $0.5 - 5$ W/cm³
**Plasma Sheath Physics**
The sheath is a thin, non-neutral region at the plasma-wafer interface that accelerates ions toward the surface, enabling anisotropic etching.
**Bohm Criterion**
Minimum ion velocity entering the sheath:
$$
u_i \geq u_B = \sqrt{\frac{k_B T_e}{M_i}}
$$
**Where:**
- $u_B$ — Bohm velocity
- $T_e$ — Electron temperature (typically 2–5 eV)
- $M_i$ — Ion mass
**Example:** For Ar⁺ ions with $T_e = 3$ eV:
$$
u_B = \sqrt{\frac{3 \times 1.6 \times 10^{-19}}{40 \times 1.67 \times 10^{-27}}} \approx 2.7 \text{ km/s}
$$
**Child-Langmuir Law**
For a collisionless sheath, the ion current density is:
$$
J = \frac{4\varepsilon_0}{9}\sqrt{\frac{2e}{M_i}} \cdot \frac{V_s^{3/2}}{d^2}
$$
**Where:**
- $\varepsilon_0 = 8.85 \times 10^{-12}$ F/m — Vacuum permittivity
- $V_s$ — Sheath voltage drop (typically 10–500 V)
- $d$ — Sheath thickness
**Sheath Thickness**
The sheath thickness scales as:
$$
d \approx \lambda_D \left(\frac{2eV_s}{k_BT_e}\right)^{3/4}
$$
**Where** the Debye length is:
$$
\lambda_D = \sqrt{\frac{\varepsilon_0 k_B T_e}{n_e e^2}}
$$
**Ion Angular Distribution**
Ions arrive at the wafer with an angular distribution:
$$
f(\theta) \propto \exp\left(-\frac{\theta^2}{2\sigma^2}\right)
$$
**Where:**
$$
\sigma \approx \arctan\left(\sqrt{\frac{k_B T_i}{eV_s}}\right)
$$
**Typical values:** $\sigma \approx 2°–5°$ for high-bias conditions.
**Electron Energy Distribution Function**
**Non-Maxwellian Distributions**
In low-pressure plasmas (1–100 mTorr), the EEDF deviates from Maxwellian.
**Two-Term Approximation**
The EEDF is expanded as:
$$
f(\varepsilon, \theta) = f_0(\varepsilon) + f_1(\varepsilon)\cos\theta
$$
The isotropic part $f_0$ satisfies:
$$
\frac{d}{d\varepsilon}\left[\varepsilon D \frac{df_0}{d\varepsilon} + \left(V + \frac{\varepsilon
u_{\text{inel}}}{
u_m}\right)f_0\right] = 0
$$
**Common Distribution Functions**
| Distribution | Functional Form | Applicability |
|-------------|-----------------|---------------|
| **Maxwellian** | $f(\varepsilon) \propto \sqrt{\varepsilon} \exp\left(-\frac{\varepsilon}{k_BT_e}\right)$ | High pressure, collisional |
| **Druyvesteyn** | $f(\varepsilon) \propto \sqrt{\varepsilon} \exp\left(-\left(\frac{\varepsilon}{k_BT_e}\right)^2\right)$ | Elastic collisions dominant |
| **Bi-Maxwellian** | Sum of two Maxwellians | Hot tail population |
**Generalized Form**
$$
f(\varepsilon) \propto \sqrt{\varepsilon} \cdot \exp\left[-\left(\frac{\varepsilon}{k_BT_e}\right)^x\right]
$$
- $x = 1$ → Maxwellian
- $x = 2$ → Druyvesteyn
**Plasma Chemistry and Reaction Kinetics**
**Species Balance Equation**
For species $i$:
$$
\frac{\partial n_i}{\partial t} +
abla \cdot \mathbf{\Gamma}_i = \sum_j R_j
$$
**Where:**
- $\mathbf{\Gamma}_i$ — Species flux
- $R_j$ — Reaction rates
**Electron-Impact Rate Coefficients**
Rate coefficients are calculated by integration over the EEDF:
$$
k = \int_0^\infty \sigma(\varepsilon) v(\varepsilon) f(\varepsilon) \, d\varepsilon = \langle \sigma v \rangle
$$
**Where:**
- $\sigma(\varepsilon)$ — Energy-dependent cross-section $[\text{m}^2]$
- $v(\varepsilon) = \sqrt{2\varepsilon/m_e}$ — Electron velocity
- $f(\varepsilon)$ — Normalized EEDF
**Heavy-Particle Reactions**
Arrhenius kinetics for neutral reactions:
$$
k = A T^n \exp\left(-\frac{E_a}{k_BT}\right)
$$
**Where:**
- $A$ — Pre-exponential factor
- $n$ — Temperature exponent
- $E_a$ — Activation energy
**Example: SF₆/O₂ Plasma Chemistry**
**Electron-Impact Reactions**
| Reaction | Type | Threshold |
|----------|------|-----------|
| $e + \text{SF}_6 \rightarrow \text{SF}_5 + \text{F} + e$ | Dissociation | ~10 eV |
| $e + \text{SF}_6 \rightarrow \text{SF}_6^-$ | Attachment | ~0 eV |
| $e + \text{SF}_6 \rightarrow \text{SF}_5^+ + \text{F} + 2e$ | Ionization | ~16 eV |
| $e + \text{O}_2 \rightarrow \text{O} + \text{O} + e$ | Dissociation | ~6 eV |
**Gas-Phase Reactions**
- $\text{F} + \text{O} \rightarrow \text{FO}$ (reduces F atom density)
- $\text{SF}_5 + \text{F} \rightarrow \text{SF}_6$ (recombination)
- $\text{O} + \text{CF}_3 \rightarrow \text{COF}_2 + \text{F}$ (polymer removal)
**Surface Reactions**
- $\text{F} + \text{Si}(s) \rightarrow \text{SiF}_{(\text{ads})}$
- $\text{SiF}_{(\text{ads})} + 3\text{F} \rightarrow \text{SiF}_4(g)$ (volatile product)
**Transport Phenomena**
**Drift-Diffusion Model**
For charged species, the flux is:
$$
\mathbf{\Gamma} = \pm \mu n \mathbf{E} - D
abla n
$$
**Where:**
- Upper sign: positive ions
- Lower sign: electrons
- $\mu$ — Mobility $[\text{m}^2/(\text{V}\cdot\text{s})]$
- $D$ — Diffusion coefficient $[\text{m}^2/\text{s}]$
**Einstein Relation**
Connects mobility and diffusion:
$$
D = \frac{\mu k_B T}{e}
$$
**Ambipolar Diffusion**
When quasi-neutrality holds ($n_e \approx n_i$):
$$
D_a = \frac{\mu_i D_e + \mu_e D_i}{\mu_i + \mu_e} \approx D_i\left(1 + \frac{T_e}{T_i}\right)
$$
Since $T_e \gg T_i$ typically: $D_a \approx D_i (1 + T_e/T_i) \approx 100 D_i$
**Neutral Transport**
For reactive neutrals (radicals), Fickian diffusion:
$$
\frac{\partial n}{\partial t} = D
abla^2 n + S - L
$$
**Surface Boundary Condition**
$$
-D\frac{\partial n}{\partial x}\bigg|_{\text{surface}} = \frac{1}{4}\gamma n v_{\text{th}}
$$
**Where:**
- $\gamma$ — Sticking/reaction coefficient (0 to 1)
- $v_{\text{th}} = \sqrt{\frac{8k_BT}{\pi m}}$ — Thermal velocity
**Knudsen Number**
Determines the appropriate transport regime:
$$
\text{Kn} = \frac{\lambda}{L}
$$
**Where:**
- $\lambda$ — Mean free path
- $L$ — Characteristic length
| Kn Range | Regime | Model |
|----------|--------|-------|
| $< 0.01$ | Continuum | Navier-Stokes |
| $0.01–0.1$ | Slip flow | Modified N-S |
| $0.1–10$ | Transition | DSMC/BGK |
| $> 10$ | Free molecular | Ballistic |
**Surface Reaction Modeling**
**Langmuir Adsorption Kinetics**
For surface coverage $\theta$:
$$
\frac{d\theta}{dt} = k_{\text{ads}}(1-\theta)P - k_{\text{des}}\theta - k_{\text{react}}\theta
$$
**At steady state:**
$$
\theta = \frac{k_{\text{ads}}P}{k_{\text{ads}}P + k_{\text{des}} + k_{\text{react}}}
$$
**Ion-Enhanced Etching**
The total etch rate combines multiple mechanisms:
$$
\text{ER} = Y_{\text{chem}} \Gamma_n + Y_{\text{phys}} \Gamma_i + Y_{\text{syn}} \Gamma_i f(\theta)
$$
**Where:**
- $Y_{\text{chem}}$ — Chemical etch yield (isotropic)
- $Y_{\text{phys}}$ — Physical sputtering yield
- $Y_{\text{syn}}$ — Ion-enhanced (synergistic) yield
- $\Gamma_n$, $\Gamma_i$ — Neutral and ion fluxes
- $f(\theta)$ — Coverage-dependent function
**Ion Sputtering Yield**
**Energy Dependence**
$$
Y(E) = A\left(\sqrt{E} - \sqrt{E_{\text{th}}}\right) \quad \text{for } E > E_{\text{th}}
$$
**Typical threshold energies:**
- Si: $E_{\text{th}} \approx 20$ eV
- SiO₂: $E_{\text{th}} \approx 30$ eV
- Si₃N₄: $E_{\text{th}} \approx 25$ eV
**Angular Dependence**
$$
Y(\theta) = Y(0) \cos^{-f}(\theta) \exp\left[-b\left(\frac{1}{\cos\theta} - 1\right)\right]
$$
**Behavior:**
- Increases from normal incidence
- Peaks at $\theta \approx 60°–70°$
- Decreases at grazing angles (reflection dominates)
**Feature-Scale Profile Evolution**
**Level Set Method**
The surface is represented as the zero contour of $\phi(\mathbf{x}, t)$:
$$
\frac{\partial \phi}{\partial t} + V_n |
abla \phi| = 0
$$
**Where:**
- $\phi > 0$ — Material
- $\phi < 0$ — Void/vacuum
- $\phi = 0$ — Surface
- $V_n$ — Local normal etch velocity
**Local Etch Rate Calculation**
The normal velocity $V_n$ depends on:
1. **Ion flux and angular distribution**
$$\Gamma_i(\mathbf{x}) = \int f(\theta, E) \, d\Omega \, dE$$
2. **Neutral flux** (with shadowing)
$$\Gamma_n(\mathbf{x}) = \Gamma_{n,0} \cdot \text{VF}(\mathbf{x})$$
where VF is the view factor
3. **Surface chemistry state**
$$V_n = f(\Gamma_i, \Gamma_n, \theta_{\text{coverage}}, T)$$
**Neutral Transport in High-Aspect-Ratio Features**
**Clausing Transmission Factor**
For a tube of aspect ratio AR:
$$
K \approx \frac{1}{1 + 0.5 \cdot \text{AR}}
$$
**View Factor Calculations**
For surface element $dA_1$ seeing $dA_2$:
$$
F_{1 \rightarrow 2} = \frac{1}{\pi} \int \frac{\cos\theta_1 \cos\theta_2}{r^2} \, dA_2
$$
**Monte Carlo Methods**
**Test-Particle Monte Carlo Algorithm**
```
1. SAMPLE incident particle from flux distribution at feature opening
- Ion: from IEDF and IADF
- Neutral: from Maxwellian
2. TRACE trajectory through feature
- Ion: ballistic, solve equation of motion
- Neutral: random walk with wall collisions
3. DETERMINE reaction at surface impact
- Sample from probability distribution
- Update surface coverage if adsorption
4. UPDATE surface geometry
- Remove material (etching)
- Add material (deposition)
5. REPEAT for statistically significant sample
```
**Ion Trajectory Integration**
Through the sheath/feature:
$$
m\frac{d^2\mathbf{r}}{dt^2} = q\mathbf{E}(\mathbf{r})
$$
**Numerical integration:** Velocity-Verlet or Boris algorithm
**Collision Sampling**
Null-collision method for efficiency:
$$
P_{\text{collision}} = 1 - \exp(-
u_{\text{max}} \Delta t)
$$
**Where** $
u_{\text{max}}$ is the maximum possible collision frequency.
**Multi-Scale Modeling Framework**
**Scale Hierarchy**
| Scale | Length | Time | Physics | Method |
|-------|--------|------|---------|--------|
| **Reactor** | cm–m | ms–s | Plasma transport, EM fields | Fluid PDE |
| **Sheath** | µm–mm | µs–ms | Ion acceleration, EEDF | Kinetic/Fluid |
| **Feature** | nm–µm | ns–ms | Profile evolution | Level set/MC |
| **Atomic** | Å–nm | ps–ns | Reaction mechanisms | MD/DFT |
**Coupling Approaches**
**Hierarchical (One-Way)**
```
Atomic scale → Surface parameters
↓
Feature scale ← Fluxes from reactor scale
↓
Reactor scale → Process outputs
```
**Concurrent (Two-Way)**
- Feature-scale results feed back to reactor scale
- Requires iterative solution
- Computationally expensive
**Numerical Methods and Challenges**
**Stiff ODE Systems**
Plasma chemistry involves timescales spanning many orders of magnitude:
| Process | Timescale |
|---------|-----------|
| Electron attachment | $\sim 10^{-10}$ s |
| Ion-molecule reactions | $\sim 10^{-6}$ s |
| Metastable decay | $\sim 10^{-3}$ s |
| Surface diffusion | $\sim 10^{-1}$ s |
**Implicit Methods Required**
**Backward Differentiation Formula (BDF):**
$$
y_{n+1} = \sum_{j=0}^{k-1} \alpha_j y_{n-j} + h\beta f(t_{n+1}, y_{n+1})
$$
**Spatial Discretization**
**Finite Volume Method**
Ensures mass conservation:
$$
\int_V \frac{\partial n}{\partial t} dV + \oint_S \mathbf{\Gamma} \cdot d\mathbf{S} = \int_V S \, dV
$$
**Mesh Requirements**
- Sheath resolution: $\Delta x < \lambda_D$
- RF skin depth: $\Delta x < \delta$
- Adaptive mesh refinement (AMR) common
**EM-Plasma Coupling**
**Iterative scheme:**
1. Solve Maxwell's equations for $\mathbf{E}$, $\mathbf{B}$
2. Update plasma transport (density, temperature)
3. Recalculate $\sigma$, $\varepsilon_{\text{plasma}}$
4. Repeat until convergence
**Advanced Topics**
**Atomic Layer Etching (ALE)**
Self-limiting reactions for atomic precision:
$$
\text{EPC} = \Theta \cdot d_{\text{ML}}
$$
**Where:**
- EPC — Etch per cycle
- $\Theta$ — Modified layer coverage fraction
- $d_{\text{ML}}$ — Monolayer thickness
**ALE Cycle**
1. **Modification step:** Reactive gas creates modified surface layer
$$\frac{d\Theta}{dt} = k_{\text{mod}}(1-\Theta)P_{\text{gas}}$$
2. **Removal step:** Ion bombardment removes modified layer only
$$\text{ER} = Y_{\text{mod}}\Gamma_i\Theta$$
**Pulsed Plasma Dynamics**
Time-modulated RF introduces:
- **Active glow:** Plasma on, high ion/radical generation
- **Afterglow:** Plasma off, selective chemistry
**Ion Energy Modulation**
By pulsing bias:
$$
\langle E_i \rangle = \frac{1}{T}\left[\int_0^{t_{\text{on}}} E_{\text{high}}dt + \int_{t_{\text{on}}}^{T} E_{\text{low}}dt\right]
$$
**High-Aspect-Ratio Etching (HAR)**
For AR > 50 (memory, 3D NAND):
**Challenges:**
- Ion angular broadening → bowing
- Neutral depletion at bottom
- Feature charging → twisting
- Mask erosion → tapering
**Ion Angular Distribution Broadening:**
$$
\sigma_{\text{effective}} = \sqrt{\sigma_{\text{sheath}}^2 + \sigma_{\text{scattering}}^2}
$$
**Neutral Flux at Bottom:**
$$
\Gamma_{\text{bottom}} \approx \Gamma_{\text{top}} \cdot K(\text{AR})
$$
**Machine Learning Integration**
**Applications:**
- Surrogate models for fast prediction
- Process optimization (Bayesian)
- Virtual metrology
- Anomaly detection
**Physics-Informed Neural Networks (PINNs):**
$$
\mathcal{L} = \mathcal{L}_{\text{data}} + \lambda \mathcal{L}_{\text{physics}}
$$
Where $\mathcal{L}_{\text{physics}}$ enforces governing equations.
**Validation and Experimental Techniques**
**Plasma Diagnostics**
| Technique | Measurement | Typical Values |
|-----------|-------------|----------------|
| **Langmuir probe** | $n_e$, $T_e$, EEDF | $10^{9}–10^{12}$ cm⁻³, 1–5 eV |
| **OES** | Relative species densities | Qualitative/semi-quantitative |
| **APMS** | Ion mass, energy | 1–500 amu, 0–500 eV |
| **LIF** | Absolute radical density | $10^{11}–10^{14}$ cm⁻³ |
| **Microwave interferometry** | $n_e$ (line-averaged) | $10^{10}–10^{12}$ cm⁻³ |
**Etch Characterization**
- **Profilometry:** Etch depth, uniformity
- **SEM/TEM:** Feature profiles, sidewall angle
- **XPS:** Surface composition
- **Ellipsometry:** Film thickness, optical properties
**Model Validation Workflow**
1. **Plasma validation:** Match $n_e$, $T_e$, species densities
2. **Flux validation:** Compare ion/neutral fluxes to wafer
3. **Etch rate validation:** Blanket wafer etch rates
4. **Profile validation:** Patterned feature cross-sections
**Key Dimensionless Numbers Summary**
| Number | Definition | Physical Meaning |
|--------|------------|------------------|
| **Knudsen** | $\text{Kn} = \lambda/L$ | Continuum vs. kinetic |
| **Damköhler** | $\text{Da} = \tau_{\text{transport}}/\tau_{\text{reaction}}$ | Transport vs. reaction limited |
| **Sticking coefficient** | $\gamma = \text{reactions}/\text{collisions}$ | Surface reactivity |
| **Aspect ratio** | $\text{AR} = \text{depth}/\text{width}$ | Feature geometry |
| **Debye number** | $N_D = n\lambda_D^3$ | Plasma ideality |
**Physical Constants**
| Constant | Symbol | Value |
|----------|--------|-------|
| Elementary charge | $e$ | $1.602 \times 10^{-19}$ C |
| Electron mass | $m_e$ | $9.109 \times 10^{-31}$ kg |
| Proton mass | $m_p$ | $1.673 \times 10^{-27}$ kg |
| Boltzmann constant | $k_B$ | $1.381 \times 10^{-23}$ J/K |
| Vacuum permittivity | $\varepsilon_0$ | $8.854 \times 10^{-12}$ F/m |
| Vacuum permeability | $\mu_0$ | $4\pi \times 10^{-7}$ H/m |