Mathematical Modeling of Diffusion and Ion Implantation in Semiconductor Manufacturing

Keywords: diffusion and ion implantation,diffusion,ion implantation,dopant diffusion,fick law,implant profile,gaussian profile,pearson distribution,ted,transient enhanced diffusion,thermal budget,semiconductor doping

Mathematical Modeling of Diffusion and Ion Implantation in Semiconductor Manufacturing

Part I: Diffusion Modeling

Fundamental Equations

Dopant redistribution in silicon at elevated temperatures is governed by Fick's Laws .

Fick's First Law

Relates flux to concentration gradient:

Where:

• $J$ — Atomic flux (atoms/cm²·s)
• $D$ — Diffusion coefficient (cm²/s)
• $C$ — Concentration (atoms/cm³)
• $x$ — Position (cm)

Fick's Second Law

The diffusion equation follows from continuity:

This parabolic PDE admits analytical solutions for idealized boundary conditions.

Temperature Dependence

The diffusion coefficient follows an Arrhenius relationship :

Parameters:

• $D_0$ — Pre-exponential factor (cm²/s)
• $E_a$ — Activation energy (eV)
• $k$ — Boltzmann's constant ($8.617 \times 10^{-5}$ eV/K)
• $T$ — Absolute temperature (K)

Typical Values for Phosphorus in Silicon:

Diffusion approximately doubles every 10–15°C near typical process temperatures (900–1100°C).

Classical Analytical Solutions

Case 1: Constant Surface Concentration (Predeposition)

Boundary Conditions:

• $C(0, t) = C_s$ (constant surface concentration)
• $C(\infty, t) = 0$ (zero at infinite depth)
• $C(x, 0) = 0$ (initially undoped)

Solution:

Complementary Error Function:

Total Incorporated Dose:

Case 2: Fixed Dose (Drive-in Diffusion)

Boundary Conditions:

• $\displaystyle\int_0^{\infty} C \, dx = Q$ (constant total dose)
• $\displaystyle\frac{\partial C}{\partial x}\bigg|_{x=0} = 0$ (no flux at surface)

Solution (Gaussian Profile):

Peak Surface Concentration:

Junction Depth Calculation

The metallurgical junction forms where dopant concentration equals background doping $C_B$.

For erfc Profile:

For Gaussian Profile:

Concentration-Dependent Diffusion

At high doping concentrations (approaching or exceeding intrinsic carrier concentration $n_i$), diffusivity becomes concentration-dependent.

Generalized Model:

Physical Interpretation:

Resulting Nonlinear PDE:

This requires numerical solution methods.

Point Defect Mediated Diffusion

Modern process modeling couples dopant diffusion to point defect dynamics.

Governing System of PDEs:

Variable Definitions:

• $C_I$ — Interstitial concentration
• $C_V$ — Vacancy concentration
• $C_A$ — Dopant atom concentration
• $k_{IV}$ — Interstitial-vacancy recombination rate
• $G$ — Generation rate
• $R$ — Surface recombination rate

Part II: Ion Implantation Modeling

Energy Loss Mechanisms

Implanted ions lose energy through two mechanisms:

Total Stopping Power:

Nuclear Stopping (Elastic Collisions)

Dominates at low energies :

Where:

• $\gamma = \displaystyle\frac{4 M_1 M_2}{(M_1 + M_2)^2}$ — Energy transfer factor
• $a$ — Screening length
• $s_n(\varepsilon)$ — Reduced nuclear stopping

Electronic Stopping (Inelastic Interactions)

Dominates at high energies :

(at intermediate energies)

LSS Theory

Lindhard, Scharff, and Schiøtt developed universal scaling using reduced units.

Reduced Energy:

Reduced Path Length:

This allows tabulation of universal range curves applicable across ion-target combinations.

Gaussian Profile Approximation

First-Order Implant Profile:

Parameters:

Peak Concentration:

Higher-Order Moment Distributions

The Gaussian approximation fails for many practical cases. The Pearson IV distribution uses four statistical moments:

Pearson IV Form:

Parameters $(a, b, m, u, K)$ are derived from the four moments through algebraic relations.

Skewness Behavior:

• Light ions (B) in heavy substrates → Negative skewness (tail toward surface)
• Heavy ions (As, Sb) in silicon → Positive skewness (tail toward bulk)

Dual Pearson Model

For channeling tails or complex profiles:

