fokker planck stochastic advection diffusion

# Fokker-Planck Equations and Stochastic Dynamics in Ion Implantation

## 1. Introduction to Stochastic Processes in Materials Engineering

The evolution of defects, dopants, and point defects in semiconductors under thermal annealing, ion implantation, and irradiation is fundamentally stochastic. Individual defects perform random walks driven by thermal diffusion and electric fields, with occasional reactions (recombination, clustering). When tracking populations of millions of particles, the Fokker-Planck equation replaces detailed stochastic trajectories with a deterministic PDE for the probability density function (PDF).

The Fokker-Planck equation is the Smoluchowski limit of the Langevin SDE, valid when friction dominates inertia—exactly the regime for point defect and dopant dynamics in crystalline solids.

## 2. Langevin Equation and Its Ensemble Behavior

Consider a single particle at position $x(t)$ in a potential $U(x)$ and thermal bath at temperature $T$. The Langevin equation:

$$m \frac{d^2x}{dt^2} = -\gamma \frac{dx}{dt} - \frac{dU}{dx} + \xi(t)$$

describes overdamped motion with:
- Friction coefficient: $\gamma$ (units: s⁻¹)
- Conservative force: $-dU/dx$
- Thermal noise: $\xi(t)$ (Gaussian, zero-mean, $\langle\xi(t)\xi(t')
angle = 2\gamma k_B T \delta(t-t')$)

The second term drives the particle downhill in potential; the noise term causes random excursions.

For semiconductors, the Smoluchowski limit ($m o 0$, $\gamma$ remains finite) gives:

$$\frac{dx}{dt} = -\frac{1}{\gamma}\frac{dU}{dx} + \sqrt{\frac{2k_B T}{\gamma}} \eta(t)$$

where $\eta(t)$ is normalized Gaussian white noise.

## 3. Derivation of the Fokker-Planck Equation

Starting from the Smoluchowski SDE, the Fokker-Planck equation for the probability density $P(x,t)$ is:

$$\frac{\partial P}{\partial t} = -\frac{\partial}{\partial x}\left[A(x)P ight] + \frac{\partial^2}{\partial x^2}\left[B(x)P ight]$$

where:
- Drift coefficient: $A(x) = -\frac{1}{\gamma}\frac{dU}{dx}$
- Diffusion coefficient: $B(x) = \frac{k_B T}{\gamma}$

This is the continuity equation for probability flux. Physically:
- The drift term $\partial_x[AP]$ represents deterministic motion along gradients
- The diffusion term $\partial_x^2[BP]$ spreads the probability cloud due to thermal noise
- The balance of these governs the kinetics of entire populations

## 4. Kramers-Moyal Expansion

The Fokker-Planck equation emerges from a systematic Kramers-Moyal expansion of the master equation for jump processes. For a particle with transition rates $W(x o x')$:

$$P(x,t+\Delta t) = \int dx' W(x' o x) P(x',t)$$

Expanding $P(x',t)$ around $x$ to second order:

$$\frac{\partial P}{\partial t} = \int dx' W(x' o x) \left[P(x,t) - (x'-x)\frac{\partial P}{\partial x} + \frac{(x'-x)^2}{2}\frac{\partial^2 P}{\partial x^2} ight]$$

Collecting terms:
- First-order moment: $m_1(x) = \int (x'-x) W(x' o x) dx'$ → drift coefficient $A(x) = -m_1(x)$
- Second-order moment: $m_2(x) = \int (x'-x)^2 W(x' o x) dx'$ → diffusion coefficient $B(x) = m_2(x)/2$

Higher-order moments are dropped in the Markovian approximation (memory-less), valid when the system correlation time is short.

## 5. Stationary Solutions and Einstein Relation

For a system approaching equilibrium (no sources/sinks), the probability converges to the Boltzmann distribution:

$$P_{eq}(x) \propto \exp\left(-\frac{U(x)}{k_B T} ight)$$

At stationarity, $\partial_t P = 0$, so the probability flux must vanish:

$$J(x) = A(x)P_{eq} - \frac{\partial}{\partial x}[B(x)P_{eq}] = 0$$

Substituting $P_{eq} = \exp(-U/k_B T)$ and requiring $J=0$:

$$\frac{A(x)}{B(x)} = \frac{1}{k_B T}\frac{dU}{dx}$$

This is the Einstein relation: $D = k_B T / \gamma$, connecting mobility $\mu = 1/\gamma$ to diffusivity $D$.

For point defects in Si:
$$D = \mu_0 \exp\left(-\frac{E_a}{k_B T} ight)$$

where $E_a$ is the activation energy (migration barrier) and $\mu_0$ is the pre-exponential factor (~10⁻⁷ cm²/Vs for interstitials).

## 6. Defect Types and Drift-Diffusion Coefficients

### Intrinsic Point Defects (Si)

Diffusivity expressions:
$$D_V = D_{V,0} \exp(-0.6 ext{ eV}/k_B T) \quad ext{[cm}^2 ext{/s]}$$
$$D_I = D_{I,0} \exp(-0.4 ext{ eV}/k_B T)$$

For charged defects in electric field $E_r$:
$$D_{ ext{eff}} = D \left[1 - \frac{qE_r \cdot l}{k_B T} ight]$$

where $l$ is the hopping distance and $qE_r$ is the potential energy shift.

## 7. Drift-Diffusion with External Fields

Under applied electric field $E_r$ (e.g., from implantation damage, doping gradients), the Fokker-Planck equation becomes:

$$\frac{\partial n}{\partial t} = abla \cdot [D(x) abla n - \mu(x) n abla \phi]$$

where $\phi(x)$ is the electric potential and $\mu = qD/k_B T$ is the drift mobility.

For a 1D channel with constant field:
$$\frac{\partial n}{\partial t} = D \frac{\partial^2 n}{\partial x^2} - \mu E_r \frac{\partial n}{\partial x}$$

The Peclet number $Pe = \mu E_r L / D$ compares drift to diffusion:
- $Pe \ll 1$: Diffusion-dominated, symmetric spreading
- $Pe \gg 1$: Drift-dominated, asymmetric toward cathode

## 8. Transient Enhanced Diffusion (TED) and Extrinsic Defect Reactions

TED is the anomalous acceleration of dopant diffusion following ion implantation. The mechanism:

1. Implantation creates supersaturation of interstitials (via Frenkel pair generation)
2. Interstitials diffuse rapidly (lower $E_m$ than vacancies)
3. Interstitials "drag along" dopants via kick-out reactions: $D_s + I o D_i^{+} + V$
4. Dopant velocity increases by factor 10-100× for ~1 hour anneal

The supersaturation $C_I(t)$ of interstitials decays via:
$$\frac{dC_I}{dt} = G - C_I C_V / au_r$$

where $G$ is implant generation rate, $ au_r \sim 1-10$ ns is recombination time, and $C_V$ is vacancy concentration.

Dopant diffusion under TED:
$$\frac{\partial C_D}{\partial t} = D_0 \frac{\partial^2 C_D}{\partial x^2} + \alpha(t) \frac{\partial C_D}{\partial x}$$

where $\alpha(t) = k_{kick} C_I(t)$ is the kick-out drift coefficient, decaying exponentially as interstitials recombine.

## 9. End-of-Range (EOR) Defect Clustering

Beyond the implanted region, end-of-range defect clusters form via:
1. Cascade-induced damage creates high local defect density
2. Vacancies agglomerate into dislocation loops (Frank loops)
3. Subsequent annealing can dissolve or grow loops

The loop growth kinetics:
$$\frac{dR}{dt} = D_V \left[\frac{C_V - C_V^{eq}}{C_V^{eq}} ight] \left[\frac{1}{R} ight]$$

Large loops (R > 10 nm) are stable; small loops shrink due to line tension.

## 10. Python Implementation: 1D Fokker-Planck Solver for Dopant Diffusion

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.ndimage import convolve1d

def fokker_planck_1d_tdma(n_prev, D, mu_E, dt, dx, x_array):
    """
    Solve 1D Fokker-Planck via Tridiagonal Matrix Algorithm (TDMA).
    
    ∂n/∂t = D ∂²n/∂x² - μE ∂n/∂x
    
    Using Crank-Nicolson implicit scheme for stability.
    
    Parameters:
    -----------
    n_prev : array
        Concentration at previous time (m^-3)
    D : float or array
        Diffusion coefficient (m²/s)
    mu_E : float or array
        Drift coefficient μE (m/s)
    dt : float
        Time step (s)
    dx : float
        Spatial step (m)
    x_array : array
        Spatial grid
    
    Returns:
    --------
    n_new : array
        Concentration at new time
    """
    
    N = len(n_prev)
    
    # Courant numbers
    r_diff = D * dt / dx**2
    r_drift = mu_E * dt / (2 * dx)
    
    # Tridiagonal coefficients (Crank-Nicolson)
    # Main diagonal: 1 + 2r_diff
    # Upper/lower: -r_diff ± r_drift
    
    a = -r_diff + r_drift  # Lower diagonal
    b = 1.0 + 2*r_diff * np.ones(N)  # Main diagonal
    c = -r_diff - r_drift  # Upper diagonal
    
    # RHS vector
    rhs = n_prev.copy()
    rhs[1:-1] += r_diff * (n_prev[2:] - 2*n_prev[1:-1] + n_prev[:-2])
    rhs[1:-1] += r_drift * (n_prev[2:] - n_prev[:-2])
    
    # TDMA forward elimination
    b_mod = b.copy()
    rhs_mod = rhs.copy()
    
    for i in range(1, N):
        factor = a / b_mod[i-1]
        b_mod[i] -= factor * c
        rhs_mod[i] -= factor * rhs_mod[i-1]
    
    # Back substitution
    n_new = np.zeros(N)
    n_new[-1] = rhs_mod[-1] / b_mod[-1]
    
    for i in range(N-2, -1, -1):
        n_new[i] = (rhs_mod[i] - c * n_new[i+1]) / b_mod[i]
    
    return n_new

def ted_interstitial_kinetics(t_array, C_I_0=1e20, tau_r=2e-9, G_rate=1e27):
    """
    Compute interstitial supersaturation decay.
    
    dC_I/dt = G - C_I² / τ_r (bimolecular recombination)
    
    Parameters:
    -----------
    t_array : array (s)
        Time array
    C_I_0 : float (m^-3)
        Initial interstitial concentration
    tau_r : float (s)
        Recombination time constant
    G_rate : float (m^-3·s^-1)
        Generation rate (during implant)
    
    Returns:
    --------
    C_I : array (m^-3)
        Interstitial concentration vs time
    """
    
    def dCdt(C_I, t):
        # Generation only during first 1 ns (implant pulse)
        G = G_rate if t < 1e-9 else 0
        recombination = C_I**2 / tau_r  # Second-order kinetics
        return G - recombination
    
    C_I = odeint(dCdt, C_I_0, t_array)
    return C_I.flatten()

def ted_dopant_diffusion_simulation(x_array, t_array, C_D_init, D_0, E_a, T=1073):
    """
    Simulate dopant diffusion under Transient Enhanced Diffusion (TED).
    
    Parameters:
    -----------
    x_array : array (m)
        Spatial grid (1D)
    t_array : array (s)
        Time grid (annealing duration)
    C_D_init : array (m^-3)
        Initial dopant profile (e.g., delta function at surface)
    D_0 : float (m²/s)
        Pre-exponential diffusion coefficient
    E_a : float (eV)
        Activation energy
    T : float (K)
        Annealing temperature
    
    Returns:
    --------
    C_D_history : array (N_space, N_time)
        Dopant concentration evolution
    """
    
    k_B = 8.617e-5  # eV/K
    dx = x_array[1] - x_array[0]
    dt = t_array[1] - t_array[0]
    
    # Diffusion coefficient
    D = D_0 * np.exp(-E_a / (k_B * T))
    
    # Interstitial kinetics (background TED)
    C_I = ted_interstitial_kinetics(t_array, C_I_0=1e21, tau_r=1e-9)
    
    # Kick-out drift coefficient (proportional to I surplus)
    k_kickout = 1e-8  # [dimensionless or m³/s in full model]
    mu_E_array = k_kickout * (C_I - 1e19)  # Relative to equilibrium
    
    # Initialize history
    N_space = len(x_array)
    N_time = len(t_array)
    C_D_history = np.zeros((N_space, N_time))
    C_D_history[:, 0] = C_D_init
    
    # Time-stepping via TDMA
    C_D_current = C_D_init.copy()
    
    for n in range(1, N_time):
        mu_E = mu_E_array[n]
        C_D_new = fokker_planck_1d_tdma(C_D_current, D, mu_E, dt, dx, x_array)
        
        # Ensure non-negativity
        C_D_new = np.maximum(C_D_new, 0)
        
        C_D_history[:, n] = C_D_new
        C_D_current = C_D_new
    
    return C_D_history

def effective_diffusion_length(D, t):
    """
    Compute diffusion length L_D = sqrt(4Dt).
    
    Parameters:
    -----------
    D : float (m²/s)
        Diffusion coefficient
    t : float (s)
        Time
    
    Returns:
    --------
    L_D : float (m)
        Characteristic diffusion length
    """
    return np.sqrt(4 * D * t)

# Simulation parameters
x = np.linspace(0, 1e-6, 500)  # 0-1 μm, 500 points
t = np.logspace(-9, -6, 100)   # 1 ns to 1 μs

# Initial condition: delta-like dopant profile at surface
C_D_init = np.exp(-(x - 50e-9)**2 / (10e-9)**2) * 1e24

# Run TED simulation
C_D_evolution = ted_dopant_diffusion_simulation(
    x, t, C_D_init, D_0=1e-13, E_a=4.0, T=1073
)

# Compute diffusion length at key times
D = 1e-13 * np.exp(-4.0 / (8.617e-5 * 1073))
L_D_1ns = effective_diffusion_length(D, 1e-9)
L_D_1us = effective_diffusion_length(D, 1e-6)

# Plotting
fig, axes = plt.subplots(2, 2, figsize=(13, 10))

# Panel A: Dopant profiles at different anneal times
time_indices = [0, 10, 30, 60, 99]
colors = plt.cm.viridis(np.linspace(0, 1, len(time_indices)))

for idx, i in enumerate(time_indices):
    axes[0, 0].semilogy(x*1e6, C_D_evolution[:, i] + 1e15, 'o-', 
                        color=colors[idx], label=f't={t[i]*1e9:.0f} ns', 
                        markersize=3, alpha=0.7)

axes[0, 0].set_xlabel('Depth (μm)', fontsize=11)
axes[0, 0].set_ylabel('Dopant Conc. (cm⁻³)', fontsize=11)
axes[0, 0].set_title('Dopant Diffusion Profiles (TED)', fontsize=12, fontweight='bold')
axes[0, 0].legend(fontsize=9)
axes[0, 0].grid(True, alpha=0.3)

# Panel B: Interstitial dynamics
C_I = ted_interstitial_kinetics(t, C_I_0=1e21)
axes[0, 1].semilogy(t*1e9, C_I, 'r-', linewidth=2)
axes[0, 1].set_xlabel('Time (ns)', fontsize=11)
axes[0, 1].set_ylabel('Interstitial Conc. (cm⁻³)', fontsize=11)
axes[0, 1].set_title('Interstitial Recombination Kinetics', fontsize=12, fontweight='bold')
axes[0, 1].grid(True, alpha=0.3, which='both')

# Panel C: Contour plot of dopant evolution
X_plot, T_plot = np.meshgrid(x*1e6, t*1e9)
levels = np.logspace(15, 24, 20)
contour = axes[1, 0].contourf(X_plot, T_plot, C_D_evolution.T, levels=levels, cmap='plasma')
axes[1, 0].set_xlabel('Depth (μm)', fontsize=11)
axes[1, 0].set_ylabel('Anneal Time (ns)', fontsize=11)
axes[1, 0].set_title('Dopant Evolution (Space-Time)', fontsize=12, fontweight='bold')
axes[1, 0].set_yscale('log')
plt.colorbar(contour, ax=axes[1, 0], label='Conc. (cm⁻³)')

# Panel D: Total integrated dopant amount
total_dopant = np.array([np.trapz(C_D_evolution[:, i], x) for i in range(len(t))])
axes[1, 1].plot(t*1e9, total_dopant/total_dopant[0], 'g-o', linewidth=2, markersize=5)
axes[1, 1].set_xlabel('Anneal Time (ns)', fontsize=11)
axes[1, 1].set_ylabel('Normalized Total Dopant', fontsize=11)
axes[1, 1].set_title('Dopant Conservation Check', fontsize=12, fontweight='bold')
axes[1, 1].set_xscale('log')
axes[1, 1].grid(True, alpha=0.3)
axes[1, 1].set_ylim([0.95, 1.05])

plt.tight_layout()
plt.savefig('ted_fokker_planck.png', dpi=150, bbox_inches='tight')
plt.show()

print("=== Fokker-Planck TED Simulation Complete ===")
print(f"Characteristic diffusion length @ 1 ns: {L_D_1ns*1e9:.2f} nm")
print(f"Characteristic diffusion length @ 1 μs: {L_D_1us*1e9:.2f} nm")
print(f"Dopant conservation: {total_dopant[-1]/total_dopant[0]:.4f}")

## 11. Multi-dimensional Fokker-Planck in 2D and 3D

For spatial processes in 2D/3D (e.g., EOR loop growth, precipitate nucleation):

$$\frac{\partial n(\vec{r},t)}{\partial t} = abla \cdot [D(\vec{r}) abla n] - abla \cdot [\vec{v}(\vec{r}) n]$$

where $\vec{v}(\vec{r}) = -
abla U / \gamma$ is the drift velocity field.

In cylindrical geometry (suitable for dislocation loops):
$$\frac{\partial n}{dt} = \frac{1}{r}\frac{\partial}{\partial r}\left[rD\frac{\partial n}{\partial r} ight] + \frac{\partial}{\partial z}\left[D\frac{\partial n}{\partial z} ight] - ext{drift terms}$$

Numerical solution via alternating direction implicit (ADI) scheme or operator splitting (fractional steps in r and z directions).

## 12. Reaction-Diffusion Systems: Clustering and Coarsening

When defects interact (e.g., V + I → void, V + V → dislocation loop), the system becomes a reaction-diffusion PDE:

$$\frac{\partial n_V}{\partial t} = D_V abla^2 n_V - k_r n_V n_I + k_f$$
$$\frac{\partial n_I}{\partial t} = D_I abla^2 n_I - k_r n_V n_I + k_f$$

where $k_r$ is bimolecular recombination rate and $k_f$ is Frenkel pair generation.

Coarsening dynamics: After implant cease ($k_f = 0$), the defect size distribution evolves via Ostwald ripening—large precipitates grow by absorbing smaller ones, following:

$$\bar{R}(t) = R_0 \left(\frac{t}{t_0} ight)^{1/3}$$

This $t^{1/3}$ growth law is universal in diffusion-limited coarsening.

## 13. Comparison with Kinetic Monte Carlo (KMC)

For implant damage: KMC handles cascade damage (ps scale); Fokker-Planck models annealing (hours scale).

## 14. Boundary Conditions and Sources/Sinks

### Dirichlet (Fixed Concentration)
$$n(x_{ ext{boundary}}) = n_0$$
Represents equilibrium with a reservoir (e.g., dopant evaporation from surface).

### Neumann (Fixed Flux)
$$-D \frac{\partial n}{\partial x}\bigg|_{ ext{boundary}} = J_0$$
Represents injection/extraction (e.g., implant ion beam).

### Mixed (Robin)
$$D \frac{\partial n}{\partial x}\bigg|_{ ext{boundary}} = k_{ ext{eff}} [n - n_{ ext{bulk}}]$$
Represents surface recombination with velocity $k_{ ext{eff}}$.

## 15. Advanced Applications

Transient Dopant Diffusion Modeling (TDD) in CMOS: Combines implant profiles (SRIM), TED/EOR kinetics, and final dopant mapping to predict threshold voltage and leakage currents.

Implant Damage Annealing: Focal Defect + Fokker-Planck couples KMC cascade damage to diffusion-driven loop dissolution and defect clustering over extended annealing schedules.

Boron Point-Defect Interactions: In B-implanted Si, kick-out reactions dominate over intrinsic diffusion, reversible B-I complex formation, and transient supersaturation.

The Fokker-Planck formalism remains essential for process simulation in semiconductor manufacturing, bridging atomistic understanding and process-level predictions.