process monitoring
**Semiconductor Manufacturing Process Parameters Monitoring: Mathematical Modeling**
**1. The Fundamental Challenge**
Modern semiconductor fabrication involves 500–1000+ sequential process steps, each with dozens of parameters requiring nanometer-scale precision.
**Key Process Types and Parameters**
- **Lithography**: exposure dose, focus, overlay alignment, resist thickness
- **Etching (dry/wet)**: etch rate, selectivity, uniformity, plasma parameters (power, pressure, gas flows)
- **Deposition (CVD, PVD, ALD)**: deposition rate, film thickness, uniformity, stress, composition
- **CMP (Chemical Mechanical Polishing)**: removal rate, within-wafer non-uniformity, dishing, erosion
- **Implantation**: dose, energy, angle, uniformity
- **Thermal processes**: temperature uniformity, ramp rates, time
**2. Statistical Process Control (SPC) — The Foundation**
**2.1 Univariate Control Charts**
For a process parameter $X$ with samples $x_1, x_2, \ldots, x_n$:
**Sample Mean:**
$$
\bar{x} = \frac{1}{n}\sum_{i=1}^{n} x_i
$$
**Sample Standard Deviation:**
$$
\sigma = \sqrt{\frac{1}{n-1}\sum_{i=1}^{n}(x_i - \bar{x})^2}
$$
**Control Limits (3-sigma):**
$$
\text{UCL} = \bar{x} + 3\sigma
$$
$$
\text{LCL} = \bar{x} - 3\sigma
$$
**2.2 Process Capability Indices**
These quantify how well a process meets specifications:
- **$C_p$ (Potential Capability):**
$$
C_p = \frac{USL - LSL}{6\sigma}
$$
- **$C_{pk}$ (Actual Capability)** — accounts for centering:
$$
C_{pk} = \min\left[\frac{USL - \mu}{3\sigma}, \frac{\mu - LSL}{3\sigma}\right]
$$
- **$C_{pm}$ (Taguchi Index)** — penalizes deviation from target $T$:
$$
C_{pm} = \frac{C_p}{\sqrt{1 + \left(\frac{\mu - T}{\sigma}\right)^2}}
$$
Semiconductor fabs typically require $C_{pk} \geq 1.67$, corresponding to defect rates below ~1 ppm.
**3. Multivariate Statistical Monitoring**
Since process parameters are highly correlated, univariate methods miss interaction effects.
**3.1 Principal Component Analysis (PCA)**
Given data matrix $\mathbf{X}$ ($n$ samples × $p$ variables), centered:
1. **Compute covariance matrix:**
$$
\mathbf{S} = \frac{1}{n-1}\mathbf{X}^T\mathbf{X}
$$
2. **Eigendecomposition:**
$$
\mathbf{S} = \mathbf{V}\mathbf{\Lambda}\mathbf{V}^T
$$
3. **Project to principal components:**
$$
\mathbf{T} = \mathbf{X}\mathbf{V}
$$
**3.2 Monitoring Statistics**
**Hotelling's $T^2$ Statistic**
Captures variation **within** the PCA model:
$$
T^2 = \sum_{i=1}^{k} \frac{t_i^2}{\lambda_i}
$$
where $k$ is the number of retained components. Under normal operation, $T^2$ follows a scaled F-distribution.
**Q-Statistic (Squared Prediction Error)**
Captures variation **outside** the model:
$$
Q = \sum_{j=1}^{p}(x_j - \hat{x}_j)^2 = \|\mathbf{x} - \mathbf{x}\mathbf{V}_k\mathbf{V}_k^T\|^2
$$
> Often more sensitive to novel faults than $T^2$.
**3.3 Partial Least Squares (PLS)**
When relating process inputs $\mathbf{X}$ to quality outputs $\mathbf{Y}$:
$$
\mathbf{Y} = \mathbf{X}\mathbf{B} + \mathbf{E}
$$
PLS finds latent variables that maximize covariance between $\mathbf{X}$ and $\mathbf{Y}$, providing both monitoring capability and a predictive model.
**4. Virtual Metrology (VM) Models**
Virtual metrology predicts physical measurement outcomes from process sensor data, enabling 100% wafer coverage without costly measurements.
**4.1 Linear Models**
For process parameters $\mathbf{x} \in \mathbb{R}^p$ and metrology target $y$:
- **Ordinary Least Squares (OLS):**
$$
\hat{\boldsymbol{\beta}} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}
$$
- **Ridge Regression** ($L_2$ regularization for collinearity):
$$
\hat{\boldsymbol{\beta}} = (\mathbf{X}^T\mathbf{X} + \lambda\mathbf{I})^{-1}\mathbf{X}^T\mathbf{y}
$$
- **LASSO** ($L_1$ regularization for sparsity/feature selection):
$$
\min_{\boldsymbol{\beta}} \|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|^2 + \lambda\|\boldsymbol{\beta}\|_1
$$
**4.2 Nonlinear Models**
**Gaussian Process Regression (GPR)**
$$
y \sim \mathcal{GP}(m(\mathbf{x}), k(\mathbf{x}, \mathbf{x}'))
$$
**Posterior predictive distribution:**
- **Mean:**
$$
\mu_* = \mathbf{K}_*^T(\mathbf{K} + \sigma_n^2\mathbf{I})^{-1}\mathbf{y}
$$
- **Variance:**
$$
\sigma_*^2 = K_{**} - \mathbf{K}_*^T(\mathbf{K} + \sigma_n^2\mathbf{I})^{-1}\mathbf{K}_*
$$
GPs provide uncertainty quantification — critical for knowing when to trigger actual metrology.
**Support Vector Regression (SVR)**
$$
\min \frac{1}{2}\|\mathbf{w}\|^2 + C\sum_i(\xi_i + \xi_i^*)
$$
Subject to $\epsilon$-insensitive tube constraints. Kernel trick enables nonlinear modeling.
**Neural Networks**
- **MLPs**: Multi-layer perceptrons for general function approximation
- **CNNs**: Convolutional neural networks for wafer map pattern recognition
- **LSTMs**: Long Short-Term Memory networks for time-series FDC traces
**5. Run-to-Run (R2R) Control**
R2R control adjusts recipe setpoints between wafers/lots to compensate for drift and disturbances.
**5.1 EWMA Controller**
For a process with model $y = a_0 + a_1 u + \epsilon$:
**Prediction update:**
$$
\hat{y}_{k+1} = \lambda y_k + (1-\lambda)\hat{y}_k
$$
**Control action:**
$$
u_{k+1} = \frac{T - \hat{y}_{k+1} + a_0}{a_1}
$$
where:
- $T$ is the target
- $\lambda \in (0,1)$ is the smoothing weight
**5.2 Double EWMA (for Linear Drift)**
When process drifts linearly:
$$
\hat{y}_{k+1} = a_k + b_k
$$
$$
a_k = \lambda y_k + (1-\lambda)(a_{k-1} + b_{k-1})
$$
$$
b_k = \gamma(a_k - a_{k-1}) + (1-\gamma)b_{k-1}
$$
**5.3 State-Space Formulation**
More general framework:
**State equation:**
$$
\mathbf{x}_{k+1} = \mathbf{A}\mathbf{x}_k + \mathbf{B}\mathbf{u}_k + \mathbf{w}_k
$$
**Observation equation:**
$$
\mathbf{y}_k = \mathbf{C}\mathbf{x}_k + \mathbf{D}\mathbf{u}_k + \mathbf{v}_k
$$
Use **Kalman filtering** for state estimation and **LQR/MPC** for optimal control.
**5.4 Model Predictive Control (MPC)**
**Objective function:**
$$
\min \sum_{i=1}^{N} \|\mathbf{y}_{k+i} - \mathbf{r}_{k+i}\|_\mathbf{Q}^2 + \sum_{j=0}^{N-1}\|\Delta\mathbf{u}_{k+j}\|_\mathbf{R}^2
$$
subject to process model and operational constraints.
> MPC handles multivariable systems with constraints naturally.
**6. Fault Detection and Classification (FDC)**
**6.1 Detection Methods**
**Mahalanobis Distance**
$$
D^2 = (\mathbf{x} - \boldsymbol{\mu})^T\mathbf{S}^{-1}(\mathbf{x} - \boldsymbol{\mu})
$$
Follows $\chi^2$ distribution under multivariate normality.
**Other Detection Methods**
- **One-Class SVM**: Learn boundary of normal operation
- **Autoencoders**: Detect anomalies via reconstruction error
**6.2 Classification Features**
For trace data (time-series from sensors), extract features:
- **Statistical moments**: mean, variance, skewness, kurtosis
- **Frequency domain**: FFT coefficients, spectral power
- **Wavelet coefficients**: Multi-resolution analysis
- **DTW distances**: Dynamic Time Warping to reference signatures
**6.3 Classification Algorithms**
- Support Vector Machines (SVM)
- Random Forest
- CNNs for pattern recognition on wafer maps
- Gradient Boosting (XGBoost, LightGBM)
**7. Spatial Modeling (Within-Wafer Variation)**
Systematic spatial patterns require explicit modeling.
**7.1 Polynomial Basis Expansion**
**Zernike Polynomials (common in lithography)**
$$
z(\rho, \theta) = \sum_{n,m} Z_n^m(\rho, \theta)
$$
These form an orthogonal basis on the unit disk, capturing radial and azimuthal variation.
**7.2 Gaussian Process Spatial Models**
$$
y(\mathbf{s}) \sim \mathcal{GP}(\mu(\mathbf{s}), k(\mathbf{s}, \mathbf{s}'))
$$
**Common Covariance Kernels**
- **Squared Exponential (RBF):**
$$
k(\mathbf{s}, \mathbf{s}') = \sigma^2 \exp\left(-\frac{\|\mathbf{s} - \mathbf{s}'\|^2}{2\ell^2}\right)
$$
- **Matérn** (more flexible smoothness):
$$
k(r) = \sigma^2 \frac{2^{1-
u}}{\Gamma(
u)}\left(\frac{\sqrt{2
u}r}{\ell}\right)^
u K_
u\left(\frac{\sqrt{2
u}r}{\ell}\right)
$$
where $K_
u$ is the modified Bessel function of the second kind.
**8. Dynamic/Time-Series Modeling**
For plasma processes, endpoint detection, and transient behavior.
**8.1 Autoregressive Models**
**AR(p) model:**
$$
x_t = \sum_{i=1}^{p} \phi_i x_{t-i} + \epsilon_t
$$
ARIMA extends this to non-stationary series.
**8.2 Dynamic PCA**
Augment data with time-lagged values:
$$
\tilde{\mathbf{X}} = [\mathbf{X}(t), \mathbf{X}(t-1), \ldots, \mathbf{X}(t-l)]
$$
Then apply standard PCA to capture temporal dynamics.
**8.3 Deep Sequence Models**
**LSTM Networks**
Gating mechanisms:
- **Forget gate:** $f_t = \sigma(W_f \cdot [h_{t-1}, x_t] + b_f)$
- **Input gate:** $i_t = \sigma(W_i \cdot [h_{t-1}, x_t] + b_i)$
- **Output gate:** $o_t = \sigma(W_o \cdot [h_{t-1}, x_t] + b_o)$
**Cell state update:**
$$
c_t = f_t \odot c_{t-1} + i_t \odot \tilde{c}_t
$$
**Hidden state:**
$$
h_t = o_t \odot \tanh(c_t)
$$
**9. Model Maintenance and Adaptation**
Semiconductor processes drift — models must adapt.
**9.1 Drift Detection Methods**
**CUSUM (Cumulative Sum)**
$$
S_k = \max(0, S_{k-1} + (x_k - \mu_0) - k)
$$
Signal when $S_k$ exceeds threshold.
**Page-Hinkley Test**
$$
m_k = \sum_{i=1}^{k}(x_i - \bar{x}_k - \delta)
$$
$$
M_k = \max_{i \leq k} m_i
$$
Alarm when $M_k - m_k > \lambda$.
**ADWIN (Adaptive Windowing)**
Automatically detects distribution changes and adjusts window size.
**9.2 Online Model Updating**
**Recursive Least Squares (RLS)**
$$
\hat{\boldsymbol{\beta}}_k = \hat{\boldsymbol{\beta}}_{k-1} + \mathbf{K}_k(y_k - \mathbf{x}_k^T\hat{\boldsymbol{\beta}}_{k-1})
$$
where $\mathbf{K}_k$ is the gain matrix updated via the Riccati equation:
$$
\mathbf{K}_k = \frac{\mathbf{P}_{k-1}\mathbf{x}_k}{\lambda + \mathbf{x}_k^T\mathbf{P}_{k-1}\mathbf{x}_k}
$$
$$
\mathbf{P}_k = \frac{1}{\lambda}(\mathbf{P}_{k-1} - \mathbf{K}_k\mathbf{x}_k^T\mathbf{P}_{k-1})
$$
**Just-in-Time (JIT) Learning**
Build local models around each new prediction point using nearest historical samples.
**10. Integrated Framework**
A complete monitoring system layers these methods:
| Layer | Methods | Purpose |
|-------|---------|---------|
| **Preprocessing** | Cleaning, synchronization, normalization | Data quality |
| **Feature Engineering** | Domain features, wavelets, PCA | Dimensionality management |
| **Monitoring** | $T^2$, Q-statistic, control charts | Detect out-of-control states |
| **Virtual Metrology** | PLS, GPR, neural networks | Predict quality without measurement |
| **FDC** | Classification models | Diagnose fault root causes |
| **Control** | R2R, MPC | Compensate for drift/disturbances |
| **Adaptation** | Online learning, drift detection | Maintain model validity |
**11. Key Mathematical Challenges**
1. **High dimensionality** — hundreds of sensors, requiring regularization and dimension reduction
2. **Collinearity** — process variables are physically coupled
3. **Non-stationarity** — drift, maintenance events, recipe changes
4. **Small sample sizes** — new recipes have limited historical data (transfer learning, Bayesian methods help)
5. **Real-time constraints** — decisions needed in seconds
6. **Rare events** — faults are infrequent, creating class imbalance
**12. Key Equations**
**Process Capability**
$$
C_{pk} = \min\left[\frac{USL - \mu}{3\sigma}, \frac{\mu - LSL}{3\sigma}\right]
$$
**Multivariate Monitoring**
$$
T^2 = \sum_{i=1}^{k} \frac{t_i^2}{\lambda_i}, \quad Q = \|\mathbf{x} - \hat{\mathbf{x}}\|^2
$$
**Virtual Metrology (Ridge Regression)**
$$
\hat{\boldsymbol{\beta}} = (\mathbf{X}^T\mathbf{X} + \lambda\mathbf{I})^{-1}\mathbf{X}^T\mathbf{y}
$$
**EWMA Control**
$$
\hat{y}_{k+1} = \lambda y_k + (1-\lambda)\hat{y}_k
$$
**Mahalanobis Distance**
$$
D^2 = (\mathbf{x} - \boldsymbol{\mu})^T\mathbf{S}^{-1}(\mathbf{x} - \boldsymbol{\mu})
$$