Where:

• $C_1(x)$, $C_2(x)$ — Two Pearson distributions with different parameters
• $f$ — Weight fraction

Lateral Distribution

Ions scatter laterally as well:

For Amorphous Targets:

Lateral straggle is critical for device scaling—it limits minimum feature sizes.

Monte Carlo Simulation (TRIM/SRIM)

For accurate profiles, especially in multilayer or crystalline structures, Monte Carlo methods track individual ion trajectories.

Algorithm:

1. Initialize ion position, direction, energy 2. Select free flight path: $\lambda = 1/(N\pi a^2)$ 3. Calculate impact parameter and scattering angle via screened Coulomb potential 4. Energy transfer to recoil: $$T = T_m \sin^2\left(\frac{\theta}{2}\right)$$ where $T_m = \gamma E$ 5. Apply electronic energy loss over path segment 6. Update ion position/direction; cascade recoils if $T > E_d$ (displacement energy) 7. Repeat until $E < E_{\text{cutoff}}$ 8. Accumulate statistics over $10^4 - 10^6$ ion histories

ZBL Interatomic Potential:

Where $\phi$ is the screening function tabulated from quantum mechanical calculations.

Channeling

In crystalline silicon, ions aligned with crystal axes experience reduced stopping.

Critical Angle for Channeling:

Where:

• $d$ — Atomic spacing along the channel
• $E$ — Ion energy

Effects:

• Channeled ions penetrate 2–10× deeper
• Creates extended tails in profiles
• Modern implants use 7° tilt or random-equivalent conditions to minimize

Damage Accumulation

Implant damage is quantified by:

Where:

• $

u(E)$ — Kinchin-Pease damage function (displaced atoms per ion)

• $F(x, E)$ — Energy deposition profile

Amorphization Threshold for Silicon:

(approximately 10–15% of atoms displaced)

Part III: Post-Implant Diffusion and Transient Enhanced Diffusion

Transient Enhanced Diffusion (TED)

After implantation, excess interstitials dramatically enhance diffusion until they anneal:

Where:

• $C_I^*$ — Equilibrium interstitial concentration

"+1" Model for Boron:

Impact: TED can cause junction depths 2–5× deeper than equilibrium diffusion would predict—critical for modern shallow junctions.

{311} Defect Dissolution Kinetics

Interstitials cluster into rod-like {311} defects that slowly dissolve:

The released interstitials sustain TED, explaining why TED persists for times much longer than point defect diffusion would suggest.

Part IV: Numerical Methods

Finite Difference Discretization

For the diffusion equation on uniform grid $(x_i, t_n)$:

Explicit (Forward Euler)

Stability Requirement (CFL Condition):

Implicit (Backward Euler)

• Unconditionally stable
• Requires solving tridiagonal system each timestep

Crank-Nicolson Method

• Average of explicit and implicit schemes
• Second-order accurate in time
• Results in tridiagonal system

Adaptive Meshing

Concentration gradients vary by orders of magnitude. Adaptive grids refine near:

• Junctions
• Surface
• Implant peaks
• Moving interfaces

Grid Spacing Scaling:

Process Simulation Flow (TCAD)

Modern simulators (Sentaurus Process, ATHENA, FLOOPS) integrate:

1. Implantation → Monte Carlo or analytical tables 2. Damage model → Amorphization, defect clustering 3. Annealing → Coupled dopant-defect PDEs 4. Oxidation → Deal-Grove kinetics, stress effects, OED 5. Silicidation, epitaxy, etc. → Specialized models

Output feeds device simulation (drift-diffusion, Monte Carlo transport).

Part V: Key Process Design Equations

Thermal Budget

The characteristic diffusion length after multiple thermal steps:

For Varying Temperature $T(t)$:

Sheet Resistance

For Uniform Mobility Approximation:

Electrical measurements to profile parameters.

Implant Dose-Energy Selection

Target Peak Concentration:

Target Depth (Empirical):

Where:

• $n \approx 0.6 - 0.8$ (depending on energy regime)
• $A$ — Ion-target dependent constant

Key Mathematical Tools:

abla \cdot (D abla C)$ | Analytical (erfc, Gaussian) or FEM/FDM |

abla \cdot (D abla C)$ in 2D/3D | Finite element methods |

Common Dopant Properties in Silicon:

Source: ChipFoundryServices — Search this topic — Ask CFSGPT