Aluminum Metal Etch Mathematical Modeling

1. Overview

1.1 Why Aluminum Etch Modeling is Complex

Aluminum etching (typically using $\text{Cl}_2/\text{BCl}_3$ plasmas) involves multiple coupled physical and chemical phenomena:

• Plasma generation and transport → determines species fluxes to wafer
• Ion-surface interactions → physical and chemical mechanisms
• Surface reactions → Langmuir-Hinshelwood kinetics
• Feature-scale evolution → profile development inside trenches/vias
• Redeposition and passivation → sidewall chemistry

1.2 Fundamental Reaction

The basic aluminum chlorination reaction:

$$ \text{Al} + 3\text{Cl} \rightarrow \text{AlCl}_3 \uparrow $$

Complications requiring sophisticated modeling:

• Breaking through native $\text{Al}_2\text{O}_3$ layer (15-30 Å)
• Maintaining profile anisotropy
• Controlling selectivity to mask and underlayers
• Managing Cu residues in Al-Cu alloys

2. Kinetic and Chemical Rate Modeling

2.1 General Etch Rate Formulation

A comprehensive etch rate model combines three primary mechanisms:

$$ ER = \underbrace{k_{th} \cdot \Gamma_{Cl} \cdot f(\theta)}_{\text{thermal chemical}} + \underbrace{Y_s \cdot \Gamma_{ion} \cdot \sqrt{E_{ion}}}_{\text{physical sputtering}} + \underbrace{\beta \cdot \Gamma_{ion}^a \cdot \Gamma_{Cl}^b \cdot E_{ion}^c}_{\text{ion-enhanced (synergistic)}} $$

Parameter Definitions:

2.2 Surface Coverage Dynamics

The reactive site balance follows Langmuir-Hinshelwood kinetics:

$$ \frac{d\theta}{dt} = k_{ads} \cdot \Gamma_{Cl} \cdot (1-\theta) - k_{des} \cdot \theta \cdot \exp\left(-\frac{E_d}{k_B T}\right) - Y_{react}(\theta, E_{ion}) \cdot \Gamma_{ion} \cdot \theta $$

Term-by-term breakdown:

• Term 1: $k_{ads} \cdot \Gamma_{Cl} \cdot (1-\theta)$ — Adsorption rate (proportional to empty sites)
• Term 2: $k_{des} \cdot \theta \cdot \exp(-E_d/k_B T)$ — Thermal desorption (Arrhenius)
• Term 3: $Y_{react} \cdot \Gamma_{ion} \cdot \theta$ — Ion-induced reaction/removal

Steady-State Solution ($d\theta/dt = 0$):

$$ \theta_{ss} = \frac{k_{ads} \cdot \Gamma_{Cl}}{k_{ads} \cdot \Gamma_{Cl} + k_{des} \cdot e^{-E_d/k_B T} + Y_{react} \cdot \Gamma_{ion}} $$

2.3 Temperature Dependence

All rate constants follow Arrhenius behavior:

$$ k_i(T) = A_i \cdot \exp\left(-\frac{E_{a,i}}{k_B T}\right) $$

Typical activation energies for aluminum etching:

• Ion-enhanced reactions: $E_a \approx 0.1 - 0.3 \text{ eV}$
• Purely thermal processes: $E_a \approx 0.5 - 1.0 \text{ eV}$
• Chlorine desorption: $E_d \approx 0.3 - 0.5 \text{ eV}$

2.4 Complete Etch Rate Expression

Combining all terms with explicit dependencies:

$$ ER(T, \Gamma_{ion}, \Gamma_{Cl}, E_{ion}) = A_1 e^{-E_1/k_B T} \Gamma_{Cl} \theta + Y_0 \Gamma_{ion} \sqrt{E_{ion}} + A_2 e^{-E_2/k_B T} \Gamma_{ion}^{0.5} \Gamma_{Cl}^{0.5} E_{ion}^{0.5} $$

3. Ion-Surface Interaction Physics

3.1 Ion Energy Distribution Function (IEDF)

For RF-biased electrodes, the IEDF is approximately bimodal:

$$ f(E) \propto \frac{1}{\sqrt{|E - E_{dc}|}} \quad \text{for } E_{dc} - E_{rf} < E < E_{dc} + E_{rf} $$

