Mathematical Modeling of Metal Deposition in Semiconductor Manufacturing

1. Overview: Metal Deposition Processes

Metal deposition is a critical step in semiconductor fabrication, creating interconnects, contacts, barrier layers, and various metallic structures. The primary deposition methods require distinct mathematical treatments:

2. Transport Phenomena Models

2.1 Gas-Phase Transport (CVD/PECVD)

The precursor concentration field follows the convection-diffusion-reaction equation:

$$ \frac{\partial C}{\partial t} + \mathbf{v} \cdot abla C = D abla^2 C + R_{gas} $$

Where:

• $C$ — precursor concentration (mol/m³)
• $\mathbf{v}$ — velocity field vector (m/s)
• $D$ — diffusion coefficient (m²/s)
• $R_{gas}$ — gas-phase reaction source term (mol/m³$\cdot$s)

2.2 Flow Field Equations

The incompressible Navier-Stokes equations govern the velocity field:

$$ \rho \left( \frac{\partial \mathbf{v}}{\partial t} + \mathbf{v} \cdot abla \mathbf{v} \right) = - abla p + \mu abla^2 \mathbf{v} $$

With continuity equation:

$$

abla \cdot \mathbf{v} = 0 $$

Where:

• $\rho$ — gas density (kg/m³)
• $p$ — pressure (Pa)
• $\mu$ — dynamic viscosity (Pa$\cdot$s)

2.3 Knudsen Number and Transport Regimes

At low pressures, the Knudsen number determines the transport regime:

$$ Kn = \frac{\lambda}{L} = \frac{k_B T}{\sqrt{2} \pi d^2 p L} $$

Where:

• $\lambda$ — mean free path (m)
• $L$ — characteristic length (m)
• $k_B$ — Boltzmann constant ($1.38 \times 10^{-23}$ J/K)
• $T$ — temperature (K)
• $d$ — molecular diameter (m)
• $p$ — pressure (Pa)

Transport regime classification:

• $Kn < 0.01$ — Continuum regime → Navier-Stokes CFD
• $0.01 < Kn < 0.1$ — Slip flow regime → Modified NS with slip boundary conditions
• $0.1 < Kn < 10$ — Transitional regime → DSMC, Boltzmann equation
• $Kn > 10$ — Free molecular regime → Ballistic/Monte Carlo methods

3. Surface Reaction Kinetics

3.1 Langmuir-Hinshelwood Mechanism

For bimolecular surface reactions (common in CVD):

$$ r = \frac{k \cdot K_A K_B \cdot p_A p_B}{(1 + K_A p_A + K_B p_B)^2} $$

Where:

• $r$ — reaction rate (mol/m²$\cdot$s)
• $k$ — surface reaction rate constant (mol/m²$\cdot$s)
• $K_A, K_B$ — adsorption equilibrium constants (Pa⁻¹)
• $p_A, p_B$ — partial pressures of reactants A and B (Pa)

3.2 Sticking Coefficient Model

The probability that an impinging molecule adsorbs on the surface:

$$ S = S_0 \exp\left( -\frac{E_a}{k_B T} \right) \cdot f(\theta) $$

Where:

• $S$ — sticking coefficient (dimensionless)
• $S_0$ — pre-exponential sticking factor
• $E_a$ — activation energy (J)
• $f(\theta) = (1 - \theta)^n$ — site blocking function
• $\theta$ — surface coverage (dimensionless, 0 to 1)
• $n$ — order of site blocking

3.3 Arrhenius Temperature Dependence

$$ k(T) = A \exp\left( -\frac{E_a}{RT} \right) $$

Where:

• $A$ — pre-exponential factor (frequency factor)
• $E_a$ — activation energy (J/mol)
• $R$ — universal gas constant (8.314 J/mol$\cdot$K)
• $T$ — absolute temperature (K)

4. Film Growth Models

4.1 Continuum Surface Evolution

Edwards-Wilkinson Equation (Linear Growth)

$$ \frac{\partial h}{\partial t} = u abla^2 h + F + \eta(\mathbf{x}, t) $$

Kardar-Parisi-Zhang (KPZ) Equation (Nonlinear Growth)

$$ \frac{\partial h}{\partial t} = u abla^2 h + \frac{\lambda}{2} | abla h|^2 + F + \eta $$

Where:

• $h(\mathbf{x}, t)$ — surface height at position $\mathbf{x}$ and time $t$
• $

u$ — surface diffusion coefficient (m²/s)

• $\lambda$ — nonlinear growth parameter
• $F$ — mean deposition flux (m/s)
• $\eta$ — stochastic noise term (Gaussian white noise)

4.2 Scaling Relations

Surface roughness evolves according to:

$$ W(L, t) = L^\alpha f\left( \frac{t}{L^z} \right) $$

