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

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

Key Dimensionless Numbers Summary

Physical Constants