Mathematical Modeling of Diffusion

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

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} $$

5.3 Coupled Reaction-Diffusion System

The full model requires solving coupled PDEs :

Dopant Equation

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

Interstitial Balance

$$ \frac{\partial C_I}{\partial t} = D_I abla^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 abla^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}} $$

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

10. Process Simulation Software

10.1 Commercial and Research Tools

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} &= abla \cdot \left(D_A \frac{C_I}{C_I^} abla C_A\right) \\[6pt] \frac{\partial C_I}{\partial t} &= D_I abla^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.