Key parameters:

• $E_{dc} = e \cdot V_{dc}$ — DC self-bias energy
• $E_{rf} = e \cdot V_{rf}$ — RF amplitude energy
• Peak separation: $\Delta E = 2 E_{rf}$

Collisional effects:

In collisional sheaths, charge-exchange collisions broaden the distribution:

$$ f(E) \propto \exp\left(-\frac{E}{\bar{E}}\right) \cdot \left[1 + \text{erf}\left(\frac{E - E_{dc}}{\sigma_E}\right)\right] $$

3.2 Ion Angular Distribution Function (IADF)

The angular spread is approximately Gaussian:

$$ f(\theta) = \frac{1}{\sqrt{2\pi}\sigma_\theta} \exp\left(-\frac{\theta^2}{2\sigma_\theta^2}\right) $$

Angular spread calculation:

$$ \sigma_\theta \approx \sqrt{\frac{k_B T_i}{e V_{sheath}}} \approx \arctan\left(\sqrt{\frac{T_i}{V_{sheath}}}\right) $$

Typical values:

• Ion temperature: $T_i \approx 0.05 - 0.5 \text{ eV}$
• Sheath voltage: $V_{sheath} \approx 50 - 500 \text{ V}$
• Angular spread: $\sigma_\theta \approx 2° - 5°$

3.3 Physical Sputtering Yield

Yamamura Formula (Angular Dependence)

$$ Y(\theta) = Y(0°) \cdot \cos^{-f}(\theta) \cdot \exp\left[b\left(1 - \frac{1}{\cos\theta}\right)\right] $$

Parameters for aluminum:

• $f \approx 1.5 - 2.0$
• $b \approx 0.1 - 0.3$ (depends on ion/target mass ratio)
• Maximum yield typically at $\theta \approx 60° - 70°$

Sigmund Theory (Energy Dependence)

$$ Y(E) = \frac{0.042 \cdot Q \cdot \alpha(M_2/M_1) \cdot S_n(E)}{U_s} $$

Where:

• $S_n(E)$ = nuclear stopping power (Thomas-Fermi)
• $U_s = 3.4 \text{ eV}$ (surface binding energy for Al)
• $Q$ = dimensionless factor ($\approx 1$ for metals)
• $\alpha$ = mass-dependent parameter
• $M_1, M_2$ = projectile and target masses

Nuclear Stopping Power

$$ S_n(\epsilon) = \frac{0.5 \ln(1 + 1.2288\epsilon)}{\epsilon + 0.1728\sqrt{\epsilon} + 0.008\epsilon^{0.1504}} $$

With reduced energy:

$$ \epsilon = \frac{M_2 E}{(M_1 + M_2) Z_1 Z_2 e^2} \cdot \frac{a_{TF}}{1} $$

3.4 Ion-Enhanced Etching Yield

The total etch yield combines mechanisms:

$$ Y_{total} = Y_{physical} + Y_{chemical} + Y_{synergistic} $$

Synergistic enhancement factor:

$$ \eta = \frac{Y_{total}}{Y_{physical} + Y_{chemical}} > 1 $$

For Al/Cl₂ systems, $\eta$ can exceed 10 under optimal conditions.

4. Plasma Modeling (Reactor Scale)

4.1 Species Continuity Equations

For each species $i$ (electrons, ions, neutrals):

$$ \frac{\partial n_i}{\partial t} + abla \cdot \vec{\Gamma}_i = S_i - L_i $$

Flux expressions:

• Drift-diffusion: $\vec{\Gamma}_i = -D_i

abla n_i + \mu_i n_i \vec{E}$

• Full momentum: $\vec{\Gamma}_i = n_i \vec{v}_i$ with momentum equation

Source/sink terms:

$$ S_i = \sum_j k_{ij} n_j n_e \quad \text{(ionization, dissociation)} $$

$$ L_i = \sum_j k_{ij}^{loss} n_i n_j \quad \text{(recombination, attachment)} $$

4.2 Electron Energy Balance

$$ \frac{\partial}{\partial t}\left(\frac{3}{2} n_e k_B T_e\right) + abla \cdot \vec{Q}_e = P_{abs} - P_{loss} $$

