multi physics coupling
**Semiconductor Manufacturing Process: Multi-Physics Coupling & Mathematical Modeling**
**1. Overview: Why Multi-Physics Coupling Matters**
Semiconductor fabrication involves hundreds of process steps where multiple physical phenomena occur simultaneously and interact nonlinearly. At the 3nm node and below, these couplings become critical—small perturbations propagate across physics domains, affecting yield, uniformity, and device performance.
**2. Key Processes and Their Coupled Physics**
**2.1 Plasma Etching (RIE, ICP, CCP)**
**Coupled domains:**
- Electromagnetics (RF field, power deposition)
- Plasma kinetics (electron/ion transport, sheath dynamics)
- Neutral gas fluid dynamics
- Gas-phase and surface chemistry
- Heat transfer
- Feature-scale transport and profile evolution
**Coupling chain:**
```svg
```
**2.2 Chemical Vapor Deposition (CVD/ALD)**
**Coupled domains:**
- Fluid dynamics (often rarefied/transitional flow)
- Heat transfer (convection, conduction, radiation)
- Multi-component mass transfer
- Gas-phase and surface reaction kinetics
- Film stress evolution
**2.3 Thermal Processing (RTP, Annealing)**
**Coupled domains:**
- Radiation heat transfer
- Solid-state diffusion (dopants)
- Defect kinetics
- Thermo-mechanical stress (slip, warpage)
**2.4 EUV Lithography**
**Coupled domains:**
- Wave optics and diffraction
- Photochemistry in resist
- Stochastic photon/electron effects
- Mask/wafer thermal-mechanical deformation
**3. Mathematical Framework: Governing Equations**
**3.1 Electromagnetics (Plasma Systems)**
For RF-driven plasma, the **time-harmonic Maxwell's equations**:
$$
abla \times \left(\mu_r^{-1}
abla \times \mathbf{E}\right) - k_0^2 \epsilon_r \mathbf{E} = -j\omega\mu_0 \mathbf{J}_{ext}
$$
The **plasma permittivity** encodes the coupling to electron density:
$$
\epsilon_r = 1 - \frac{\omega_{pe}^2}{\omega(\omega + j
u_m)}
$$
Where the **plasma frequency** is:
$$
\omega_{pe} = \sqrt{\frac{n_e e^2}{m_e \epsilon_0}}
$$
**Key parameters:**
- $n_e$ — electron density
- $e$ — electron charge
- $m_e$ — electron mass
- $\epsilon_0$ — permittivity of free space
- $
u_m$ — electron-neutral collision frequency
- $\omega$ — angular frequency of RF excitation
> **Note:** This creates a **strong nonlinear coupling**: the EM field depends on plasma density, which in turn depends on power absorption from the EM field.
**3.2 Plasma Transport (Drift-Diffusion Approximation)**
**Electron continuity equation:**
$$
\frac{\partial n_e}{\partial t} +
abla \cdot \boldsymbol{\Gamma}_e = S_e
$$
**Electron flux:**
$$
\boldsymbol{\Gamma}_e = -\mu_e n_e \mathbf{E} - D_e
abla n_e
$$
**Electron energy density equation:**
$$
\frac{\partial n_\epsilon}{\partial t} +
abla \cdot \boldsymbol{\Gamma}_\epsilon + \mathbf{E} \cdot \boldsymbol{\Gamma}_e = S_\epsilon - \sum_j \varepsilon_j R_j
$$
**Where:**
- $n_e$ — electron density
- $\boldsymbol{\Gamma}_e$ — electron flux vector
- $\mu_e$ — electron mobility
- $D_e$ — electron diffusion coefficient
- $S_e$ — electron source term (ionization, attachment, recombination)
- $n_\epsilon$ — electron energy density
- $\varepsilon_j$ — energy loss per reaction $j$
- $R_j$ — reaction rate for process $j$
**Ion transport** (for multiple species $i$):
$$
\frac{\partial n_i}{\partial t} +
abla \cdot \boldsymbol{\Gamma}_i = S_i
$$
**3.3 Neutral Gas Flow (Navier-Stokes Equations)**
**Continuity equation:**
$$
\frac{\partial \rho}{\partial t} +
abla \cdot (\rho \mathbf{u}) = 0
$$
**Momentum equation:**
$$
\rho \frac{D\mathbf{u}}{Dt} = -
abla p +
abla \cdot \boldsymbol{\tau} + \mathbf{F}_{body}
$$
**Where:**
- $\rho$ — gas density
- $\mathbf{u}$ — velocity vector
- $p$ — pressure
- $\boldsymbol{\tau}$ — viscous stress tensor
- $\mathbf{F}_{body}$ — body forces
**Low-pressure corrections (Knudsen effects):**
At low pressures where Knudsen number $Kn = \lambda/L > 0.01$, slip boundary conditions are required:
$$
u_{slip} = \frac{2-\sigma}{\sigma} \lambda \left.\frac{\partial u}{\partial n}\right|_{wall}
$$
Where:
- $\lambda$ — mean free path
- $L$ — characteristic length
- $\sigma$ — tangential momentum accommodation coefficient
**3.4 Species Transport and Chemistry**
**Convection-diffusion-reaction equation:**
$$
\frac{\partial c_k}{\partial t} +
abla \cdot (c_k \mathbf{u}) =
abla \cdot (D_k
abla c_k) + R_k
$$
**Gas-phase reaction rates:**
$$
R_k = \sum_j
u_{kj} \, k_j(T) \prod_l c_l^{a_{lj}}
$$
**Where:**
- $c_k$ — concentration of species $k$
- $D_k$ — diffusion coefficient
- $R_k$ — net production rate
- $
u_{kj}$ — stoichiometric coefficient
- $k_j(T)$ — temperature-dependent rate constant
- $a_{lj}$ — reaction order
**Surface reactions (Langmuir-Hinshelwood kinetics):**
$$
r_s = k_s \theta_A \theta_B
$$
**Surface coverage:**
$$
\theta_i = \frac{K_i c_i}{1 + \sum_j K_j c_j}
$$
**3.5 Heat Transfer**
**Energy equation:**
$$
\rho c_p \frac{\partial T}{\partial t} + \rho c_p \mathbf{u} \cdot
abla T =
abla \cdot (k
abla T) + Q
$$
**Heat sources in plasma systems:**
$$
Q = Q_{Joule} + Q_{ion} + Q_{reaction} + Q_{radiation}
$$
**Joule heating (time-averaged):**
$$
Q_{Joule} = \frac{1}{2} \text{Re}(\mathbf{J}^* \cdot \mathbf{E})
$$
**Where:**
- $\rho$ — density
- $c_p$ — specific heat capacity
- $k$ — thermal conductivity
- $Q$ — volumetric heat source
- $\mathbf{J}^*$ — complex conjugate of current density
**3.6 Solid Mechanics (Film Stress)**
**Equilibrium equation:**
$$
abla \cdot \boldsymbol{\sigma} = 0
$$
**Constitutive relation with thermal strain:**
$$
\boldsymbol{\sigma} = \mathbf{C} : (\boldsymbol{\epsilon} - \boldsymbol{\epsilon}_{th} - \boldsymbol{\epsilon}_{intrinsic})
$$
**Thermal strain tensor:**
$$
\boldsymbol{\epsilon}_{th} = \alpha(T - T_0)\mathbf{I}
$$
**Where:**
- $\boldsymbol{\sigma}$ — stress tensor
- $\mathbf{C}$ — stiffness tensor
- $\boldsymbol{\epsilon}$ — total strain tensor
- $\alpha$ — coefficient of thermal expansion
- $T_0$ — reference temperature
- $\mathbf{I}$ — identity tensor
**Stoney equation** (wafer curvature from film stress):
$$
\sigma_f = \frac{E_s h_s^2}{6(1-
u_s)h_f}\kappa
$$
**Where:**
- $\sigma_f$ — film stress
- $E_s$ — substrate Young's modulus
- $
u_s$ — substrate Poisson's ratio
- $h_s$ — substrate thickness
- $h_f$ — film thickness
- $\kappa$ — wafer curvature
**4. Feature-Scale Modeling**
At the nanometer scale within etched features, continuum assumptions break down.
**4.1 Profile Evolution (Level Set Method)**
The etch front $\phi(\mathbf{x},t) = 0$ evolves according to:
$$
\frac{\partial \phi}{\partial t} + V_n |
abla \phi| = 0
$$
**Local etch rate** depends on coupled physics:
$$
V_n = \Gamma_{ion}(E,\theta) \cdot Y_{phys}(E,\theta) + \Gamma_{rad} \cdot Y_{chem}(T) + \Gamma_{ion} \cdot \Gamma_{rad} \cdot Y_{synergy}
$$
**Where:**
- $\phi$ — level set function (zero at interface)
- $V_n$ — normal velocity of interface
- $\Gamma_{ion}$ — ion flux (from sheath model)
- $\Gamma_{rad}$ — radical flux (from feature-scale transport)
- $Y_{phys}$ — physical sputtering yield
- $Y_{chem}$ — chemical etch yield
- $Y_{synergy}$ — ion-enhanced chemical yield
- $\theta$ — local incidence angle
- $E$ — ion energy
**4.2 Feature-Scale Transport**
Within high-aspect-ratio features, **Knudsen diffusion** dominates:
$$
D_{Kn} = \frac{d}{3}\sqrt{\frac{8k_BT}{\pi m}}
$$
**Where:**
- $d$ — feature diameter/width
- $k_B$ — Boltzmann constant
- $T$ — temperature
- $m$ — molecular mass
**View factor calculations** for flux at the bottom of features:
$$
\Gamma_{bottom} = \Gamma_{top} \cdot \int_{\Omega} f(\theta) \cos\theta \, d\Omega
$$
**4.3 Ion Angular and Energy Distribution**
At the sheath-feature interface:
$$
f(E, \theta) = f_E(E) \cdot f_\theta(\theta)
$$
**Angular distribution** (from sheath collisionality):
$$
f_\theta(\theta) \propto \cos^n(\theta) \exp\left(-\frac{\theta^2}{2\sigma_\theta^2}\right)
$$
**Where:**
- $f_E(E)$ — ion energy distribution function
- $f_\theta(\theta)$ — ion angular distribution function
- $n$ — exponent (depends on sheath collisionality)
- $\sigma_\theta$ — angular spread parameter
**5. Multi-Scale Coupling Strategy**
```svg
```
**Scale hierarchy:**
1. **Reactor scale (cm–m)**
- Continuum fluid dynamics
- Maxwell's equations for EM fields
- Drift-diffusion for charged species
- Numerical methods: FEM, FVM
2. **Feature scale (nm–μm)**
- Knudsen transport in high-aspect-ratio structures
- Direct Simulation Monte Carlo (DSMC)
- Level set methods for profile evolution
3. **Atomic scale (Å–nm)**
- Density Functional Theory (DFT) for reaction barriers
- Molecular Dynamics (MD) for sputtering yields
- Kinetic Monte Carlo (KMC) for surface evolution
**6. Coupled System Structure**
The full system can be written abstractly as:
$$
\mathbf{M}(\mathbf{u})\frac{\partial \mathbf{u}}{\partial t} = \mathbf{F}(\mathbf{u},
abla\mathbf{u},
abla^2\mathbf{u}, t)
$$
**State vector:**
$$
\mathbf{u} = \begin{bmatrix} n_e \\ n_\epsilon \\ n_{i,k} \\ c_j \\ T \\ \mathbf{E} \\ \mathbf{u}_{gas} \\ p \\ \boldsymbol{\sigma} \\ \phi_{profile} \\ \vdots \end{bmatrix}
$$
**Jacobian structure reveals coupling:**
$$
\mathbf{J} = \frac{\partial \mathbf{F}}{\partial \mathbf{u}} = \begin{pmatrix}
J_{ee} & J_{e\epsilon} & J_{ei} & J_{ec} & \cdots \\
J_{\epsilon e} & J_{\epsilon\epsilon} & J_{\epsilon i} & & \\
J_{ie} & J_{i\epsilon} & J_{ii} & & \\
J_{ce} & & & J_{cc} & \\
\vdots & & & & \ddots
\end{pmatrix}
$$
**Off-diagonal blocks** represent inter-physics coupling strengths.
**7. Numerical Solution Strategies**
**7.1 Coupling Approaches**
**Monolithic (fully coupled):**
- Solve all physics simultaneously
- Newton iteration on full Jacobian
- Robust but computationally expensive
- Required for strongly coupled physics (plasma + EM)
**Partitioned (sequential):**
- Solve each physics domain separately
- Iterate between domains until convergence
- More efficient for weakly coupled physics
- Risk of convergence issues
**Hybrid approach:**
- Group strongly coupled physics into blocks
- Sequential coupling between blocks
**7.2 Spatial Discretization**
**Finite Element Method (FEM)** — weak form for species transport:
$$
\int_\Omega w \frac{\partial c}{\partial t} \, d\Omega + \int_\Omega w (\mathbf{u} \cdot
abla c) \, d\Omega + \int_\Omega
abla w \cdot (D
abla c) \, d\Omega = \int_\Omega w R \, d\Omega
$$
**SUPG Stabilization** for convection-dominated problems:
$$
w \rightarrow w + \tau_{SUPG} \, \mathbf{u} \cdot
abla w
$$
**Where:**
- $w$ — test function
- $c$ — concentration field
- $\tau_{SUPG}$ — stabilization parameter
**7.3 Time Integration**
**Stiff systems** require implicit methods:
- **BDF** (Backward Differentiation Formulas)
- **ESDIRK** (Explicit Singly Diagonally Implicit Runge-Kutta)
**Operator splitting** for multi-physics:
$$
\mathbf{u}^{n+1} = \mathcal{L}_1(\Delta t) \circ \mathcal{L}_2(\Delta t) \circ \mathcal{L}_3(\Delta t) \, \mathbf{u}^n
$$
**Where:**
- $\mathcal{L}_i$ — solution operator for physics domain $i$
- $\Delta t$ — time step
- $\circ$ — composition of operators
**8. Specific Application: ICP Etch Model**
**Complete coupled system summary:**
| Physics Domain | Governing Equations | Key Coupling Variables |
|----------------|---------------------|------------------------|
| EM (inductive) | $
abla \times (
abla \times \mathbf{E}) + k^2\epsilon_p \mathbf{E} = 0$ | $n_e \rightarrow \epsilon_p$ |
| Electron transport | $
abla \cdot \Gamma_e = S_e$ | $\mathbf{E}_{dc}, n_e, T_e$ |
| Electron energy | $
abla \cdot \Gamma_\epsilon = Q_{EM} - Q_{loss}$ | $T_e \rightarrow$ rate coefficients |
| Ion transport | $
abla \cdot \Gamma_i = S_i$ | $n_e, \mathbf{E}_{dc}$ |
| Neutral chemistry | $
abla \cdot (c_k \mathbf{u} - D_k
abla c_k) = R_k$ | $T_e \rightarrow k_{diss}$ |
| Gas flow | Navier-Stokes | $T_{gas}$ |
| Heat transfer | $
abla \cdot (k
abla T) + Q = 0$ | $Q_{plasma}$ |
| Sheath | Child-Langmuir / PIC | $n_e, T_e, V_{dc}$ |
| Feature transport | Knudsen + angular | $\Gamma_{ion}, \Gamma_{rad}$ from reactor |
| Profile evolution | Level set | $V_n$ from surface kinetics |
**9. EUV Lithography: Stochastic Multi-Physics**
At EUV wavelength (13.5 nm), photon shot noise becomes significant.
**9.1 Aerial Image Formation**
$$
I(\mathbf{r}) = \left|\mathcal{F}^{-1}\left[\tilde{M}(\mathbf{f}) \cdot H(\mathbf{f})\right]\right|^2
$$
**Where:**
- $I(\mathbf{r})$ — intensity at position $\mathbf{r}$
- $\tilde{M}(\mathbf{f})$ — mask spectrum (Fourier transform of mask pattern)
- $H(\mathbf{f})$ — pupil function (includes aberrations, partial coherence)
- $\mathcal{F}^{-1}$ — inverse Fourier transform
**9.2 Photon Statistics**
$$
N \sim \text{Poisson}(\bar{N})
$$
$$
\sigma_N = \sqrt{\bar{N}}
$$
**Where:**
- $N$ — number of photons absorbed
- $\bar{N}$ — expected number of photons
- $\sigma_N$ — standard deviation (shot noise)
**9.3 Resist Exposure (Stochastic Dill Model)**
$$
\frac{\partial [PAG]}{\partial t} = -C \cdot I \cdot [PAG] + \xi(t)
$$
**Where:**
- $[PAG]$ — photoactive compound concentration
- $C$ — exposure rate constant
- $I$ — local intensity
- $\xi(t)$ — stochastic noise term
**9.4 Line Edge Roughness (LER)**
$$
\sigma_{LER} \propto \sqrt{\frac{1}{\text{dose}}} \cdot \frac{1}{\text{image contrast}}
$$
> **Note:** This requires **Kinetic Monte Carlo** or **Gillespie algorithm** rather than continuum PDEs.
**10. Process Optimization (Inverse Problem)**
**10.1 Problem Formulation**
**Objective:** Minimize profile deviation from target
$$
\min_{\mathbf{p}} J = \int_\Gamma \left|\phi(\mathbf{x}; \mathbf{p}) - \phi_{target}\right|^2 \, d\Gamma
$$
**Subject to physics constraints:**
$$
\mathbf{F}(\mathbf{u}, \mathbf{p}) = 0
$$
**Control parameters** $\mathbf{p}$:
- RF power
- Chamber pressure
- Gas flow rates
- Substrate temperature
- Process time
**10.2 Adjoint Method for Efficient Gradients**
**Gradient computation:**
$$
\frac{dJ}{d\mathbf{p}} = \frac{\partial J}{\partial \mathbf{p}} - \boldsymbol{\lambda}^T \frac{\partial \mathbf{F}}{\partial \mathbf{p}}
$$
**Adjoint equation:**
$$
\left(\frac{\partial \mathbf{F}}{\partial \mathbf{u}}\right)^T \boldsymbol{\lambda} = \left(\frac{\partial J}{\partial \mathbf{u}}\right)^T
$$
**Where:**
- $\boldsymbol{\lambda}$ — adjoint variable (Lagrange multiplier)
- $\mathbf{u}$ — state variables
- $\mathbf{p}$ — control parameters
**11. Emerging Approaches**
**11.1 Physics-Informed Neural Networks (PINNs)**
**Loss function:**
$$
\mathcal{L} = \mathcal{L}_{data} + \lambda \mathcal{L}_{PDE}
$$
**Where:**
- $\mathcal{L}_{data}$ — data fitting loss
- $\mathcal{L}_{PDE}$ — PDE residual loss at collocation points
- $\lambda$ — regularization parameter
**11.2 Digital Twins**
**Key features:**
- Real-time reduced-order models calibrated to equipment sensors
- Combine physics-based models with ML for fast prediction
- Enable predictive maintenance and process control
**11.3 Uncertainty Quantification**
**Methods:**
- **Polynomial Chaos Expansion (PCE)** — for parametric uncertainty propagation
- **Bayesian Inference** — for model calibration with experimental data
- **Monte Carlo Sampling** — for statistical analysis of outputs
**12. Mathematical Structure**
The semiconductor manufacturing multi-physics problem has a characteristic mathematical structure:
1. **Hierarchy of scales** (atomic → feature → reactor)
- Requires multi-scale methods
- Information passing between scales via homogenization
2. **Nonlinear coupling** between physics domains
- Varying coupling strengths
- Both explicit and implicit dependencies
3. **Stiff ODEs/DAEs**
- Disparate time scales (electron dynamics ~ ns, thermal ~ s)
- Requires implicit time integration
4. **Moving boundaries**
- Etch/deposition fronts
- Requires interface tracking (level set, phase field)
5. **Rarefied gas effects**
- At low pressures ($Kn > 0.01$)
- Requires kinetic corrections or DSMC
6. **Stochastic effects**
- At nanometer scales (EUV, atomic-scale roughness)
- Requires Monte Carlo methods
**Key Physical Constants**
| Symbol | Value | Description |
|--------|-------|-------------|
| $e$ | $1.602 \times 10^{-19}$ C | Elementary charge |
| $m_e$ | $9.109 \times 10^{-31}$ kg | Electron mass |
| $\epsilon_0$ | $8.854 \times 10^{-12}$ F/m | Permittivity of free space |
| $\mu_0$ | $4\pi \times 10^{-7}$ H/m | Permeability of free space |
| $k_B$ | $1.381 \times 10^{-23}$ J/K | Boltzmann constant |
| $N_A$ | $6.022 \times 10^{23}$ mol$^{-1}$ | Avogadro's number |
**Common Dimensionless Numbers**
| Number | Definition | Physical Meaning |
|--------|------------|------------------|
| Knudsen ($Kn$) | $\lambda / L$ | Mean free path / characteristic length |
| Reynolds ($Re$) | $\rho u L / \mu$ | Inertia / viscous forces |
| Péclet ($Pe$) | $u L / D$ | Convection / diffusion |
| Damköhler ($Da$) | $k L / u$ | Reaction / convection rate |
| Biot ($Bi$) | $h L / k$ | Surface / bulk heat transfer |