Semiconductor Manufacturing Process: Numerical Methods, Mathematics & Modeling

Keywords: numerical methods, FEM FDM FVM, finite element, finite difference, conjugate gradient, monte carlo, level set, TCAD simulation, computational methods

Semiconductor Manufacturing Process: Numerical Methods, Mathematics & Modeling

A comprehensive guide covering the mathematical foundations, numerical methods, and computational modeling approaches used in semiconductor fabrication processes.

1. Manufacturing Processes and Their Physics

Semiconductor fabrication involves sequential processes, each governed by different physics:

1.1 Lithography

• Optical projection through reduction lens system
• Photoresist chemistry: exposure, bake, development
• Resolution limit: $R = k_1 \frac{\lambda}{NA}$
• Depth of focus: $DOF = k_2 \frac{\lambda}{NA^2}$

1.2 Plasma Etching

• Plasma generation: RF/microwave excitation
• Ion bombardment: directional etching
• Chemical reactions: isotropic component
• Selectivity: differential etch rates between materials

1.3 Chemical Vapor Deposition (CVD)

• Gas-phase transport: convection and diffusion
• Surface reactions: adsorption, reaction, desorption
• Film conformality: step coverage in features
• Temperature dependence: Arrhenius kinetics

1.4 Ion Implantation

• Ion acceleration: keV to MeV energies
• Stopping mechanisms: electronic and nuclear
• Damage formation: vacancy-interstitial pairs
• Channeling effects: crystallographic orientation dependence

2. Core Mathematical Frameworks

2.1 Partial Differential Equations

Nearly every process involves PDEs of different types:

Parabolic (Diffusion/Heat Transport)

• Application: Dopant diffusion, thermal processing, resist chemistry
• Characteristics: Smoothing behavior, infinite propagation speed
• Diffusion coefficient: $D = D_0 \exp\left(-\frac{E_a}{k_B T}\right)$

Elliptic (Steady-State Fields)

• Application: Electrostatics, plasma sheaths, device simulation
• Boundary conditions: Dirichlet, Neumann, or mixed
• Properties: Maximum principle, smoothness

Hyperbolic (Wave Propagation)

• Application: Light propagation in lithography
• Characteristics: Finite propagation speed
• Dispersion: wavelength-dependent phase velocity

2.2 Transport Theory

The Boltzmann transport equation underpins plasma modeling and carrier transport:

Where:

• $f(\mathbf{r}, \mathbf{v}, t)$ = distribution function (6D phase space)
• $\mathbf{F}$ = external force (electric field, etc.)
• RHS = collision integral

Solution approaches:

• Moment methods: Fluid approximations (continuity, momentum, energy)
• Monte Carlo sampling: Stochastic particle tracking
• Deterministic discretization: Spherical harmonics expansion

2.3 Reaction-Diffusion Systems

Coupled species with chemical reactions:

Examples:

• Dopant-defect interactions: Transient enhanced diffusion
• Dopants: $\frac{\partial C_D}{\partial t} =

abla \cdot (D_D abla C_D) + k_{DI} C_D C_I$

• Interstitials: $\frac{\partial C_I}{\partial t} =

abla \cdot (D_I abla C_I) - k_{IV} C_I C_V + G$

• Vacancies: $\frac{\partial C_V}{\partial t} =

abla \cdot (D_V abla C_V) - k_{IV} C_I C_V + G$

• Resist chemistry:
• Photoacid generation: $\frac{\partial [PAG]}{\partial t} = -C \cdot I \cdot [PAG]$
• Acid diffusion: $\frac{\partial [H^+]}{\partial t} = D_{acid}

abla^2 [H^+]$

• Deprotection: $\frac{\partial M}{\partial t} = -k_{amp} [H^+] M$

2.4 Semiconductor Device Equations

The drift-diffusion model for carrier transport:

Current densities:

Einstein relation: $D = \frac{k_B T}{q} \mu$

3. Numerical Methods by Category

3.1 Spatial Discretization

Finite Difference Method (FDM)

Central difference (second derivative):

Forward difference (first derivative):

Characteristics:

• Simple implementation on regular grids
• Truncation error: $O(\Delta x^2)$ for central differences
• Challenges with complex geometries
• Stability requires careful time step selection

Finite Element Method (FEM)

Variational formulation - find $u$ minimizing:

Weak form - find $u \in V$ such that for all $v \in V$:

Implementation steps:

1. Mesh generation: Divide domain into elements (triangles, tetrahedra) 2. Shape functions: Local polynomial basis $N_i(\mathbf{x})$ 3. Assembly: Build global stiffness matrix $\mathbf{K}$ and load vector $\mathbf{f}$ 4. Solution: Solve $\mathbf{K} \mathbf{u} = \mathbf{f}$

