Mathematical Modeling of Diffusion in…

Home› Knowledge Base› Mathematical Modeling of Diffusion in Semiconductor Manufacturing

Mathematical Modeling of Diffusion in Semiconductor Manufacturing

1. Fundamental Governing Equations

1.1 Fick's Laws of Diffusion

The foundation of diffusion modeling in semiconductor manufacturing rests on Fick's laws:

Fick's First Law

The flux is proportional to the concentration gradient:

J = -D \frac{\partial C}{\partial x}

Where:

$J$ = flux (atoms/cm²·s)
$D$ = diffusion coefficient (cm²/s)
$C$ = concentration (atoms/cm³)
$x$ = position (cm)

Note: The negative sign indicates diffusion occurs from high to low concentration regions.

Fick's Second Law

Derived from the continuity equation combined with Fick's first law:

\frac{\partial C}{\partial t} = D \frac{\partial^2 C}{\partial x^2}

Key characteristics:

This is a parabolic partial differential equation
Mathematically identical to the heat equation
Assumes constant diffusion coefficient $D$

1.2 Temperature Dependence (Arrhenius Relationship)

The diffusion coefficient follows the Arrhenius relationship:

D(T) = D_0 \exp\left(-\frac{E_a}{kT}\right)

Where:

$D_0$ = pre-exponential factor (cm²/s)
$E_a$ = activation energy (eV)
$k$ = Boltzmann constant ($8.617 \times 10^{-5}$ eV/K)
$T$ = absolute temperature (K)

1.3 Typical Dopant Parameters in Silicon

Dopant	$D_0$ (cm²/s)	$E_a$ (eV)	$D$ at 1100°C (cm²/s)
Boron (B)	~10.5	~3.69	~$10^{-13}$
Phosphorus (P)	~10.5	~3.69	~$10^{-13}$
Arsenic (As)	~0.32	~3.56	~$10^{-14}$
Antimony (Sb)	~5.6	~3.95	~$10^{-14}$

2. Analytical Solutions for Standard Boundary Conditions

2.1 Constant Surface Concentration (Predeposition)

Boundary and Initial Conditions

$C(0,t) = C_s$ — surface held at solid solubility
$C(x,0) = 0$ — initially undoped wafer
$C(\infty,t) = 0$ — semi-infinite substrate

Solution: Complementary Error Function Profile

C(x,t) = C_s \cdot \text{erfc}\left(\frac{x}{2\sqrt{Dt}}\right)

Where the complementary error function is defined as:

\text{erfc}(\eta) = 1 - \text{erf}(\eta) = 1 - \frac{2}{\sqrt{\pi}}\int_0^\eta e^{-u^2} \, du

Total Dose Introduced

Q = \int_0^\infty C(x,t) \, dx = \frac{2 C_s \sqrt{Dt}}{\sqrt{\pi}} \approx 1.13 \, C_s \sqrt{Dt}

Key Properties

Surface concentration remains constant at $C_s$
Profile penetrates deeper with increasing $\sqrt{Dt}$
Characteristic diffusion length: $L_D = 2\sqrt{Dt}$

2.2 Fixed Dose / Gaussian Drive-in

Boundary and Initial Conditions

Total dose $Q$ is conserved (no dopant enters or leaves)
Zero flux at surface: $\left.\frac{\partial C}{\partial x}\right|_{x=0} = 0$
Delta-function or thin layer initial condition

Solution: Gaussian Profile

C(x,t) = \frac{Q}{\sqrt{\pi Dt}} \exp\left(-\frac{x^2}{4Dt}\right)

Time-Dependent Surface Concentration

C_s(t) = C(0,t) = \frac{Q}{\sqrt{\pi Dt}}

Key characteristics:

Surface concentration decreases with time as $t^{-1/2}$
Profile broadens while maintaining total dose
Peak always at surface ($x = 0$)

2.3 Junction Depth Calculation

The junction depth $x_j$ is the position where dopant concentration equals background concentration $C_B$:

For erfc Profile

x_j = 2\sqrt{Dt} \cdot \text{erfc}^{-1}\left(\frac{C_B}{C_s}\right)

For Gaussian Profile

x_j = 2\sqrt{Dt \cdot \ln\left(\frac{Q}{C_B \sqrt{\pi Dt}}\right)}

3. Green's Function Method

3.1 General Solution for Arbitrary Initial Conditions

