Plasma Etch Modeling

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

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

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

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

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:

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

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

Dimensionless Numbers Summary

Key Physical Constants