Advantages:

• Handles complex geometries naturally
• Systematic error estimation
• Adaptive refinement possible

Finite Volume Method (FVM)

Conservation form:

Discrete form (cell $i$):

Characteristics:

• Conserves quantities exactly by construction
• Natural for fluid dynamics
• Upwinding for convection-dominated problems

3.2 Time Integration

Explicit Methods

Forward Euler:

Runge-Kutta 4th order (RK4):

Where:

• $k_1 = f(t^n, u^n)$
• $k_2 = f(t^n + \frac{\Delta t}{2}, u^n + \frac{\Delta t}{2} k_1)$
• $k_3 = f(t^n + \frac{\Delta t}{2}, u^n + \frac{\Delta t}{2} k_2)$
• $k_4 = f(t^n + \Delta t, u^n + \Delta t \cdot k_3)$

Stability constraint (CFL condition for diffusion):

Implicit Methods

Backward Euler:

Crank-Nicolson (second-order accurate):

BDF Methods (Backward Differentiation Formulas):

• BDF1: Backward Euler (1st order)
• BDF2: $\frac{3}{2}u^{n+1} - 2u^n + \frac{1}{2}u^{n-1} = \Delta t \cdot f^{n+1}$ (2nd order)

Characteristics:

• Unconditionally stable (A-stable)
• Requires nonlinear solver per time step
• Essential for stiff systems

Operator Splitting

Strang splitting for $\frac{\partial u}{\partial t} = Lu + Nu$ (linear + nonlinear):

Applications:

• Separate diffusion and reaction
• Different time scales for different physics
• Preserves second-order accuracy

3.3 Linear Algebra

Direct Methods

LU Factorization: $\mathbf{A} = \mathbf{L}\mathbf{U}$

Sparse direct solvers:

• PARDISO (Intel MKL)
• SuperLU
• MUMPS
• UMFPACK

Complexity: $O(N^\alpha)$ where $\alpha \approx 1.5-2$ for 3D problems

Iterative Methods

Conjugate Gradient (CG) for symmetric positive definite:

┌─────────────────────────────────────────────────────┐
│ r_0 = b - Ax_0                                      │
│ p_0 = r_0                                           │
│ for k = 0, 1, 2, ...                                │
│     α_k = (r_k^T r_k) / (p_k^T A p_k)               │
│     x_{k+1} = x_k + α_k p_k                         │
│     r_{k+1} = r_k - α_k A p_k                       │
│     β_k = (r_{k+1}^T r_{k+1}) / (r_k^T r_k)         │
│     p_{k+1} = r_{k+1} + β_k p_k                     │
└─────────────────────────────────────────────────────┘

GMRES (Generalized Minimal Residual) for non-symmetric systems

BiCGSTAB (Bi-Conjugate Gradient Stabilized)

Preconditioning

Purpose: Transform $\mathbf{A}\mathbf{x} = \mathbf{b}$ to $\mathbf{M}^{-1}\mathbf{A}\mathbf{x} = \mathbf{M}^{-1}\mathbf{b}$

Common preconditioners:

• ILU (Incomplete LU): Approximate factorization
• Multigrid: Hierarchical coarse-grid correction
• Domain decomposition: Parallel-friendly

Multigrid V-cycle:

3.4 Monte Carlo Methods

Particle-in-Cell (PIC) for Plasmas

Algorithm:

1. Push particles: $\mathbf{x}^{n+1} = \mathbf{x}^n + \mathbf{v}^n \Delta t$ 2. Weight to grid: $\rho_j = \sum_p q_p W(\mathbf{x}_p - \mathbf{x}_j)$ 3. Solve fields: $ abla^2 \phi = -\rho/\varepsilon_0$ 4. Interpolate to particles: $\mathbf{E}_p = \sum_j \mathbf{E}_j W(\mathbf{x}_p - \mathbf{x}_j)$ 5. Accelerate: $\mathbf{v}^{n+1} = \mathbf{v}^n + (q/m)\mathbf{E}_p \Delta t$

Monte Carlo Collisions: Null-collision method for efficiency

Direct Simulation Monte Carlo (DSMC)

For rarefied gas dynamics (high Knudsen number):

Algorithm:

1. Move particles (ballistic) 2. Index/sort particles into cells 3. Select collision pairs probabilistically 4. Perform collisions (conserve momentum, energy) 5. Sample macroscopic properties

Kinetic Monte Carlo (KMC)

For atomic-scale processes:

Rate calculation: $k_i = u_0 \exp\left(-\frac{E_a}{k_B T}\right)$

Event selection (BKL algorithm):

1. Calculate total rate: $R_{tot} = \sum_i k_i$ 2. Select event $j$ with probability $k_j / R_{tot}$ 3. Advance time: $\Delta t = -\ln(r) / R_{tot}$ where $r \in (0,1)$ 4. Execute event 5. Update rates

3.5 Interface Tracking

Level Set Methods

Interface = zero contour of $\phi(\mathbf{x}, t)$

Evolution equation:

Signed distance property: $| abla \phi| = 1$

Reinitialization (maintain distance property):

Advantages:

• Handles topological changes naturally
• Curvature: $\kappa =

abla \cdot \left( \frac{ abla \phi}{| abla \phi|} \right)$

• Normal: $\mathbf{n} = \frac{

abla \phi}{| abla \phi|}$

Fast Marching Method

For static Hamilton-Jacobi equations:

Complexity: $O(N \log N)$ using heap data structure

Application: Arrival time problems, distance computation

4. Key Application Areas

4.1 Lithography Simulation

Simulation Chain

┌─────────────────────────────────────────────────────┐
│ Mask (GDS) → Optical Simulation → Aerial Image →    │
│ → Resist Exposure → PEB Diffusion → Development →   │
│ → Final Profile                                     │
└─────────────────────────────────────────────────────┘

Hopkins Formulation (Partially Coherent Imaging)

Where:

• $J(f,g)$ = source intensity distribution
• $H(f,g)$ = pupil function
• $O(f,g)$ = mask spectrum

SOCS Decomposition

Sum of Coherent Systems:

• $\lambda_k$ = eigenvalues (decreasing)
• $h_k$ = eigenkernels
• Typically $N \sim 10-30$ sufficient

Rigorous Electromagnetic Methods

RCWA (Rigorous Coupled Wave Analysis):

• Fourier expansion of fields and permittivity
• Matrix eigenvalue problem per layer
• S-matrix or T-matrix propagation

FDTD (Finite Difference Time Domain):

• Yee grid staggering
• PML absorbing boundaries
• Handles arbitrary 3D structures

Resist Models

Dill exposure model:

Enhanced Fujita-Doolittle development:

4.2 Plasma Process Modeling

Multi-Scale Framework

┌─────────────────────────────────────────────────────┐
│ Reactor Scale (cm)     Feature Scale (nm)           │
│       ↓                      ↑                      │
│   Plasma Model    →    Flux/Distributions           │
│       ↓                      ↑                      │
│   Surface Fluxes   →   Profile Evolution            │
└─────────────────────────────────────────────────────┘

Fluid Plasma Model

Continuity:

Momentum (drift-diffusion):

Energy:

Poisson:

Feature-Scale Model

Surface advancement:

Where:

• $\Gamma_{ion}$ = ion flux
• $Y_{ion}$ = ion-enhanced yield (angle, energy dependent)
• $S_{chem}$ = chemical sticking coefficient
• $\Gamma_{dep}$ = deposition flux

4.3 TCAD Device Simulation

Scharfetter-Gummel Discretization

Current between nodes $i$ and $j$:

Bernoulli function:

Properties:

• Exact for constant field
• Numerically stable for large bias
• Preserves current continuity

Quantum Corrections

Density gradient model:

Schrödinger-Poisson (1D slice):

5. Multi-Scale and Multi-Physics Coupling

5.1 Length Scale Hierarchy

┌─────────────────────────────────────────────────────┐
│ Atomic     Feature    Device      Die       Wafer   │
│ (0.1 nm)   (10 nm)   (100 nm)   (1 mm)   (300 mm)   │
│    │          │          │         │         │      │
│    └────┬─────┴────┬─────┴────┬────┴────┬────┘      │
│         │          │          │         │           │
│       Ab initio   KMC    Continuum   Pattern        │
│        DFT       MD        PDE       Effects        │
└─────────────────────────────────────────────────────┘

5.2 Coupling Approaches

Sequential (Parameter Passing)

┌─────────────────────────────────────────────────────┐
│ Lower Scale → Parameters → Higher Scale             │
└─────────────────────────────────────────────────────┘

Examples:

• DFT → activation energies → KMC rates
• MD → surface diffusion coefficients → continuum
• Feature-scale → pattern density → wafer-scale

Concurrent (Domain Decomposition)

Different physics in different regions, coupled at interfaces:

Handshaking region:

Force matching or energy-based coupling

Homogenization

Effective properties from microstructure:

Application: Pattern-density effects in CMP