Heat flux:

$$ \vec{Q}_e = \frac{5}{2} k_B T_e \vec{\Gamma}_e - \kappa_e abla T_e $$

Power absorption (ICP):

$$ P_{abs} = \frac{1}{2} \text{Re}(\sigma_p) |E|^2 $$

Collisional losses:

$$ P_{loss} = \sum_j n_e n_j k_j \varepsilon_j $$

Where $\varepsilon_j$ is the energy loss per collision event $j$.

4.3 Plasma Conductivity

$$ \sigma_p = \frac{n_e e^2}{m_e( u_m + i\omega)} $$

Skin depth:

$$ \delta = \sqrt{\frac{2}{\omega \mu_0 \text{Re}(\sigma_p)}} $$

4.4 Electromagnetic Field Equations

Maxwell's equations (frequency domain):

$$

abla \times \vec{E} = -i\omega \vec{B} $$

$$

abla \times \vec{B} = \mu_0 \sigma_p \vec{E} + i\omega \mu_0 \epsilon_0 \vec{E} $$

Wave equation:

$$

abla^2 \vec{E} + \left(\frac{\omega^2}{c^2} - i\omega\mu_0\sigma_p\right)\vec{E} = 0 $$

4.5 Sheath Physics

Child-Langmuir Law (Collisionless Sheath)

$$ J_{ion} = \frac{4\epsilon_0}{9}\sqrt{\frac{2e}{M}} \cdot \frac{V_s^{3/2}}{s^2} $$

Where:

• $J_{ion}$ = ion current density
• $V_s$ = sheath voltage
• $s$ = sheath thickness
• $M$ = ion mass

Bohm Criterion

Ions must enter sheath with velocity:

$$ v_{Bohm} = \sqrt{\frac{k_B T_e}{M}} $$

Ion flux at sheath edge:

$$ \Gamma_{ion} = n_s \cdot v_{Bohm} = 0.61 \cdot n_0 \sqrt{\frac{k_B T_e}{M}} $$

Sheath Thickness

$$ s \approx \lambda_D \cdot \left(\frac{2 e V_s}{k_B T_e}\right)^{3/4} $$

Debye length:

$$ \lambda_D = \sqrt{\frac{\epsilon_0 k_B T_e}{n_e e^2}} $$

5. Feature-Scale Profile Evolution

5.1 Level Set Method

The surface is represented implicitly by $\phi(\vec{r}, t) = 0$:

$$ \frac{\partial \phi}{\partial t} + V_n | abla \phi| = 0 $$

Normal velocity calculation:

$$ V_n(\vec{r}) = \int_0^{E_{max}} \int_0^{\theta_{max}} Y(E, \theta_{local}) \cdot f_{IEDF}(E) \cdot f_{IADF}(\theta) \cdot \Gamma_{ion}(\vec{r}) \, dE \, d\theta $$

Plus contributions from:

• Neutral chemical etching
• Redeposition
• Surface diffusion

5.2 Hamilton-Jacobi Formulation

$$ \frac{\partial \phi}{\partial t} + H( abla \phi, \vec{r}, t) = 0 $$

Hamiltonian for etch:

$$ H = V_n \sqrt{\phi_x^2 + \phi_y^2 + \phi_z^2} $$

With $V_n$ dependent on:

• Local surface normal: $\hat{n} = -

abla\phi / | abla\phi|$

• Local fluxes: $\Gamma(\vec{r})$
• Local angles: $\theta = \arccos(\hat{n} \cdot \hat{z})$

5.3 Visibility and View Factors

Direct Flux

The flux reaching a point inside a feature depends on solid angle visibility:

$$ \Gamma_{direct}(\vec{r}) = \int_{\Omega_{visible}} \Gamma_0 \cdot \cos\theta \cdot \frac{d\Omega}{\pi} $$

Reflected/Reemitted Flux

For neutrals with sticking coefficient $s$:

$$ \Gamma_{total}(\vec{r}) = \Gamma_{direct}(\vec{r}) + (1-s) \cdot \Gamma_{reflected}(\vec{r}) $$

This leads to coupled integral equations:

