Part 6, Chapter 6

Plasma-Wall Interaction

Sputtering, sheath physics, Bohm criterion, and tritium retention

6.1 Physical Sputtering

When energetic ions strike a solid surface, they transfer momentum to target atoms through a collision cascade. If a surface atom receives enough energy to overcome its surface binding energy U_s, it is ejected -- this is physical sputtering. The sputtering yield Y(E) (atoms removed per incident ion) is described by the Bohdansky formula:

$$Y(E) = Q\,\frac{s_n(\epsilon)}{1 + \Gamma\,k_e\,\epsilon^{0.3}}\left(1 - \left(\frac{E_{th}}{E}\right)^{2/3}\right)\left(1 - \frac{E_{th}}{E}\right)^2$$

where epsilon = E / E_TF is the reduced energy (E_TF is the Thomas-Fermi energy), s_n(epsilon) is the nuclear stopping cross-section, k_e is the Lindhard electronic stopping coefficient, Q is a yield factor, and Gamma is a dimensionless constant. The threshold energy E_th below which no sputtering occurs is approximately:

$$E_{th} = \frac{U_s}{\gamma(1 - \gamma)} \quad \text{where} \quad \gamma = \frac{4 M_1 M_2}{(M_1 + M_2)^2}$$

The energy transfer factor gamma determines the maximum fraction of kinetic energy transferred in a head-on collision. For light ions (D, T) on heavy targets (W), gamma is small, making E_th large and Y small. This is why tungsten (W) is the primary choice for divertor plasma-facing components: its high mass (M = 184) gives a high sputtering threshold of about 200 eV for deuterium and its high melting point (3422 C) provides resilience against transient heat loads.

At typical divertor conditions (D+ at 50-300 eV on W), the sputtering yield is Y approximately 10^-3 to 10^-2, meaning roughly one tungsten atom is eroded per 100-1000 incident ions. For carbon (formerly used in many tokamaks), Y is 10-100x higher, plus chemical sputtering by hydrogen produces volatile hydrocarbons even at room temperature.

6.2 The Plasma Sheath

When a plasma contacts a material surface, electrons (being lighter and faster) initially escape more rapidly than ions, charging the surface negatively. A positively charged sheath region of a few Debye lengths forms to equalize the fluxes. The sheath potential drop is:

$$\phi_s = -\frac{T_e}{2e}\ln\!\left(\frac{2\pi m_e}{m_i}\right)$$

For a deuterium plasma, this gives phi_s approximately -3.0 T_e/e. The sheath accelerates ions to energies of about 3 T_e before they hit the surface, while repelling all but the most energetic electrons. The energy flux transmitted through the sheath to the surface is:

$$q_{wall} = \gamma_{sh}\,n_s\,c_s\,T_e \quad \text{with} \quad \gamma_{sh} = \frac{2}{1-\delta_{se}} + \frac{2T_i}{T_e} - \frac{e\phi_s}{T_e}$$

where delta_se is the secondary electron emission coefficient and gamma_sh approximately 7-8 for typical conditions. The sheath is collisionless (thickness approximately 5 lambda_D approximately 0.05 mm) and represents the boundary condition for all plasma-surface interaction calculations.

6.3 The Bohm Sheath Criterion

For a stable sheath to exist, ions must enter the sheath with a minimum velocity. This fundamental requirement is the Bohm sheath criterion:

$$v_i \geq c_s = \sqrt{\frac{k_B(T_e + T_i)}{m_i}}$$

The derivation starts from Poisson's equation in the sheath. For a monotonically decreasing potential (repelling electrons, accelerating ions), the ion density must decrease slower than the electron density as the potential drops. Assuming Boltzmann electrons n_e = n_0 exp(e phi / T_e) and cold ions accelerated from rest through potential phi, the ion density is n_i = n_0 (1 - 2e phi / m_i v_0^2)^(-1/2). Requiring d^2 phi / dx^2 less than 0 at the sheath edge yields:

