Mathematical Modeling of Epitaxy in Semiconductor Front-End Processing (FEP)
1. Overview
Epitaxy is a critical Front-End Process (FEP) step where crystalline films are grown on crystalline substrates with precise control of:
- Thickness
- Composition
- Doping concentration
- Defect density
Mathematical modeling enables:
- Process optimization
- Defect prediction
- Virtual fabrication
- Equipment design
1.1 Types of Epitaxy
- Homoepitaxy: Same material as substrate (e.g., Si on Si)
- Heteroepitaxy: Different material from substrate (e.g., GaAs on Si, SiGe on Si)
1.2 Epitaxy Methods
- Vapor Phase Epitaxy (VPE) / Chemical Vapor Deposition (CVD)
- Atmospheric Pressure CVD (APCVD)
- Low Pressure CVD (LPCVD)
- Metal-Organic CVD (MOCVD)
- Molecular Beam Epitaxy (MBE)
- Liquid Phase Epitaxy (LPE)
- Solid Phase Epitaxy (SPE)
2. Fundamental Thermodynamic Framework
2.1 Driving Force for Growth
The supersaturation provides the thermodynamic driving force:
$$
\Delta \mu = k_B T \ln\left(\frac{P}{P_{eq}}\right)
$$
Where:
- $\Delta \mu$ = chemical potential difference (driving force)
- $k_B$ = Boltzmann's constant ($1.38 \times 10^{-23}$ J/K)
- $T$ = absolute temperature (K)
- $P$ = actual partial pressure of precursor
- $P_{eq}$ = equilibrium vapor pressure
2.2 Free Energy of Mixing (Multi-component Systems)
For systems like SiGe alloys:
$$
\Delta G_{mix} = RT\left(x \ln x + (1-x) \ln(1-x)\right) + \Omega x(1-x)
$$
Where:
- $R$ = universal gas constant (8.314 J/mol·K)
- $x$ = mole fraction of component
- $\Omega$ = interaction parameter (regular solution model)
2.3 Gibbs Free Energy of Formation
$$
\Delta G = \Delta H - T\Delta S
$$
For spontaneous growth: $\Delta G < 0$
3. Growth Rate Kinetics
3.1 The Two-Regime Model
Epitaxial growth rate is governed by two competing mechanisms:
Overall growth rate equation:
$$
G = \frac{k_s \cdot h_g \cdot C_g}{k_s + h_g}
$$
Where:
- $G$ = growth rate (nm/min or μm/min)
- $k_s$ = surface reaction rate constant
- $h_g$ = gas-phase mass transfer coefficient
- $C_g$ = gas-phase reactant concentration
3.2 Temperature Dependence
The surface reaction rate follows Arrhenius behavior:
$$
k_s = A \exp\left(-\frac{E_a}{k_B T}\right)
$$
Where:
- $A$ = pre-exponential factor (frequency factor)
- $E_a$ = activation energy (eV or J/mol)
3.3 Growth Rate Regimes
| Temperature Regime | Limiting Factor | Growth Rate Expression | Temperature Dependence |
|:-------------------|:----------------|:-----------------------|:-----------------------|
| Low T | Surface reaction | $G \approx k_s \cdot C_g$ | Strong (exponential) |
| High T | Mass transport | $G \approx h_g \cdot C_g$ | Weak (~$T^{1.5-2}$) |
3.4 Boundary Layer Analysis
For horizontal CVD reactors, the boundary layer thickness evolves as:
$$
\delta(x) = \sqrt{\frac{
u \cdot x}{v_{\infty}}}
$$
Where:
- $\delta(x)$ = boundary layer thickness at position $x$
- $
u$ = kinematic viscosity (m²/s)
- $x$ = distance from gas inlet (m)
- $v_{\infty}$ = free stream gas velocity (m/s)
The mass transfer coefficient:
$$
h_g = \frac{D_{gas}}{\delta}
$$
Where $D_{gas}$ is the gas-phase diffusion coefficient.
4. Surface Kinetics: BCF Theory
The Burton-Cabrera-Frank (BCF) model describes atomic-scale growth mechanisms.
4.1 Surface Diffusion Equation
$$
D_s
abla^2 n_s - \frac{n_s - n_{eq}}{\tau_s} + J_{ads} = 0
$$
Where:
- $n_s$ = adatom surface density (atoms/cm²)
- $D_s$ = surface diffusion coefficient (cm²/s)
- $n_{eq}$ = equilibrium adatom density
- $\tau_s$ = mean adatom lifetime before desorption (s)
- $J_{ads}$ = adsorption flux (atoms/cm²·s)
4.2 Characteristic Diffusion Length
$$
\lambda_s = \sqrt{D_s \tau_s}
$$
This parameter determines the growth mode:
- Step-flow growth: $\lambda_s > L$ (terrace width)
- 2D nucleation growth: $\lambda_s < L$
4.3 Surface Diffusion Coefficient
$$
D_s = D_0 \exp\left(-\frac{E_m}{k_B T}\right)
$$
Where:
- $D_0$ = pre-exponential factor (~$10^{-3}$ cm²/s)
- $E_m$ = migration energy barrier (eV)
4.4 Step Velocity
$$
v_{step} = \frac{2 D_s (n_s - n_{eq})}{\lambda_s} \tanh\left(\frac{L}{2\lambda_s}\right)
$$
Where $L$ is the inter-step spacing (terrace width).
4.5 Growth Rate from Step Flow
$$
G = \frac{v_{step} \cdot h_{step}}{L}
$$
Where $h_{step}$ is the step height (monolayer thickness).
5. Heteroepitaxy and Strain Modeling
5.1 Lattice Mismatch
$$
f = \frac{a_{film} - a_{substrate}}{a_{substrate}}
$$
Where:
- $f$ = lattice mismatch (dimensionless, often expressed as %)
- $a_{film}$ = lattice constant of film material
- $a_{substrate}$ = lattice constant of substrate
Example values:
| System | Lattice Mismatch |
|:-------|:-----------------|
| Si₀.₇Ge₀.₃ on Si | ~1.2% |
| Ge on Si | ~4.2% |
| GaAs on Si | ~4.0% |
| InAs on GaAs | ~7.2% |
| GaN on Sapphire | ~16% |
5.2 Strain Components
For biaxial strain in (001) films:
$$
\varepsilon_{xx} = \varepsilon_{yy} = \varepsilon_{\parallel} = \frac{a_s - a_f}{a_f} \approx -f
$$
$$
\varepsilon_{zz} = \varepsilon_{\perp} = -\frac{2C_{12}}{C_{11}} \varepsilon_{\parallel}
$$
Where $C_{11}$ and $C_{12}$ are elastic constants.
5.3 Elastic Energy
For a coherently strained film:
$$
E_{elastic} = \frac{2G(1+
u)}{1-
u} f^2 h = M f^2 h
$$
Where:
- $G$ = shear modulus (Pa)
- $
u$ = Poisson's ratio
- $h$ = film thickness
- $M$ = biaxial modulus = $\frac{2G(1+
u)}{1-
u}$
5.4 Critical Thickness (Matthews-Blakeslee)
$$
h_c = \frac{b}{8\pi f(1+
u)} \left[\ln\left(\frac{h_c}{b}\right) + 1\right]
$$
Where:
- $h_c$ = critical thickness for dislocation formation
- $b$ = Burgers vector magnitude
- $f$ = lattice mismatch
- $
u$ = Poisson's ratio
5.5 People-Bean Approximation (for SiGe)
Empirical formula:
$$
h_c \approx \frac{0.55}{f^2} \text{ (nm, with } f \text{ as a decimal)}
$$
Or equivalently:
$$
h_c \approx \frac{5500}{x^2} \text{ (nm, for Si}_{1-x}\text{Ge}_x\text{)}
$$
5.6 Threading Dislocation Density
Above critical thickness, dislocation density evolves:
$$
\rho_{TD}(h) = \rho_0 \exp\left(-\frac{h}{h_0}\right) + \rho_{\infty}
$$
Where:
- $\rho_{TD}$ = threading dislocation density (cm⁻²)
- $\rho_0$ = initial density
- $h_0$ = characteristic decay length
- $\rho_{\infty}$ = residual density
6. Reactor-Scale Modeling
6.1 Coupled Transport Equations
6.1.1 Momentum Conservation (Navier-Stokes)
$$
\rho\left(\frac{\partial \mathbf{v}}{\partial t} + \mathbf{v} \cdot
abla \mathbf{v}\right) = -
abla p + \mu
abla^2 \mathbf{v} + \rho \mathbf{g}
$$
Where:
- $\rho$ = gas density (kg/m³)
- $\mathbf{v}$ = velocity vector (m/s)
- $p$ = pressure (Pa)
- $\mu$ = dynamic viscosity (Pa·s)
- $\mathbf{g}$ = gravitational acceleration
6.1.2 Continuity Equation
$$
\frac{\partial \rho}{\partial t} +
abla \cdot (\rho \mathbf{v}) = 0
$$
6.1.3 Species Transport
$$
\frac{\partial C_i}{\partial t} + \mathbf{v} \cdot
abla C_i = D_i
abla^2 C_i + R_i
$$
Where:
- $C_i$ = concentration of species $i$ (mol/m³)
- $D_i$ = diffusion coefficient of species $i$ (m²/s)
- $R_i$ = net reaction rate (mol/m³·s)
6.1.4 Energy Conservation
$$
\rho c_p \left(\frac{\partial T}{\partial t} + \mathbf{v} \cdot
abla T\right) = k
abla^2 T + \sum_j \Delta H_j r_j
$$
Where:
- $c_p$ = specific heat capacity (J/kg·K)
- $k$ = thermal conductivity (W/m·K)
- $\Delta H_j$ = enthalpy of reaction $j$ (J/mol)
- $r_j$ = rate of reaction $j$ (mol/m³·s)
6.2 Silicon CVD Chemistry
6.2.1 From Silane (SiH₄)
Gas phase decomposition:
$$
\text{SiH}_4 \xrightarrow{k_1} \text{SiH}_2 + \text{H}_2
$$
Surface reaction:
$$
\text{SiH}_2(g) + * \xrightarrow{k_2} \text{Si}(s) + \text{H}_2(g)
$$
Where $*$ denotes a surface site.
6.2.2 From Dichlorosilane (DCS)
$$
\text{SiH}_2\text{Cl}_2 \rightarrow \text{SiCl}_2 + \text{H}_2
$$
$$
\text{SiCl}_2 + \text{H}_2 \rightarrow \text{Si}(s) + 2\text{HCl}
$$
6.2.3 Rate Law
$$
r_{dep} = k_2 P_{SiH_2} (1 - \theta)
$$
Where:
- $P_{SiH_2}$ = partial pressure of SiH₂
- $\theta$ = surface site coverage
6.3 Dimensionless Numbers
| Number | Definition | Physical Meaning |
|:-------|:-----------|:-----------------|
| Reynolds | $Re = \frac{\rho v L}{\mu}$ | Inertia vs. viscous forces |
| Prandtl | $Pr = \frac{\mu c_p}{k}$ | Momentum vs. thermal diffusivity |
| Schmidt | $Sc = \frac{\mu}{\rho D}$ | Momentum vs. mass diffusivity |
| Damköhler | $Da = \frac{k_s L}{D}$ | Reaction rate vs. diffusion rate |
| Grashof | $Gr = \frac{g \beta \Delta T L^3}{
u^2}$ | Buoyancy vs. viscous forces |
7. Selective Epitaxial Growth (SEG) Modeling
7.1 Overview
In SEG, growth occurs on exposed Si but not on dielectric (SiO₂/Si₃N₄).
7.2 Loading Effect Model
$$
G_{local} = G_0 \left(1 + \alpha \cdot \frac{A_{mask}}{A_{Si}}\right)
$$
Where:
- $G_{local}$ = local growth rate
- $G_0$ = baseline growth rate
- $\alpha$ = pattern sensitivity factor
- $A_{mask}$ = dielectric (mask) area
- $A_{Si}$ = exposed silicon area
7.3 Pattern-Dependent Growth
Sources of non-uniformity:
- Local depletion of reactants over Si regions
- Species reflected/desorbed from mask contribute to nearby Si
- Gas-phase diffusion length effects
7.4 Selectivity Condition
For selective growth on Si vs. oxide:
$$
r_{deposition,Si} > 0 \quad \text{and} \quad r_{deposition,oxide} < r_{etching,oxide}
$$
Achieved by adding HCl:
$$
\text{Si}(nuclei) + 2\text{HCl} \rightarrow \text{SiCl}_2 + \text{H}_2
$$
Nuclei on oxide are etched before they can grow, maintaining selectivity.
7.5 Faceting Model
Growth rate depends on crystallographic orientation:
$$
G_{(hkl)} = G_0 \cdot f(hkl) \cdot \exp\left(-\frac{E_{a,(hkl)}}{k_B T}\right)
$$
Typical growth rate hierarchy:
$$
G_{(100)} > G_{(110)} > G_{(111)}
$$
8. Dopant Incorporation
8.1 Segregation Coefficient
Equilibrium segregation coefficient:
$$
k_0 = \frac{C_{solid}}{C_{liquid/gas}}
$$
Effective segregation coefficient:
$$
k_{eff} = \frac{k_0}{k_0 + (1-k_0)\exp\left(-\frac{G\delta}{D_l}\right)}
$$
Where:
- $k_0$ = equilibrium segregation coefficient
- $G$ = growth rate
- $\delta$ = boundary layer thickness
- $D_l$ = diffusivity in liquid/gas phase
8.2 Dopant Concentration in Film
$$
C_{film} = k_{eff} \cdot C_{gas}
$$
8.3 Dopant Profile Abruptness
The transition width is limited by:
- Surface segregation length: $\lambda_{seg}$
- Diffusion during growth: $L_D = \sqrt{D \cdot t}$
- Autodoping from substrate
$$
\Delta z_{transition} \approx \sqrt{\lambda_{seg}^2 + L_D^2}
$$
8.4 Common Dopants for Si Epitaxy
| Dopant | Type | Precursor | Segregation Behavior |
|:-------|:-----|:----------|:---------------------|
| B | p-type | B₂H₆, BCl₃ | Low segregation |
| P | n-type | PH₃, PCl₃ | Moderate segregation |
| As | n-type | AsH₃ | Strong segregation |
| Sb | n-type | SbH₃ | Very strong segregation |
9. Atomistic Simulation Methods
9.1 Kinetic Monte Carlo (KMC)
9.1.1 Event Rates
Each atomic event has a rate following Arrhenius:
$$
\Gamma_i =
u_0 \exp\left(-\frac{E_i}{k_B T}\right)
$$
Where:
- $\Gamma_i$ = rate of event $i$ (s⁻¹)
- $
u_0$ = attempt frequency (~10¹²-10¹³ s⁻¹)
- $E_i$ = activation energy for event $i$
9.1.2 Events Modeled
- Adsorption: $\Gamma_{ads} = \frac{P}{\sqrt{2\pi m k_B T}} \cdot s$
- Desorption: $\Gamma_{des} =
u_0 \exp(-E_{des}/k_B T)$
- Surface diffusion: $\Gamma_{diff} =
u_0 \exp(-E_m/k_B T)$
- Step attachment: $\Gamma_{attach}$
- Step detachment: $\Gamma_{detach}$
9.1.3 Time Advancement
$$
\Delta t = -\frac{\ln(r)}{\Gamma_{total}} = -\frac{\ln(r)}{\sum_i \Gamma_i}
$$
Where $r$ is a uniform random number in $(0,1]$.
9.2 Density Functional Theory (DFT)
Provides input parameters for KMC:
- Adsorption energies
- Migration barriers
- Surface reconstruction energetics
- Reaction pathways
Kohn-Sham equation:
$$
\left[-\frac{\hbar^2}{2m}
abla^2 + V_{eff}(\mathbf{r})\right]\psi_i(\mathbf{r}) = \varepsilon_i \psi_i(\mathbf{r})
$$
9.3 Molecular Dynamics (MD)
Newton's equations:
$$
m_i \frac{d^2 \mathbf{r}_i}{dt^2} = -
abla_i U(\mathbf{r}_1, \mathbf{r}_2, ..., \mathbf{r}_N)
$$
Where $U$ is the interatomic potential (e.g., Stillinger-Weber, Tersoff for Si).
10. Nucleation Theory
10.1 Classical Nucleation Theory (CNT)
10.1.1 Gibbs Free Energy Change
$$
\Delta G(r) = -\frac{4}{3}\pi r^3 \cdot \frac{\Delta \mu}{\Omega} + 4\pi r^2 \gamma
$$
Where:
- $r$ = nucleus radius
- $\Delta \mu$ = supersaturation (driving force)
- $\Omega$ = atomic volume
- $\gamma$ = surface energy
10.1.2 Critical Nucleus Radius
Setting $\frac{d(\Delta G)}{dr} = 0$:
$$
r^* = \frac{2\gamma \Omega}{\Delta \mu}
$$
10.1.3 Free Energy Barrier
$$
\Delta G^* = \frac{16 \pi \gamma^3 \Omega^2}{3 (\Delta \mu)^2}
$$
10.1.4 Nucleation Rate
$$
J = Z \beta^ N_s \exp\left(-\frac{\Delta G^}{k_B T}\right)
$$
Where:
- $J$ = nucleation rate (nuclei/cm²·s)
- $Z$ = Zeldovich factor (~0.01-0.1)
- $\beta^*$ = attachment rate to critical nucleus
- $N_s$ = surface site density
10.2 Growth Modes
| Mode | Surface Energy Condition | Growth Behavior | Example |
|:-----|:-------------------------|:----------------|:--------|
| Frank-van der Merwe | $\gamma_s \geq \gamma_f + \gamma_{int}$ | Layer-by-layer (2D) | Si on Si |
| Volmer-Weber | $\gamma_s < \gamma_f + \gamma_{int}$ | Island (3D) | Metals on oxides |
| Stranski-Krastanov | Intermediate | 2D then 3D islands | InAs/GaAs QDs |
10.3 2D Nucleation
Critical island size (atoms):
$$
i^* = \frac{\pi \gamma_{step}^2 \Omega}{(\Delta \mu)^2 k_B T}
$$
11. TCAD Process Simulation
11.1 Overview
Tools: Synopsys Sentaurus Process, Silvaco Victory Process
11.2 Diffusion-Reaction System
$$
\frac{\partial C_i}{\partial t} =
abla \cdot (D_i
abla C_i - \mu_i C_i
abla \phi) + G_i - R_i
$$
Where:
- First term: Fickian diffusion
- Second term: Drift in electric field (for charged species)
- $G_i$ = generation rate
- $R_i$ = recombination rate
11.3 Point Defect Dynamics
Vacancy concentration:
$$
\frac{\partial C_V}{\partial t} = D_V
abla^2 C_V + G_V - k_{IV} C_I C_V
$$
Interstitial concentration:
$$
\frac{\partial C_I}{\partial t} = D_I
abla^2 C_I + G_I - k_{IV} C_I C_V
$$
Where $k_{IV}$ is the recombination rate constant.
11.4 Stress Evolution
Equilibrium equation:
$$
abla \cdot \boldsymbol{\sigma} = 0
$$
Constitutive relation:
$$
\boldsymbol{\sigma} = \mathbf{C} : (\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^{thermal} - \boldsymbol{\varepsilon}^{intrinsic})
$$
Where:
- $\boldsymbol{\sigma}$ = stress tensor
- $\mathbf{C}$ = elastic stiffness tensor
- $\boldsymbol{\varepsilon}$ = total strain
- $\boldsymbol{\varepsilon}^{thermal}$ = thermal strain = $\alpha \Delta T$
- $\boldsymbol{\varepsilon}^{intrinsic}$ = intrinsic strain (lattice mismatch)
11.5 Level Set Method for Interface Tracking
$$
\frac{\partial \phi}{\partial t} + v_n |
abla \phi| = 0
$$
Where:
- $\phi$ = level set function (interface at $\phi = 0$)
- $v_n$ = interface normal velocity
12. Advanced Topics
12.1 Atomic Layer Epitaxy (ALE) / Atomic Layer Deposition (ALD)
Self-limiting surface reactions modeled as Langmuir kinetics:
$$
\theta = \frac{K \cdot P \cdot t}{1 + K \cdot P \cdot t} \rightarrow 1 \quad \text{as } t \rightarrow \infty
$$
Growth per cycle (GPC):
$$
GPC = \theta_{sat} \cdot d_{monolayer}
$$
Typical GPC values: 0.5-1.5 Å/cycle
12.2 III-V on Silicon Integration
Challenges and models:
- Anti-phase boundaries (APBs): Form at single-step terraces
- Threading dislocations: $\rho_{TD} \propto f^2$ initially
- Thermal mismatch stress: $\sigma_{thermal} = \frac{E \Delta \alpha \Delta T}{1-
u}$
12.3 Quantum Dot Formation (Stranski-Krastanov)
Critical thickness for islanding:
$$
h_{SK} \approx \frac{\gamma}{M f^2}
$$
Island density:
$$
n_{island} \propto \exp\left(-\frac{E_{island}}{k_B T}\right) \cdot F^{1/3}
$$
Where $F$ is the deposition flux.
12.4 Machine Learning in Epitaxy Modeling
Physics-Informed Neural Networks (PINNs):
$$
\mathcal{L}_{total} = \mathcal{L}_{data} + \lambda_{PDE}\mathcal{L}_{physics} + \lambda_{BC}\mathcal{L}_{boundary}
$$
Where:
- $\mathcal{L}_{data}$ = data fitting loss
- $\mathcal{L}_{physics}$ = PDE residual loss
- $\mathcal{L}_{boundary}$ = boundary condition loss
- $\lambda$ = weighting parameters
Applications:
- Surrogate models for reactor optimization
- Inverse problems (parameter extraction)
- Process window optimization
- Defect prediction
13. Key Equations
| Phenomenon | Key Equation | Primary Parameters |
|:-----------|:-------------|:-------------------|
| Growth rate (dual regime) | $G = \frac{k_s h_g C_g}{k_s + h_g}$ | Temperature, pressure, flow |
| Surface diffusion length | $\lambda_s = \sqrt{D_s \tau_s}$ | Temperature |
| Lattice mismatch | $f = \frac{a_f - a_s}{a_s}$ | Material system |
| Critical thickness | $h_c = \frac{b}{8\pi f(1+
u)}\left[\ln\frac{h_c}{b}+1\right]$ | Mismatch, Burgers vector |
| Elastic strain energy | $E = M f^2 h$ | Mismatch, thickness, modulus |
| Nucleation rate | $J \propto \exp(-\Delta G^*/k_BT)$ | Supersaturation, surface energy |
| Species transport | $\frac{\partial C}{\partial t} + \mathbf{v}\cdot
abla C = D
abla^2 C + R$ | Diffusivity, velocity, reactions |
| KMC event rate | $\Gamma =
u_0 \exp(-E_a/k_BT)$ | Activation energy, temperature |
Physical Constants
| Constant | Symbol | Value |
|:---------|:-------|:------|
| Boltzmann constant | $k_B$ | $1.38 \times 10^{-23}$ J/K |
| Gas constant | $R$ | 8.314 J/mol·K |
| Planck constant | $h$ | $6.63 \times 10^{-34}$ J·s |
| Electron charge | $e$ | $1.60 \times 10^{-19}$ C |
| Si lattice constant | $a_{Si}$ | 5.431 Å |
| Ge lattice constant | $a_{Ge}$ | 5.658 Å |
| GaAs lattice constant | $a_{GaAs}$ | 5.653 Å |