Metal Deposition

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:

ProcessPhysics DomainKey Mathematics
PVD (Sputtering)Ballistic transport, plasma physicsBoltzmann transport, Monte Carlo
CVD/PECVDGas-phase transport, surface reactionsNavier-Stokes, reaction-diffusion
ALDSelf-limiting surface chemistrySite-balance kinetics
Electroplating (ECD)Electrochemistry, mass transportButler-Volmer, Nernst-Planck

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 \nabla C = D \nabla^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 \nabla \mathbf{v} \right) = -\nabla p + \mu \nabla^2 \mathbf{v} $$

With continuity equation:

$$ \nabla \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 \nabla^2 h + F + \eta(\mathbf{x}, t) $$

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

$$ \frac{\partial h}{\partial t} = u \nabla^2 h + \frac{\lambda}{2} |\nabla 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 \nabla f + \frac{e \mathbf{E}}{m_e} \cdot \nabla_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

MethodApplicationCharacteristics
Q-LearningDiscrete parameter optimizationModel-free, tabular
Deep Q-Network (DQN)High-dimensional state spacesNeural network approximation
Policy GradientContinuous process controlDirect policy optimization
Actor-Critic (A2C/PPO)Complex control tasksCombined value and policy
Model-Based RLPhysics-informed controlSample efficient

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 \nabla C_i - \frac{z_i F D_i}{RT} C_i \nabla \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

Metal Deposition Technologies in Sub-3nm Interconnects PVD Sputtering, CVD, ALD, and Electrochemical Plating (ECP) for Cu / Ru / Mo 1. PVD (Sputtering) Metal Target (Ta/Ti) Physical Momentum Transfer Line-of-Sight Deposit High Deposition Rate Poor Step Coverage Used for: Seed Layer, Adhesion / Liners Ionized PVD (iPVD) 2. CVD Metal Chemical Gas Reaction Conformal Deposition Medium-High Temp Good Aspect Ratio Fill Used for: Tungsten (W) Plugs & Co contacts MOCVD / PECVD 3. ALD Metal Self-Limiting Monolayers 100% Conformal Fill Atomic-Level Precision Slow Deposition Rate Used for: Ru / Mo Liners & Ultra-thin Barriers Sub-nm Control 4. ECP (Plating) Superfill Cu Electrochemical Plating Bottom-Up Superfill Void-Free Trench Fill Accelerator & Suppressor Used for: Dual Damascene Copper (Cu) Interconnect Chemical Additives Co-Optimization of Barrier, Seed, and Bulk Metal Deposition Sequence for Modern Dual Damascene BEOL

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} = -\nabla_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

DomainKey EquationsApplication
TransportNavier-Stokes, Convection-DiffusionGas flow, precursor delivery
KineticsArrhenius, Langmuir-HinshelwoodReaction rates
Surface EvolutionKPZ, Level-set, Edwards-WilkinsonFilm morphology
PlasmaBoltzmann, Child-LangmuirIon/electron dynamics
ElectrochemistryButler-Volmer, Nernst-PlanckCopper plating
ControlBellman, MDP, RL algorithmsRecipe optimization
StatisticsGPR, PLS, PCAVirtual metrology
MultiscaleMD, kMC, ContinuumIntegrated simulation

15. Physical Constants

ConstantSymbolValueUnits
Boltzmann constant$k_B$$1.38 \times 10^{-23}$J/K
Gas constant$R$$8.314$J/(mol$\cdot$K)
Faraday constant$F$$96,485$C/mol
Elementary charge$e$$1.60 \times 10^{-19}$C
Vacuum permittivity$\varepsilon_0$$8.85 \times 10^{-12}$F/m
Avogadro's number$N_A$$6.02 \times 10^{23}$mol⁻¹
Electron mass$m_e$$9.11 \times 10^{-31}$kg

Take metal deposition further

Ask the copilot about this term, or have our engineers assess it against your process.