$$\frac{1}{2}m_i v_0^2 \geq \frac{1}{2}k_B T_e \quad \Rightarrow \quad v_0 \geq \sqrt{\frac{k_B T_e}{m_i}}$$

This means there must be a "presheath" region where ions are accelerated from thermal velocities up to the sound speed. The presheath extends over a scale comparable to the ion mean free path or the plasma dimension, and the density drops by a factor of about 2 across it (n_sheath approximately 0.5 n_bulk). The Bohm criterion is universal -- it applies to any plasma-surface interface and sets the particle flux to:

$$\Gamma_i = \frac{1}{2}\,n_0\,c_s$$

6.4 Recycling and Tritium Retention

When hydrogen isotope ions strike a surface, they can be reflected, implanted, or absorbed. The recycling coefficient R is the fraction of incident particles that return to the plasma (as atoms or molecules). For saturated surfaces, R approaches 1.0, meaning nearly all particles are recycled. The reflection coefficient depends on the ion energy and target material -- for deuterium on tungsten at 100 eV, the reflection coefficient is about 0.5.

Tritium retention is a critical safety issue for fusion reactors. Tritium can be retained by: (1) implantation into the near-surface layer, (2) diffusion into the bulk material, (3) trapping at lattice defects, grain boundaries, and voids, and (4) co-deposition with eroded material (especially problematic for carbon). The retention fraction depends on the material temperature, fluence, and surface condition.

For tungsten, the retention fraction is relatively low (approximately 10^-5 to 10^-4 of the incident fluence at 500 K) because tungsten has low hydrogen solubility and diffusivity. ITER has a safety limit of 700 g of in-vessel tritium, which constrains the choice of plasma-facing materials and the required bake-out frequency. Carbon was eliminated from ITER's divertor design precisely because carbon co-deposits with tritium at unacceptable rates.

Interactive Simulations

Sputtering Yield vs Ion Energy for W, C, Be Targets

Python

script.py115 lines

import numpy as np

# Physical Sputtering Yield: Bohdansky/Eckstein Formula
print("Physical Sputtering Yield: D+ on W, C, Be")
print("=" * 55)

e_charge = 1.602e-19
a_Bohr = 5.292e-11  # m

# Material parameters: [Z, M (amu), U_s (eV), Q, lambda]
materials = {
    'W':  {'Z2': 74, 'M2': 183.84, 'U_s': 8.68, 'Q': 0.042, 'label': 'Tungsten'},
    'C':  {'Z2': 6,  'M2': 12.011, 'U_s': 7.41, 'Q': 0.035, 'label': 'Carbon'},
    'Be': {'Z2': 4,  'M2': 9.012,  'U_s': 3.38, 'Q': 0.042, 'label': 'Beryllium'},
}

# Projectile: Deuterium
Z1 = 1
M1 = 2.014

def thomas_fermi_energy(Z1, Z2, M1, M2):
    """Thomas-Fermi screening length and energy"""
    a_TF = 0.4685 * a_Bohr / (Z1**(2/3) + Z2**(2/3))**0.5  # A
    E_TF = Z1 * Z2 * e_charge**2 / (4 * np.pi * 8.854e-12 * a_TF) / e_charge  # eV
    return a_TF, E_TF

def sputtering_yield(E, Z1, M1, Z2, M2, U_s, Q):
    """Bohdansky formula for sputtering yield"""
    # Energy transfer factor
    gamma = 4 * M1 * M2 / (M1 + M2)**2

# Threshold energy
    if M1 / M2 <= 0.2:
        E_th = U_s / (gamma * (1 - gamma))
    else:
        E_th = 8 * U_s * (M1 / M2)**0.4

# Thomas-Fermi reduced energy
    a_TF, E_TF = thomas_fermi_energy(Z1, Z2, M1, M2)
    epsilon = E / E_TF

# Nuclear stopping (Lindhard universal function, KrC approximation)
    s_n = 0.5 * np.log(1 + 1.2288 * epsilon) /           (epsilon + 0.1728 * epsilon**0.5 + 0.008 * epsilon**0.1504)