For an arbitrary initial profile $C_0(x')$, the solution is a convolution with the Gaussian kernel (Green's function):

C(x,t) = \int_{-\infty}^{\infty} C_0(x') \cdot \frac{1}{2\sqrt{\pi Dt}} \exp\left(-\frac{(x-x')^2}{4Dt}\right) dx'

Physical interpretation:

Each point in the initial distribution spreads as a Gaussian
The final profile is the superposition of all spreading contributions

3.2 Application: Ion-Implanted Gaussian Profile

Initial Implant Profile

C_0(x) = \frac{Q}{\sqrt{2\pi} \, \Delta R_p} \exp\left(-\frac{(x - R_p)^2}{2 \Delta R_p^2}\right)

Where:

$Q$ = implanted dose (atoms/cm²)
$R_p$ = projected range (mean depth)
$\Delta R_p$ = straggle (standard deviation)

Profile After Diffusion

C(x,t) = \frac{Q}{\sqrt{2\pi \, \sigma_{eff}^2}} \exp\left(-\frac{(x - R_p)^2}{2 \sigma_{eff}^2}\right)

Effective Straggle

\sigma_{eff} = \sqrt{\Delta R_p^2 + 2Dt}

Key observations:

Peak remains at $R_p$ (no shift in position)
Peak concentration decreases
Profile broadens symmetrically

4. Concentration-Dependent Diffusion

4.1 Nonlinear Diffusion Equation

At high dopant concentrations (above intrinsic carrier concentration $n_i$), diffusion becomes concentration-dependent:

\frac{\partial C}{\partial t} = \frac{\partial}{\partial x}\left(D(C) \frac{\partial C}{\partial x}\right)

4.2 Concentration-Dependent Diffusivity Models

Simple Power Law Model

D(C) = D^i \left(1 + \left(\frac{C}{n_i}\right)^r\right)