$$ \Gamma(\vec{r}) = \Gamma_{plasma}(\vec{r}) + (1-s) \int_{S'} K(\vec{r}, \vec{r'}) \Gamma(\vec{r'}) dS' $$

Kernel function:

$$ K(\vec{r}, \vec{r'}) = \frac{\cos\theta \cos\theta'}{\pi |\vec{r} - \vec{r'}|^2} \cdot V(\vec{r}, \vec{r'}) $$

Where $V(\vec{r}, \vec{r'})$ is the visibility function (1 if visible, 0 otherwise).

5.4 Aspect Ratio Dependent Etching (ARDE)

Empirical model:

$$ \frac{ER(AR)}{ER_0} = \frac{1}{1 + (AR/AR_c)^n} $$

Where:

• $AR = \text{depth}/\text{width}$ (aspect ratio)
• $AR_c$ = critical aspect ratio (process-dependent)
• $n \approx 1 - 2$

Knudsen transport model:

$$ \Gamma_{neutral}(z) = \Gamma_0 \cdot \frac{W}{W + \alpha \cdot z} $$

Where:

• $z$ = feature depth
• $W$ = feature width
• $\alpha$ = Clausing factor (depends on geometry and sticking)

Clausing factor for cylinder:

$$ \alpha = \frac{8}{3} \cdot \frac{1 - s}{s} $$

6. Aluminum-Specific Phenomena

6.1 Native Oxide Breakthrough

$\text{Al}_2\text{O}_3$ (15-30 Å native oxide) requires physical sputtering:

$$ ER_{oxide} \approx Y_{\text{BCl}_3^+}(E) \cdot \Gamma_{ion} $$

Why BCl₃ is critical:

1. Heavy $\text{BCl}_3^+$ ions provide efficient momentum transfer 2. BCl₃ scavenges oxygen chemically:

$$ 2\text{BCl}_3 + \text{Al}_2\text{O}_3 \rightarrow 2\text{AlCl}_3 \uparrow + \text{B}_2\text{O}_3 $$

Breakthrough time:

$$ t_{breakthrough} = \frac{d_{oxide}}{ER_{oxide}} = \frac{d_{oxide}}{Y_{BCl_3^+} \cdot \Gamma_{ion}} $$

6.2 Sidewall Passivation Dynamics

Anisotropic profiles require passivation of sidewalls:

$$ \frac{d\tau_{pass}}{dt} = R_{dep}(\Gamma_{redeposition}, s_{stick}) - R_{removal}(\Gamma_{ion}, \theta_{sidewall}) $$

Deposition sources:

• $\text{AlCl}_x$ redeposition from etch products
• Photoresist erosion products (C, H, O, N)
• Intentional additives: $\text{N}_2 \rightarrow \text{AlN}$ formation

Why sidewalls are protected:

At grazing incidence ($\theta \approx 85° - 90°$):

• Ion flux geometric factor: $\Gamma_{sidewall} = \Gamma_0 \cdot \cos(90° - \alpha) \approx \Gamma_0 \cdot \sin\alpha$
• For $\alpha = 5°$: $\Gamma_{sidewall} \approx 0.09 \cdot \Gamma_0$

• Sputtering yield at grazing incidence approaches zero
• Net passivation accumulates → blocks lateral etching

6.3 Notching and Charging Effects

At dielectric interfaces, differential charging causes ion deflection:

Surface charge evolution:

$$ \frac{d\sigma}{dt} = J_{ion} - J_{electron} $$

Where:

• $\sigma$ = surface charge density (C/cm²)
• $J_{ion}$ = ion current (always positive)
• $J_{electron}$ = electron current (depends on local potential)

Local electric field:

$$ \vec{E}_{charging} = - abla V_{charging} $$

Laplace equation in feature:

$$

abla^2 V = -\frac{\rho}{\epsilon_0} \quad \text{(with } \rho = 0 \text{ in vacuum)} $$

Modified ion trajectory:

$$ m \frac{d^2\vec{r}}{dt^2} = e\left(\vec{E}_{sheath} + \vec{E}_{charging}\right) $$

Result: Ions deflect toward charged surfaces → notching at feature bottom.

Mitigation strategies:

• Pulsed plasmas (allow electron neutralization)
• Low-frequency bias (time for charge equilibration)
• Conductive underlayers

6.4 Copper Residue Formation (Al-Cu Alloys)

Al-Cu alloys (0.5-4% Cu) leave Cu residues because Cu chlorides are less volatile:

Volatility comparison:

Residue accumulation rate:

$$ \frac{d[\text{Cu}]_{surface}}{dt} = x_{Cu} \cdot ER_{Al} - ER_{Cu} $$

Where:

• $x_{Cu}$ = Cu atomic fraction in alloy
• At low temperature: $ER_{Cu} \ll x_{Cu} \cdot ER_{Al}$

Solutions:

• Elevated substrate temperature ($>$150°C)
• Increased BCl₃ fraction
• Post-etch treatments

7. Numerical Methods

7.1 Level Set Discretization

Upwind Finite Differences

Using Hamilton-Jacobi ENO (Essentially Non-Oscillatory) schemes:

$$ \phi_i^{n+1} = \phi_i^n - \Delta t \cdot H(\phi_x^-, \phi_x^+, \phi_y^-, \phi_y^+) $$

One-sided derivatives:

$$ \phi_x^- = \frac{\phi_i - \phi_{i-1}}{\Delta x}, \quad \phi_x^+ = \frac{\phi_{i+1} - \phi_i}{\Delta x} $$

Godunov flux for $H = V_n | abla\phi|$:

$$ H^{Godunov} = \begin{cases} V_n \sqrt{\max(\phi_x^{-,+},0)^2 + \max(\phi_y^{-,+},0)^2} & \text{if } V_n > 0 \\ V_n \sqrt{\max(\phi_x^{+,-},0)^2 + \max(\phi_y^{+,-},0)^2} & \text{if } V_n < 0 \end{cases} $$

Reinitialization

Maintain $| abla\phi| = 1$ using:

$$ \frac{\partial \phi}{\partial \tau} = \text{sign}(\phi_0)(1 - | abla\phi|) $$

Iterate in pseudo-time $\tau$ until convergence.

7.2 Monte Carlo Feature-Scale Simulation

Algorithm:

1. INITIALIZE surface mesh 2. FOR each time step: a. FOR i = 1 to N_particles:

• Sample particle from IEDF, IADF
• Launch from plasma boundary
• TRACE trajectory until surface hit
• APPLY reaction probability:
• Etch (remove cell) with probability P_etch
• Reflect with probability P_reflect
• Deposit with probability P_deposit

b. UPDATE surface mesh c. CHECK for convergence 3. OUTPUT final profile

Variance reduction techniques:

• Importance sampling: Weight particles toward features of interest
• Particle splitting: Increase statistics in critical regions
• Russian roulette: Terminate low-weight particles probabilistically

7.3 Coupled Multi-Scale Modeling

Coupling strategy:

$$ \text{Reactor} \xrightarrow{\Gamma_i, f(E), f(\theta)} \text{Feature} \xrightarrow{ER(\vec{r})} \text{Reactor} $$

7.4 Plasma Solver Discretization

Finite element for Poisson's equation:

$$

abla \cdot (\epsilon abla V) = -\rho $$

Weak form:

$$ \int_\Omega \epsilon abla V \cdot abla w \, d\Omega = \int_\Omega \rho \, w \, d\Omega $$

Finite volume for transport:

$$ \frac{d(n_i V_j)}{dt} = -\sum_{faces} \Gamma_i \cdot \hat{n} \cdot A + S_i V_j $$

8. Process Window and Optimization

8.1 Response Surface Modeling

Quadratic response surface:

$$ ER = \beta_0 + \sum_{i=1}^{k} \beta_i x_i + \sum_{i=1}^{k} \beta_{ii} x_i^2 + \sum_{i

Key process variables ($x_i$):

• Pressure (mTorr)
• RF source power (W)
• RF bias power (W)
• Cl₂ flow (sccm)
• BCl₃ flow (sccm)
• Temperature (°C)

Matrix formulation:

$$ \vec{y} = X\vec{\beta} + \vec{\epsilon} $$

Least squares solution:

$$ \hat{\vec{\beta}} = (X^T X)^{-1} X^T \vec{y} $$

8.2 Multi-Objective Optimization

Desirability function approach:

$$ D = \left(\prod_{i=1}^{n} d_i^{w_i}\right)^{1/\sum w_i} $$

Individual desirabilities:

$$ d_i = \begin{cases} 0 & \text{if } y_i < L_i \\ \left(\frac{y_i - L_i}{T_i - L_i}\right)^s & \text{if } L_i \leq y_i \leq T_i \\ 1 & \text{if } y_i > T_i \end{cases} $$

Optimization problem:

$$ \max_{\vec{x}} D(\vec{x}) $$

Subject to:

• $85° < \text{sidewall angle} < 90°$
• $\text{Selectivity}_{Al:resist} > 3:1$
• $\text{Selectivity}_{Al:TiN} > 10:1$
• $\text{Uniformity} < 3\%$ (1σ)

8.3 Virtual Metrology

Prediction model:

$$ \vec{y}_{etch} = f_{ML}\left(\vec{x}_{recipe}, \vec{x}_{OES}, \vec{x}_{chamber}\right) $$

Input features:

• Recipe: Power, pressure, flows, time
• OES: Emission line intensities (e.g., Al 396nm, Cl 837nm)
• Chamber: Impedance, temperature, previous wafer history

Machine learning approaches:

• Neural networks (for complex nonlinear relationships)
• Gaussian processes (with uncertainty quantification)
• Partial least squares (for high-dimensional, correlated inputs)

8.4 Run-to-Run Control

EWMA (Exponentially Weighted Moving Average) controller:

$$ \vec{x}_{k+1} = \vec{x}_k + \Lambda G^{-1}(\vec{y}_{target} - \vec{y}_k) $$

Where:

• $\Lambda$ = diagonal weighting matrix (0 < λ < 1)
• $G$ = process gain matrix ($\partial y / \partial x$)

Drift compensation:

$$ \vec{x}_{k+1} = \vec{x}_k + \Lambda_1 G^{-1}(\vec{y}_{target} - \vec{y}_k) + \Lambda_2 (\vec{x}_{k} - \vec{x}_{k-1}) $$

9. Equations:

abla\phi| = 0$ |

abla \cdot \vec{\Gamma}_i = S_i - L_i$ |

10. Modern Developments

10.1 Machine Learning Integration

Applications:

• Yield prediction: Neural networks trained on MD simulation data
• Surrogate models: Replace expensive PDE solvers for real-time optimization
• Process control: Reinforcement learning for adaptive recipes

Example: Gaussian Process for Etch Rate:

$$ ER(\vec{x}) \sim \mathcal{GP}\left(m(\vec{x}), k(\vec{x}, \vec{x}')\right) $$

With squared exponential kernel:

$$ k(\vec{x}, \vec{x}') = \sigma_f^2 \exp\left(-\frac{|\vec{x} - \vec{x}'|^2}{2\ell^2}\right) $$

10.2 Atomistic-Continuum Bridging

ReaxFF molecular dynamics:

• Reactive force fields for Al-Cl-O systems
• Calculate fundamental yields and sticking coefficients
• Feed into continuum models

DFT calculations:

• Adsorption energies: $E_{ads} = E_{surface+adsorbate} - E_{surface} - E_{adsorbate}$
• Activation barriers via NEB (Nudged Elastic Band)
• Electronic structure effects on reactivity

10.3 Digital Twins

Components:

• Real-time sensor data ingestion
• Physics-based + ML hybrid models
• Predictive maintenance algorithms
• Virtual process development

Update equation:

$$ \vec{\theta}_{model}^{(k+1)} = \vec{\theta}_{model}^{(k)} + K_k \left(\vec{y}_{measured} - \vec{y}_{predicted}\right) $$

10.4 Uncertainty Quantification

Bayesian calibration:

$$ p(\vec{\theta}|\vec{y}) \propto p(\vec{y}|\vec{\theta}) \cdot p(\vec{\theta}) $$

Propagation through models:

$$ \text{Var}(y) \approx \sum_i \left(\frac{\partial y}{\partial \theta_i}\right)^2 \text{Var}(\theta_i) $$

Monte Carlo uncertainty:

$$ \bar{y} \pm t_{\alpha/2} \cdot \frac{s}{\sqrt{N}} $$

Physical Constants

Process Conditions