Where:

• $W$ — interface width (roughness)
• $L$ — system size
• $\alpha$ — roughness exponent
• $z$ — dynamic exponent
• $f$ — scaling function

5. Step Coverage and Conformality

5.1 Thiele Modulus

For high-aspect-ratio features, the Thiele modulus determines conformality:

$$ \phi = L \sqrt{\frac{k_s}{D_{eff}}} $$

Where:

• $\phi$ — Thiele modulus (dimensionless)
• $L$ — feature depth (m)
• $k_s$ — surface reaction rate constant (m/s)
• $D_{eff}$ — effective diffusivity (m²/s)

Step coverage regimes:

• $\phi \ll 1$ — Reaction-limited → Excellent conformality
• $\phi \gg 1$ — Transport-limited → Poor step coverage (bread-loafing)

5.2 Knudsen Diffusion in Trenches

$$ D_K = \frac{w}{3} \sqrt{\frac{8 R T}{\pi M}} $$

Where:

• $D_K$ — Knudsen diffusion coefficient (m²/s)
• $w$ — trench width (m)
• $R$ — universal gas constant (J/mol$\cdot$K)
• $T$ — temperature (K)
• $M$ — molecular weight (kg/mol)

5.3 Feature-Scale Concentration Profile

Solving for concentration in a trench with reactive walls:

$$ D_{eff} \frac{d^2 C}{dy^2} = \frac{2 k_s C}{w} $$

General solution:

$$ C(y) = C_0 \frac{\cosh\left( \phi \frac{L - y}{L} \right)}{\cosh(\phi)} $$

6. Atomic Layer Deposition (ALD) Models

6.1 Self-Limiting Surface Kinetics

Surface site balance equation:

$$ \frac{d\theta}{dt} = k_a C (1 - \theta) - k_d \theta $$

Where:

• $\theta$ — fractional surface coverage
• $k_a$ — adsorption rate constant (m³/mol$\cdot$s)
• $k_d$ — desorption rate constant (s⁻¹)
• $C$ — gas-phase precursor concentration (mol/m³)

At equilibrium saturation:

$$ \theta_{eq} = \frac{k_a C}{k_a C + k_d} \approx 1 \quad \text{(for strong chemisorption)} $$

6.2 Growth Per Cycle (GPC)

$$ \text{GPC} = \Gamma_0 \cdot \Omega \cdot \eta $$

Where:

• $\Gamma_0$ — surface site density (sites/m²)
• $\Omega$ — volume per deposited atom (m³)
• $\eta$ — reaction efficiency (dimensionless)

6.3 Saturation Dose-Time Relationship

$$ \theta(t) = 1 - \exp\left( -\frac{S \cdot \Phi \cdot t}{\Gamma_0} \right) $$

Impingement flux from kinetic theory:

$$ \Phi = \frac{p}{\sqrt{2 \pi m k_B T}} $$

Where:

• $\Phi$ — molecular impingement flux (molecules/m²$\cdot$s)
• $p$ — precursor partial pressure (Pa)
• $m$ — molecular mass (kg)

7. Plasma Modeling (PVD/PECVD)

7.1 Plasma Sheath Physics

Child-Langmuir law for ion current density:

$$ J_{ion} = \frac{4 \varepsilon_0}{9} \sqrt{\frac{2e}{M_i}} \frac{V_s^{3/2}}{d_s^2} $$

Where:

• $J_{ion}$ — ion current density (A/m²)
• $\varepsilon_0$ — vacuum permittivity ($8.85 \times 10^{-12}$ F/m)
• $e$ — elementary charge ($1.6 \times 10^{-19}$ C)
• $M_i$ — ion mass (kg)
• $V_s$ — sheath voltage (V)
• $d_s$ — sheath thickness (m)

7.2 Ion Energy at Substrate

$$ \varepsilon_{ion} \approx e V_s + \frac{1}{2} M_i v_{Bohm}^2 $$

Bohm velocity:

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

Where:

• $T_e$ — electron temperature (K or eV)

7.3 Sputtering Yield (Sigmund Formula)

$$ Y(E) = \frac{3 \alpha}{4 \pi^2} \cdot \frac{4 M_1 M_2}{(M_1 + M_2)^2} \cdot \frac{E}{U_0} $$

Where:

• $Y$ — sputtering yield (atoms/ion)
• $\alpha$ — dimensionless factor (~0.2–0.4)
• $M_1$ — incident ion mass
• $M_2$ — target atom mass
• $E$ — incident ion energy (eV)
• $U_0$ — surface binding energy (eV)

7.4 Electron Energy Distribution Function (EEDF)

The Boltzmann equation in energy space:

$$ \frac{\partial f}{\partial t} + \mathbf{v} \cdot abla f + \frac{e \mathbf{E}}{m_e} \cdot abla_v f = C[f] $$

Where:

• $f$ — electron energy distribution function
• $\mathbf{E}$ — electric field
• $m_e$ — electron mass
• $C[f]$ — collision integral

8. MDP: Markov Decision Process for Process Control

8.1 MDP Formulation

A Markov Decision Process is defined by the tuple:

$$ \mathcal{M} = (S, A, P, R, \gamma) $$

Components in semiconductor context:

• State space $S$: Film thickness, resistivity, uniformity, equipment state, wafer position
• Action space $A$: Temperature, pressure, flow rates, RF power, deposition time
• Transition probability $P(s' | s, a)$: Stochastic process model
• Reward function $R(s, a)$: Yield, uniformity, throughput, quality metrics
• Discount factor $\gamma$: Time preference (typically 0.9–0.99)

8.2 Bellman Optimality Equation

$$ V^(s) = \max_{a \in A} \left[ R(s, a) + \gamma \sum_{s'} P(s' | s, a) V^(s') \right] $$

Q-function formulation:

$$ Q^(s, a) = R(s, a) + \gamma \sum_{s'} P(s' | s, a) \max_{a'} Q^(s', a') $$

8.3 Run-to-Run (R2R) Control

Optimal recipe adjustment after each wafer:

$$ \mathbf{u}_{k+1} = \mathbf{u}_k + \mathbf{K} (\mathbf{y}_{target} - \mathbf{y}_k) $$

Where:

• $\mathbf{u}_k$ — process recipe parameters at run $k$
• $\mathbf{y}_k$ — measured output at run $k$
• $\mathbf{K}$ — controller gain matrix (from MDP policy optimization)

8.4 Reinforcement Learning Approaches

9. Electrochemical Deposition (Copper Damascene)

9.1 Butler-Volmer Equation

$$ i = i_0 \left[ \exp\left( \frac{\alpha_a F \eta}{RT} \right) - \exp\left( -\frac{\alpha_c F \eta}{RT} \right) \right] $$

Where:

• $i$ — current density (A/m²)
• $i_0$ — exchange current density (A/m²)
• $\alpha_a, \alpha_c$ — anodic and cathodic transfer coefficients
• $F$ — Faraday constant (96,485 C/mol)
• $\eta = E - E_{eq}$ — overpotential (V)
• $R$ — gas constant (J/mol$\cdot$K)
• $T$ — temperature (K)

9.2 Mass Transport Limited Current

$$ i_L = \frac{n F D C_b}{\delta} $$

Where:

• $i_L$ — limiting current density (A/m²)
• $n$ — number of electrons transferred
• $D$ — diffusion coefficient of Cu²⁺ (m²/s)
• $C_b$ — bulk concentration (mol/m³)
• $\delta$ — diffusion layer thickness (m)

9.3 Nernst-Planck Equation

$$ \mathbf{J}_i = -D_i abla C_i - \frac{z_i F D_i}{RT} C_i abla \phi + C_i \mathbf{v} $$

Where:

• $\mathbf{J}_i$ — flux of species $i$
• $z_i$ — charge number
• $\phi$ — electric potential

9.4 Superfilling (Bottom-Up Fill)

The curvature-enhanced accelerator mechanism:

$$ v_n = v_0 (1 + \kappa \cdot \Gamma_{acc}) $$

Where:

• $v_n$ — local growth velocity normal to surface
• $v_0$ — baseline growth velocity
• $\kappa$ — local surface curvature (1/m)
• $\Gamma_{acc}$ — accelerator surface concentration

10. Multiscale Modeling Framework

10.1 Hierarchical Scale Integration

-
┌──────────────────────────────────────────────────────────────┐
│                     REACTOR SCALE                            │
│            CFD: Flow, temperature, concentration             │
│                   Time: seconds | Length: cm                 │
└─────────────────────────┬────────────────────────────────────┘
                          │ Boundary fluxes
                          ▼
┌──────────────────────────────────────────────────────────────┐
│                     FEATURE SCALE                            │
│        Level-set / String method for surface evolution       │
│                   Time: seconds | Length: $\mu$m             │
└─────────────────────────┬────────────────────────────────────┘
                          │ Local rates
                          ▼
┌──────────────────────────────────────────────────────────────┐
│                    MESOSCALE (kMC)                           │
│          Kinetic Monte Carlo: nucleation, island growth      │
│                    Time: ms | Length: nm                     │
└─────────────────────────┬────────────────────────────────────┘
                          │ Rate parameters
                          ▼
┌──────────────────────────────────────────────────────────────┐
│                  ATOMISTIC (MD/DFT)                          │
│       Molecular dynamics, ab initio: binding energies,       │
│              diffusion barriers, reaction paths              │
│                    Time: ps | Length: Å                      │
└──────────────────────────────────────────────────────────────┘

10.2 Kinetic Monte Carlo (kMC)

Event rate from transition state theory:

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

Total rate and time step:

$$ k_{total} = \sum_i k_i, \quad \Delta t = -\frac{\ln(r)}{k_{total}} $$

Where $r \in (0, 1]$ is a uniform random number.

10.3 Molecular Dynamics

Newton's equations of motion:

$$ m_i \frac{d^2 \mathbf{r}_i}{dt^2} = - abla_i U(\mathbf{r}_1, \mathbf{r}_2, \ldots, \mathbf{r}_N) $$

Lennard-Jones potential:

$$ U_{LJ}(r) = 4\varepsilon \left[ \left( \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r} \right)^6 \right] $$

Embedded Atom Method (EAM) for metals:

$$ U = \sum_i F_i(\rho_i) + \frac{1}{2} \sum_{i eq j} \phi_{ij}(r_{ij}) $$

Where $\rho_i = \sum_{j eq i} f_j(r_{ij})$ is the electron density at atom $i$.

11. Uniformity Modeling

11.1 Wafer-Scale Thickness Distribution (Sputtering)

For a circular magnetron target:

$$ t(r) = \int_{target} \frac{Y \cdot J_{ion} \cdot \cos\theta_t \cdot \cos\theta_w}{\pi R^2} \, dA $$

Where:

• $t(r)$ — thickness at radial position $r$
• $\theta_t$ — emission angle from target
• $\theta_w$ — incidence angle at wafer

11.2 Uniformity Metrics

Within-Wafer Uniformity (WIW):

$$ \sigma_{WIW} = \frac{1}{\bar{t}} \sqrt{\frac{1}{N} \sum_{i=1}^{N} (t_i - \bar{t})^2} \times 100\% $$

Wafer-to-Wafer Uniformity (WTW):

$$ \sigma_{WTW} = \frac{1}{\bar{t}_{avg}} \sqrt{\frac{1}{M} \sum_{j=1}^{M} (\bar{t}_j - \bar{t}_{avg})^2} \times 100\% $$

Target specifications:

• $\sigma_{WIW} < 1\%$ for advanced nodes (≤7 nm)
• $\sigma_{WTW} < 0.5\%$ for high-volume manufacturing

12. Virtual Metrology and Statistical Models

12.1 Gaussian Process Regression (GPR)

$$ f(\mathbf{x}) \sim \mathcal{GP}(m(\mathbf{x}), k(\mathbf{x}, \mathbf{x}')) $$

Squared exponential (RBF) kernel:

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

Predictive distribution:

$$ f_ | \mathbf{X}, \mathbf{y}, \mathbf{x}_ \sim \mathcal{N}(\bar{f}_, \text{var}(f_)) $$

12.2 Partial Least Squares (PLS)

$$ \mathbf{Y} = \mathbf{X} \mathbf{B} + \mathbf{E} $$

Where:

• $\mathbf{X}$ — process parameter matrix
• $\mathbf{Y}$ — quality outcome matrix
• $\mathbf{B}$ — regression coefficient matrix
• $\mathbf{E}$ — residual matrix

12.3 Principal Component Analysis (PCA)

$$ \mathbf{X} = \mathbf{T} \mathbf{P}^T + \mathbf{E} $$

Hotelling's $T^2$ statistic for fault detection:

$$ T^2 = \sum_{i=1}^{k} \frac{t_i^2}{\lambda_i} $$

13. Process Optimization

13.1 Response Surface Methodology (RSM)

Second-order polynomial model:

$$ y = \beta_0 + \sum_{i=1}^{k} \beta_i x_i + \sum_{i=1}^{k} \beta_{ii} x_i^2 + \sum_{i < j} \beta_{ij} x_i x_j + \varepsilon $$

13.2 Constrained Optimization

$$ \min_{\mathbf{x}} f(\mathbf{x}) \quad \text{subject to} \quad g_i(\mathbf{x}) \leq 0, \quad h_j(\mathbf{x}) = 0 $$

Example constraints:

• $g_1$: Non-uniformity ≤ 3%
• $g_2$: Resistivity within spec
• $g_3$: Throughput ≥ target
• $h_1$: Total film thickness = target

13.3 Pareto Multi-Objective Optimization

$$ \min_{\mathbf{x}} \left[ f_1(\mathbf{x}), f_2(\mathbf{x}), \ldots, f_m(\mathbf{x}) \right] $$

Common trade-offs:

• Uniformity vs. throughput
• Film quality vs. cost
• Conformality vs. deposition rate

14. Mathematical Toolkit

15. Physical Constants