Skip to content
plasmacoll

Ionization and attachment

Download notebookView on GitHubRun when the site was built.

A zero-dimensional reactor has no field, push or walls. Only collisions change the species, so it is the place to check rates.

With Maxwell-molecule cross sections the ionization frequency νiz\nu_{iz} and the attachment frequency νa\nu_a are constants, and the electron density follows ne(t)=n0 e(νiz−νa)tn_e(t) = n_0\,e^{(\nu_{iz} - \nu_a)t}.

import numpy as np


import plasmacoll
from plasmacoll.constants import ATOMIC_MASS, ELECTRON_MASS, ev_to_kelvin


M = 2 * 1.008 * ATOMIC_MASS
density = 1e20
gas = plasmacoll.NeutralBackground("H2", density, 300.0, M)




def maxwell(rate):
    return plasmacoll.CrossSection.maxwell_molecule(rate, ELECTRON_MASS)




def build(k_iz, k_a, seed=0):
    processes = [
        plasmacoll.CollisionProcess(
            "ionization",
            "H2",
            maxwell(k_iz),
            name="ionization",
            energy_frame="lab",
            products={"electron": "e", "ion": "H2+"},
        ),
        plasmacoll.CollisionProcess(
            "attachment",
            "H2",
            maxwell(k_a),
            name="attachment",
            energy_frame="lab",
            products={"negative_ion": "H-"},
        ),
    ]
    mcc = plasmacoll.MonteCarloCollisions(
        {"e": ELECTRON_MASS, "H2+": M, "H-": M / 2}, [gas], {"e": processes}, seed=seed
    )
    electrons = plasmacoll.ParticleArrays.maxwellian(
        20_000, ELECTRON_MASS, ev_to_kelvin(3.0), seed=seed + 1
    )
    return plasmacoll.ZeroDReactor(mcc, {"e": electrons}, volume=1e-6)

Three species from one reaction set

Section titled “Three species from one reaction set”

Ionization at νiz=1.5×105\nu_{iz} = 1.5\times10^5 s⁻¹ competes with attachment at νa=0.5×105\nu_a = 0.5\times10^5 s⁻¹. Every ionization adds an electron and an H₂⁺, and every attachment turns an electron into an H⁻. Charge is conserved: ne+nH−=nH2++n0n_e + n_{H^-} = n_{H_2^+} + n_0.

k_iz, k_a = 1.5e-15, 0.5e-15
nu_iz, nu_a = density * k_iz, density * k_a
dt = 0.01 / (nu_iz + nu_a)
history = build(k_iz, k_a).run(dt, 1000, every=25)
t = history.time * 1e6
n0 = history.density["e"][0]


fig, ax = plt.subplots()
for color, name in zip(COLORS, ("e", "H2+", "H-"), strict=False):
    ax.plot(t, history.density[name], "o", color=color, label=name)
ax.plot(
    t,
    n0 * np.exp((nu_iz - nu_a) * history.time),
    "--",
    color=INK,
    linewidth=1.5,
    label="n₀ exp((ν_iz − ν_a) t)",
)
ax.set_yscale("log")
ax.set_xlabel("time (µs)")
ax.set_ylabel("density (m⁻³)")
ax.set_title("Electron growth with ionization and attachment", loc="left")
ax.legend()


charge = history.density["e"] + history.density["H-"] - history.density["H2+"]
print("charge balance n_e + n_H- - n_H2+ - n0:", np.max(np.abs(charge - n0)))
charge balance n_e + n_H- - n_H2+ - n0: 0.0
Notebook output

The time-step error

Section titled “The time-step error”

The null-collision method lets a marker collide at most once per step, and the new electrons do not collide in the step that created them. The electron count therefore grows by (1+p)(1 + p) per step, with p≈νΔtp \approx \nu\Delta t, instead of eνΔte^{\nu\Delta t}. Over a fixed time the relative deficit is first order in νΔt\nu\Delta t. In practice the step is chosen by the field solver and the plasma frequency, and it is usually already small enough.

final_time = 2.0 / (density * k_iz)
steps_per_run = np.array([5, 10, 20, 40, 80, 160])
errors = []
for steps in steps_per_run:
    reactor = build(k_iz, 0.0, seed=int(steps))
    history = reactor.run(final_time / steps, int(steps), every=int(steps))
    errors.append(1 - history.density["e"][-1] / history.density["e"][0] / np.exp(2.0))


# the discrete expectation: (1 + p) per step, with the operator's default bound safety 1.02
nu_dt = 2.0 / steps_per_run
p = (1 - np.exp(-1.02 * nu_dt)) / 1.02
expected = 1 - (1 + p) ** steps_per_run / np.exp(2.0)


fig, ax = plt.subplots()
ax.loglog(
    nu_dt, errors, "o", color=COLORS[0], label="simulation, 20 000 initial markers"
)
ax.loglog(nu_dt, expected, "--", color=INK, linewidth=1.5, label="1 − (1 + p)ⁿ / e²")
ax.set_xticks([0.02, 0.05, 0.1, 0.2, 0.4], labels=["0.02", "0.05", "0.1", "0.2", "0.4"])
ax.set_yticks([0.03, 0.1, 0.3], labels=["3 %", "10 %", "30 %"])
ax.minorticks_off()
ax.set_xlabel("ν_iz Δt")
ax.set_ylabel("relative deficit of n_e at ν_iz t = 2")
ax.set_title("The error is first order in the time step", loc="left")
ax.legend();
/home/runner/work/plasmacoll/plasmacoll/src/plasmacoll/mcc.py:846: RuntimeWarning: e collision probability per step 0.335 exceeds 0.1 (nu_max dt = 0.408); markers collide at most once per step, so collisions are lost. Reduce dt.
  result = self.collide(
/home/runner/work/plasmacoll/plasmacoll/src/plasmacoll/mcc.py:846: RuntimeWarning: e collision probability per step 0.185 exceeds 0.1 (nu_max dt = 0.204); markers collide at most once per step, so collisions are lost. Reduce dt.
  result = self.collide(
Notebook output