5.3 Multi-Physics Coupling

Monolithic vs. Partitioned

Monolithic: Solve all physics simultaneously

• Strong coupling
• Large, often ill-conditioned systems

Partitioned: Iterate between physics

while not converged:
    Solve Physics 1 with fixed Physics 2 variables
    Solve Physics 2 with fixed Physics 1 variables
    Check convergence

• Reuse existing solvers
• May have stability issues

6. Uncertainty Quantification

6.1 Sources of Uncertainty

• Process variations: Dose, focus, temperature, pressure
• Material variations: Film thickness, composition, defect density
• Model uncertainty: Parameter calibration, structural assumptions
• Measurement noise: Metrology errors

6.2 Polynomial Chaos Expansion

Expansion:

Where:

• $\boldsymbol{\xi}$ = random variables (inputs)
• $\Psi_k$ = orthogonal polynomial basis
• $u_k$ = deterministic coefficients

Basis selection:

Statistics from coefficients:

• Mean: $\mathbb{E}[u] = u_0$
• Variance: $\text{Var}[u] = \sum_{k=1}^{P} u_k^2 \langle \Psi_k^2 \rangle$

6.3 Stochastic Collocation

Algorithm:

1. Select collocation points $\boldsymbol{\xi}^{(q)}$ (Gauss quadrature, sparse grids) 2. Solve deterministic problem at each point 3. Construct interpolant/response surface 4. Compute statistics by integration

Advantages:

• Non-intrusive (uses existing solvers)
• Flexible basis
• Good for smooth dependence on parameters

6.4 Sensitivity Analysis

Sobol indices (variance decomposition):

First-order index:

Total effect index:

Interpretation:

• High $S_i$ → Parameter $i$ is important
• $S_{T_i} - S_i$ → Interaction effects

7. Emerging Methods

7.1 Physics-Informed Neural Networks (PINNs)

Loss function:

Physics loss (PDE residual):

Where $\mathcal{N}$ is the differential operator.

Architecture:

Input (x,t) → Hidden Layers → Output u
                 ↓
         Automatic Differentiation
                 ↓
         Physics Residual Loss

Applications:

• Inverse problems
• Surrogate models
• Data assimilation

7.2 Machine Learning Surrogates

Gaussian Process Regression

Prior: $u(\mathbf{x}) \sim \mathcal{GP}(m(\mathbf{x}), k(\mathbf{x}, \mathbf{x}'))$

Common kernels:

• Squared exponential: $k(\mathbf{x}, \mathbf{x}') = \sigma^2 \exp\left(-\frac{|\mathbf{x} - \mathbf{x}'|^2}{2\ell^2}\right)$
• Matérn

Prediction with uncertainty:

Deep Neural Networks

Convolutional Neural Networks (CNNs):

• Image-to-image mapping
• Lithography: mask → resist profile
• Defect detection

Graph Neural Networks (GNNs):

• Irregular mesh data
• Molecular property prediction

7.3 Reduced Order Models

Proper Orthogonal Decomposition (POD)

Snapshot matrix: $\mathbf{X} = [\mathbf{u}_1, \mathbf{u}_2, \ldots, \mathbf{u}_M]$

SVD: $\mathbf{X} = \mathbf{U} \boldsymbol{\Sigma} \mathbf{V}^T$

Reduced basis: First $r$ columns of $\mathbf{U}$

Projection: $\mathbf{u} \approx \mathbf{U}_r \mathbf{a}$ where $\mathbf{a} \in \mathbb{R}^r$

Reduced system:

8. Computational Considerations

8.1 Challenges

• Memory: 3D simulations with fine grids → TB scale
• Time: Full process simulation → days to weeks
• Coupling: Different physics require different solvers
• Scales: 10+ orders of magnitude in length
• Stiffness: Vastly different time scales

8.2 High-Performance Computing Solutions

Parallel Algorithms

Domain decomposition:

Communication patterns:

• Ghost cell exchange (explicit methods)
• Iterative subdomain coupling (implicit)

Scalability:

• Strong scaling: Fixed problem, more processors
• Weak scaling: Problem grows with processors

GPU Acceleration

Suitable algorithms:

• Matrix-vector products (SpMV)
• FFT-based convolutions
• Monte Carlo sampling
• Stencil operations

Programming models: CUDA, OpenCL, HIP, SYCL

Adaptive Methods

Mesh refinement (AMR):

• Refine where gradients are large
• Coarsen where solution is smooth
• Quadtree/octree data structures

Time adaptivity:

8.3 Software Ecosystem

Summary Table

Source: ChipFoundryServices — Search this topic — Ask CFSGPT