Charged Defect Model (Fair's Equation)

D = D^0 + D^- \frac{n}{n_i} + D^{=} \left(\frac{n}{n_i}\right)^2 + D^+ \frac{p}{n_i}

Where:

$D^0$ = neutral defect contribution
$D^-$ = singly negative defect contribution
$D^{=}$ = doubly negative defect contribution
$D^+$ = positive defect contribution
$n, p$ = electron and hole concentrations

4.3 Electric Field Enhancement

High concentration gradients create internal electric fields that enhance diffusion:

J = -D \frac{\partial C}{\partial x} - \mu C \mathcal{E}

For extrinsic conditions with a single dopant species:

J = -hD \frac{\partial C}{\partial x}

Field enhancement factor:

h = 1 + \frac{C}{n + p}

For fully ionized n-type dopant at high concentration: $h \approx 2$
Results in approximately 2× faster effective diffusion

4.4 Resulting Profile Shapes

Phosphorus: "Kink-and-tail" profile at high concentrations
Arsenic: Box-like profiles due to clustering
Boron: Enhanced tail diffusion in oxidizing ambient

5. Point Defect-Mediated Diffusion

5.1 Diffusion Mechanisms

Dopants don't diffuse as isolated atoms—they move via defect complexes:

Vacancy Mechanism

A + V \rightleftharpoons AV \quad \text{(dopant-vacancy pair forms, diffuses, dissociates)}

Interstitial Mechanism

A + I \rightleftharpoons AI \quad \text{(dopant-interstitial pair)}

Kick-out Mechanism

A_s + I \rightleftharpoons A_i \quad \text{(substitutional ↔ interstitial)}

5.2 Effective Diffusivity

D_{eff} = D_V \frac{C_V}{C_V^*} + D_I \frac{C_I}{C_I^*}

Where:

$D_V, D_I$ = diffusivity via vacancy/interstitial mechanism
$C_V, C_I$ = actual vacancy/interstitial concentrations
$C_V^, C_I^$ = equilibrium concentrations

Fractional interstitialcy:

f_I = \frac{D_I}{D_V + D_I}

Dopant	$f_I$	Dominant Mechanism
Boron	~1.0	Interstitial
Phosphorus	~0.9	Interstitial
Arsenic	~0.4	Mixed
Antimony	~0.02	Vacancy

5.3 Coupled Reaction-Diffusion System

The full model requires solving coupled PDEs:

Dopant Equation

\frac{\partial C_A}{\partial t} = \nabla \cdot \left(D_A \frac{C_I}{C_I^*} \nabla C_A\right)

Interstitial Balance

\frac{\partial C_I}{\partial t} = D_I \nabla^2 C_I + G - k_{IV}\left(C_I C_V - C_I^* C_V^*\right)

Vacancy Balance

\frac{\partial C_V}{\partial t} = D_V \nabla^2 C_V + G - k_{IV}\left(C_I C_V - C_I^* C_V^*\right)

Where:

$G$ = defect generation rate
$k_{IV}$ = bulk recombination rate constant

5.4 Transient Enhanced Diffusion (TED)

After ion implantation, excess interstitials cause anomalously rapid diffusion:

The "+1" Model:

\int_0^\infty (C_I - C_I^*) \, dx \approx \Phi \quad \text{(implant dose)}

Enhancement factor:

\frac{D_{eff}}{D^*} = \frac{C_I}{C_I^*} \gg 1 \quad \text{(transient)}

Key characteristics:

Enhancement decays as interstitials recombine
Time constant: typically 10-100 seconds at 1000°C
Critical for shallow junction formation

6. Oxidation Effects

6.1 Oxidation-Enhanced Diffusion (OED)

During thermal oxidation, silicon interstitials are injected into the substrate:

\frac{C_I}{C_I^*} = 1 + A \left(\frac{dx_{ox}}{dt}\right)^n

Effective diffusivity:

D_{eff} = D^* \left[1 + f_I \left(\frac{C_I}{C_I^*} - 1\right)\right]

Dopants enhanced by oxidation:

Boron (high $f_I$)
Phosphorus (high $f_I$)

6.2 Oxidation-Retarded Diffusion (ORD)

Growing oxide absorbs vacancies, reducing vacancy concentration:

\frac{C_V}{C_V^*} < 1

Dopants retarded by oxidation:

Antimony (low $f_I$, primarily vacancy-mediated)

6.3 Segregation at SiO₂/Si Interface

Dopants redistribute at the interface according to the segregation coefficient:

m = \frac{C_{Si}}{C_{SiO_2}}\bigg|_{\text{interface}}

Dopant	Segregation Coefficient $m$	Behavior
Boron	~0.3	Pile-down (into oxide)
Phosphorus	~10	Pile-up (into silicon)
Arsenic	~10	Pile-up

7. Numerical Methods

7.1 Finite Difference Method

Discretize space and time on grid $(x_i, t^n)$:

Explicit Scheme (FTCS)

\frac{C_i^{n+1} - C_i^n}{\Delta t} = D \frac{C_{i+1}^n - 2C_i^n + C_{i-1}^n}{(\Delta x)^2}

Rearranged:

C_i^{n+1} = C_i^n + \alpha \left(C_{i+1}^n - 2C_i^n + C_{i-1}^n\right)

Where Fourier number:

\alpha = \frac{D \Delta t}{(\Delta x)^2}

Stability requirement (von Neumann analysis):

\alpha \leq \frac{1}{2}

Implicit Scheme (BTCS)

\frac{C_i^{n+1} - C_i^n}{\Delta t} = D \frac{C_{i+1}^{n+1} - 2C_i^{n+1} + C_{i-1}^{n+1}}{(\Delta x)^2}

Unconditionally stable (no restriction on $\alpha$)
Requires solving tridiagonal system at each time step

Crank-Nicolson Scheme (Second-Order Accurate)

C_i^{n+1} - C_i^n = \frac{\alpha}{2}\left[(C_{i+1}^{n+1} - 2C_i^{n+1} + C_{i-1}^{n+1}) + (C_{i+1}^n - 2C_i^n + C_{i-1}^n)\right]

Properties:

Unconditionally stable
Second-order accurate in both space and time
Results in tridiagonal system: solved by Thomas algorithm

7.2 Handling Concentration-Dependent Diffusion

Use iterative methods:

1. Estimate $D^{(k)}$ from current concentration $C^{(k)}$ 2. Solve linear diffusion equation for $C^{(k+1)}$ 3. Update diffusivity: $D^{(k+1)} = D(C^{(k+1)})$ 4. Iterate until $\|C^{(k+1)} - C^{(k)}\| < \epsilon$

7.3 Moving Boundary Problems

For oxidation with moving Si/SiO₂ interface:

Approaches:

Coordinate transformation: Map to fixed domain via $\xi = x/s(t)$
Front-tracking methods: Explicitly track interface position
Level-set methods: Implicit interface representation
Phase-field methods: Diffuse interface approximation

8. Thermal Budget Concept

8.1 The Dt Product

Diffusion profiles scale with $\sqrt{Dt}$. The thermal budget quantifies total diffusion:

(Dt)_{total} = \sum_i D(T_i) \cdot t_i

8.2 Continuous Temperature Profile

For time-varying temperature:

(Dt)_{eff} = \int_0^{t_{total}} D(T(\tau)) \, d\tau

8.3 Equivalent Time at Reference Temperature

t_{eq} = \sum_i t_i \exp\left(\frac{E_a}{k}\left(\frac{1}{T_{ref}} - \frac{1}{T_i}\right)\right)

8.4 Combining Multiple Diffusion Steps

For sequential Gaussian redistributions:

\sigma_{final} = \sqrt{\sum_i 2D_i t_i}

For erfc profiles, use effective $(Dt)_{total}$:

C(x) = C_s \cdot \text{erfc}\left(\frac{x}{2\sqrt{(Dt)_{total}}}\right)

9. Key Dimensionless Parameters

Parameter	Definition	Physical Meaning
Fourier Number	$Fo = \dfrac{Dt}{L^2}$	Diffusion time vs. characteristic length
Damköhler Number	$Da = \dfrac{kL^2}{D}$	Reaction rate vs. diffusion rate
Péclet Number	$Pe = \dfrac{vL}{D}$	Advection (drift) vs. diffusion
Biot Number	$Bi = \dfrac{hL}{D}$	Surface transfer vs. bulk diffusion

10. Process Simulation Software

10.1 Commercial and Research Tools

Simulator	Developer	Key Capabilities
Sentaurus Process	Synopsys	Full 3D, atomistic KMC, advanced models
Athena	Silvaco	Integrated with device simulation (Atlas)
SUPREM-IV	Stanford	Classic 1D/2D, widely validated
FLOOPS	U. Florida	Research-oriented, extensible
Victory Process	Silvaco	Modern 3D process simulation

10.2 Physical Models Incorporated

Multiple coupled dopant species
Full point-defect dynamics (I, V, clusters)
Stress-dependent diffusion
Cluster nucleation and dissolution
Atomistic kinetic Monte Carlo (KMC) options
Quantum corrections for ultra-shallow junctions

Mathematical Modeling Hierarchy

Level 1: Simple Analytical Models

\frac{\partial C}{\partial t} = D \frac{\partial^2 C}{\partial x^2}

Constant $D$
erfc and Gaussian solutions
Junction depth calculations

Level 2: Intermediate Complexity

\frac{\partial C}{\partial t} = \frac{\partial}{\partial x}\left(D(C) \frac{\partial C}{\partial x}\right)

Concentration-dependent $D$
Electric field effects
Nonlinear PDEs requiring numerical methods

Level 3: Advanced Coupled Models

\begin{aligned} \frac{\partial C_A}{\partial t} &= \nabla \cdot \left(D_A \frac{C_I}{C_I^*} \nabla C_A\right) \\[6pt] \frac{\partial C_I}{\partial t} &= D_I \nabla^2 C_I + G - k_{IV}(C_I C_V - C_I^* C_V^*) \end{aligned}

Coupled dopant-defect systems
TED, OED/ORD effects
Process simulators required

Level 4: State-of-the-Art

Atomistic kinetic Monte Carlo
Molecular dynamics for interface phenomena
Ab initio calculations for defect properties
Essential for sub-10nm technology nodes

Key Insight

The fundamental scaling of semiconductor diffusion is governed by $\sqrt{Dt}$, but the effective diffusion coefficient $D$ depends on:

Temperature (Arrhenius)
Concentration (charged defects)
Point defect supersaturation (TED)
Processing ambient (oxidation)
Mechanical stress

This complexity requires sophisticated physical models for modern nanometer-scale devices.

diffusion modelingdiffusion modelfick law modelingdopant diffusion modelsemiconductor diffusion modelthermal diffusion modeldiffusion coefficient calculationdiffusion simulationdiffusion mathematics

Explore 500+ Semiconductor & AI Topics

From EUV lithography to CUDA optimization — search the full knowledge base or chat with our AI assistant.

🔍 Search Topics 💬 Ask CFSGPT 📚 Browse All

Related Topics

Explore 500+ Semiconductor & AI Topics