# Electronic stopping coefficient
    k_e = 0.0793 * (Z1**(2/3) * Z2**(1/2) * (M1 + M2)**(3/2)) /           ((Z1**(2/3) + Z2**(2/3))**(3/4) * M1**(3/2) * M2**(1/2))

# Yield
    Y = np.where(E > E_th,
        Q * s_n / (1 + 0.5 * k_e * epsilon**0.3) *         (1 - (E_th / E)**(2/3)) * (1 - E_th / E)**2,
        0.0)

return Y, E_th

# Energy scan
E = np.logspace(0, 4, 500)  # 1 eV to 10 keV

print()
for mat_name, mat in materials.items():
    Y, E_th = sputtering_yield(E, Z1, M1, mat['Z2'], mat['M2'], mat['U_s'], mat['Q'])
    Y_max = np.max(Y)
    E_max = E[np.argmax(Y)]
    print(f"{mat['label']} ({mat_name}):")
    print(f"  Z = {mat['Z2']}, M = {mat['M2']:.1f} amu, U_s = {mat['U_s']:.2f} eV")
    print(f"  Threshold energy: E_th = {E_th:.1f} eV")
    print(f"  Peak yield: Y_max = {Y_max:.4f} at E = {E_max:.0f} eV")
    print()

# Detailed table
print("Sputtering yield Y(E) for D+ projectiles:")
print(f"{'E (eV)':>8} {'Y(W)':>12} {'Y(C)':>12} {'Y(Be)':>12}")
print("-" * 46)

E_table = [10, 20, 30, 50, 100, 200, 300, 500, 1000, 2000, 5000, 10000]
for Ei in E_table:
    E_arr = np.array([float(Ei)])
    Y_W, _ = sputtering_yield(E_arr, Z1, M1, 74, 183.84, 8.68, 0.042)
    Y_C, _ = sputtering_yield(E_arr, Z1, M1, 6, 12.011, 7.41, 0.035)
    Y_Be, _ = sputtering_yield(E_arr, Z1, M1, 4, 9.012, 3.38, 0.042)
    print(f"{Ei:8d} {Y_W[0]:12.6f} {Y_C[0]:12.6f} {Y_Be[0]:12.6f}")

print()
print("=" * 55)
print("Sheath Potential and Ion Impact Energy")
print("=" * 55)
print()

# Sheath potential for D plasma on different materials
T_e_arr = np.array([5, 10, 20, 50, 100, 200])
m_D = 2.014 * 1.661e-27
m_e = 9.109e-31

print(f"{'T_e (eV)':>10} {'phi_s (V)':>10} {'E_impact (eV)':>14} {'Y(W)':>10} {'Y(C)':>10}")
print("-" * 56)

for T_e in T_e_arr:
    phi_s = -T_e / 2 * np.log(2 * np.pi * m_e / m_D)
    E_impact = -e_charge * phi_s / e_charge + 2 * T_e  # sheath + thermal
    # In eV units
    phi_s_V = T_e / 2 * np.log(2 * np.pi * m_e / m_D)
    E_imp = phi_s_V + 2 * T_e  # total impact energy approx 3*T_e + 2*T_i

E_arr = np.array([E_imp])
    Y_W, _ = sputtering_yield(E_arr, Z1, M1, 74, 183.84, 8.68, 0.042)
    Y_C, _ = sputtering_yield(E_arr, Z1, M1, 6, 12.011, 7.41, 0.035)

print(f"{T_e:10.0f} {phi_s_V:10.1f} {E_imp:14.1f} {Y_W[0]:10.6f} {Y_C[0]:10.6f}")

print()
print("Key findings:")
print("- W sputtering threshold (~200 eV) much higher than C (~30 eV) or Be (~15 eV)")
print("- At detached conditions (T_e < 5 eV, E_impact < 20 eV): W yield ~ 0")
print("- Carbon has additional chemical sputtering (not modeled here) of Y ~ 0.01-0.05")
print("- ITER chose W divertor + Be first wall to minimize erosion and T retention")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Previous Next →