Electron energy distribution

Download notebookView on GitHubRun when the site was built.

Hot electrons in a gas lose energy slowly in elastic collisions and quickly by excitation and ionization, until no inelastic channel is open. This tutorial reads a cross-section set in LXCat format and follows the electron energy distribution function (EEDF) in a 0D reactor.

The LXCat file is written here with synthetic model numbers. With real data downloaded from LXCat, only the file name changes.

from pathlib import Path

import cunumpy as xp
import numpy as np

import plasmacoll
from plasmacoll.constants import (
    ATOMIC_MASS,
    ELECTRON_MASS,
    ELEMENTARY_CHARGE,
    ev_to_kelvin,
)

energy = np.round(np.geomspace(1e-2, 1e3, 60), 4)


def block(kind, target, parameter, sigma, process):
    lines = [kind, target]
    if parameter is not None:
        lines.append(f" {parameter:.4e}")
    lines += [f"PROCESS: {process}", "-----------------------------"]
    lines += [f" {e:.4e}\t{s:.4e}" for e, s in zip(energy, sigma, strict=True)]
    lines.append("-----------------------------")
    return "\n".join(lines)


elastic = 8e-20 / (1 + (energy / 10) ** 0.8)
excitation = np.where(
    energy > 11.5, 2e-21 * (energy - 11.5) / (1 + (energy / 30) ** 1.5), 0
)
x = np.clip(energy / 15.8, 1, None)
ionization = np.where(energy > 15.8, 3e-20 * np.log(x) / x, 0)
Path("model_gas.txt").write_text(
    "\n\n".join(
        [
            "SYNTHETIC MODEL GAS: these numbers are not measured data.",
            block("ELASTIC", "Gas", 1.37e-5, elastic, "E + Gas -> E + Gas, Elastic"),
            block(
                "EXCITATION",
                "Gas -> Gas*",
                11.5,
                excitation,
                "E + Gas -> E + Gas*, Excitation",
            ),
            block(
                "IONIZATION",
                "Gas -> Gas^+",
                15.8,
                ionization,
                "E + Gas -> E + E + Gas+, Ionization",
            ),
        ]
    )
)

table = plasmacoll.read_lxcat("model_gas.txt")
for process in table:
    print(
        f"{process.kind:11s} {process.target:>4s} -> {process.product:5s} threshold {process.threshold:5.1f} eV"
    )

ELASTIC      Gas ->       threshold   0.0 eV
EXCITATION   Gas -> Gas*  threshold  11.5 eV
IONIZATION   Gas -> Gas^+ threshold  15.8 eV

find_lxcat_process picks a block by kind or by its PROCESS comment, and its cross_section goes straight into a CollisionProcess. Electron cross-section sets are tabulated against the electron energy in the gas frame, so the processes use energy_frame="lab".

M = 40 * ATOMIC_MASS
gas = plasmacoll.NeutralBackground("gas", 1e22, 300.0, M)
products = {"electron": "e", "ion": "ion"}
processes = [
    plasmacoll.CollisionProcess(
        "elastic",
        "gas",
        plasmacoll.find_lxcat_process(table, "ELASTIC").cross_section,
        name="elastic",
        energy_frame="lab",
    ),
    plasmacoll.CollisionProcess(
        "excitation",
        "gas",
        plasmacoll.find_lxcat_process(table, "EXCITATION").cross_section,
        name="excitation",
        energy_frame="lab",
    ),
    plasmacoll.CollisionProcess(
        "ionization",
        "gas",
        plasmacoll.find_lxcat_process(table, "IONIZATION").cross_section,
        name="ionization",
        energy_frame="lab",
        products=products,
    ),
]
mcc = plasmacoll.MonteCarloCollisions(
    {"e": ELECTRON_MASS, "ion": M}, [gas], {"e": processes}, seed=7
)
electrons = plasmacoll.ParticleArrays.maxwellian(
    50_000, ELECTRON_MASS, ev_to_kelvin(20.0), seed=8
)
reactor = plasmacoll.ZeroDReactor(mcc, {"e": electrons}, volume=1e-9)

dt = 0.05 / mcc.frequency_bound("e", 1.0e7)
snapshots = {}
for steps in (0, 10, 40, 160, 640):
    while reactor.iteration < steps:
        reactor.step(dt)
    e = reactor.species["e"]
    snapshots[reactor.time] = xp.to_numpy(
        0.5 * ELECTRON_MASS * xp.sum(e.velocities**2, axis=1) / ELEMENTARY_CHARGE
    )

bins = np.linspace(0, 80, 81)
centers = 0.5 * (bins[1:] + bins[:-1])
fig, ax = plt.subplots()
for color, (time, energies) in zip(COLORS, snapshots.items(), strict=False):
    eedf, _ = np.histogram(energies, bins=bins, density=True)
    ax.semilogy(
        centers,
        np.where(eedf > 0, eedf, np.nan),
        color=color,
        label=f"t = {time * 1e9:.1f} ns",
    )
ax.axvline(11.5, color=MUTED, linewidth=1, linestyle=":")
ax.annotate(
    "excitation threshold",
    (11.5, 0.2),
    xytext=(4, -12),
    textcoords="offset points",
    color=MUTED,
)
ax.set_ylim(1e-5, 0.3)
ax.set_xlabel("electron energy (eV)")
ax.set_ylabel("EEDF (1/eV)")
ax.set_title("Inelastic collisions deplete the tail above threshold", loc="left")
ax.legend(loc="upper right");

history = reactor.state()
print("collisions per process:", history.collisions)
print(f"electrons {history.num_markers['e']}, ions {history.num_markers['ion']}")
print(f"mean electron energy {history.mean_energy_ev['e']:.2f} eV")

collisions per process: {'elastic': 514816, 'excitation': 56732, 'ionization': 22754}
electrons 72754, ions 22754
mean electron energy 6.73 eV

Path("model_gas.txt").unlink()