Stopping and shine-through
Attenuate a beam along its path through the plasma: where its atoms are ionised (the birth profile), how much shines through onto the wall, and how both depend on the density and the beam energy.
A beam stream decays as dI/ds = -n_e sigma_s I, with the effective stopping cross section sigma_s of every ionising channel: ionisation by electrons and ions and charge exchange, including multi-step paths through excited states. deposit traces each marker on a straight line in steps ds, adds n_e sigma_s ds to the optical depth at each step inside the last closed surface, and tallies the atoms each step ionises on flux shells.
import matplotlib.pyplot as pltimport numpy as np
import freeofcharge as focfrom freeofcharge.profiles import IonSpecies
stopping = foc.SuzukiStopping()energy_per_amu = np.geomspace(20e3, 2e6, 60) # eVfig, ax = plt.subplots(figsize=(6, 3.8), layout="constrained")for ne, style in ((1e19, "-"), (1e20, "--")): pure = np.asarray(stopping(energy_per_amu, ne, 1e4, [ne], [2], [1])) helium = np.asarray( stopping(energy_per_amu, ne, 1e4, [0.9 * ne, 0.05 * ne], [2, 4], [1, 2]) ) ax.loglog(energy_per_amu / 1e3, pure, style, color="C0", label=f"D, n_e = {ne:.0e}") ax.loglog( energy_per_amu / 1e3, helium, style, color="C1", label=f"D + 5 % He, n_e = {ne:.0e}", ) assert np.all(np.diff(pure[energy_per_amu > 100e3]) < 0) # falls with energy assert np.all(helium > pure) # impurities stop the beam fasterax.set_xlabel("beam energy [keV/amu]")ax.set_ylabel("sigma_s [m^2]")ax.legend(fontsize=8)plt.show()
low, high = ( np.asarray(stopping(500e3, ne, 1e4, [ne], [2], [1])) for ne in (1e19, 1e20))assert high > low # multi-step ionisation through excited states grows with density
SuzukiStopping is the fit of Suzuki et al. (Plasma Phys. Control. Fusion 40, 2097, 1998) for H, D and T beams in hydrogen plasmas with He, B, C, N, O and Fe impurities, from 10 keV/amu to 10 MeV/amu, as implemented in ASCOT5; the tests reproduce ASCOT5’s values. The cross section falls with the beam energy, so faster beams reach deeper, and rises with the density and with impurities. The paper fits 10 to 100 and 100 to 10000 keV/amu separately, and the two fits differ by a few per cent where they meet; SuzukiStopping blends them across 100/1.1 to 110 keV/amu so the cross section is continuous (SuzukiStopping(blend=1.0) gives the published fits with their step).
A check against Beer-Lambert
Section titled “A check against Beer-Lambert”With a constant cross section (ConstantStopping) in a plasma of constant density, the shine-through of a pencil beam is exp(-n sigma L), with L the length of its path inside the last closed surface. The path length below is measured independently, on a much finer grid than the deposition’s step:
equilibrium = foc.Equilibrium.solovev(6.2, 1.8, 1.7, 5.3, wall_psi_n=1.3)origin, direction = foc.tangential_injection(5.3, 12.0)pencil = foc.NeutralBeam.gaussian(origin, direction, 1.0e6, 1.0e6, markers=1)
s = np.arange(0.0, 30.0, 1e-4)points = origin[None, :] + s[:, None] * direction[None, :]length = ( np.asarray( equilibrium.inside(np.hypot(points[:, 0], points[:, 1]), points[:, 2]) ).sum() * 1e-4)print(f"path length in the plasma: {length:.3f} m")
n = 5e19flat = foc.Plasma( rho=np.array([0.0, 1.0]), ne=np.full(2, n), te=np.full(2, 1e4), ti=np.full(2, 1e4), ions={"D": IonSpecies(2, 1, np.full(2, n))},)depth, shine = [], []for sigma in np.geomspace(1e-22, 1e-20, 9): result = foc.deposit( pencil, equilibrium, flat, foc.ConstantStopping(sigma), step=0.005, births_per_marker=0, ) depth.append(n * sigma * length) shine.append(result.shine_through_fraction)depth, shine = np.array(depth), np.array(shine)np.testing.assert_allclose(shine, np.exp(-depth), rtol=1e-3)
x = np.linspace(0, depth.max(), 200)fig, ax = plt.subplots(figsize=(5, 3.5), layout="constrained")ax.semilogy(x, np.exp(-x), "k-", label="exp(-n sigma L)")ax.semilogy(depth, shine, "o", label="deposit")ax.set_xlabel("optical depth n sigma L")ax.set_ylabel("shine-through fraction")ax.legend()plt.show()path length in the plasma: 11.985 m
Density and energy
Section titled “Density and energy”With Suzuki’s cross sections and parabolic profiles, the shine-through falls steeply with the density: the reason a beam may only be switched on above a minimum density. A faster beam shines through more but deposits deeper. The power balance closes for every run: what is not absorbed shines through.
ions = {"D": (2, 1, 0.96), "He": (4, 2, 0.02)}rho_edges = np.linspace(0.0, 1.0, 21)densities = np.array([1e19, 2e19, 3e19, 5e19, 7e19, 1e20])energies = (0.5e6, 1.0e6)shine = {}for energy in energies: beam = foc.NeutralBeam.gaussian( origin, direction, energy, 16.5e6, width=0.3, divergence=0.005, markers=400, seed=1, ) shine[energy] = [] for ne0 in densities: plasma = foc.Plasma.parabolic(ne0, 1e4, ions=ions) result = foc.deposit( beam, equilibrium, plasma, step=0.03, births_per_marker=0, rho_edges=rho_edges, ) balance = ( result.absorbed_power + result.shine_through_power - result.injected_power ) assert abs(balance) < 1e-9 * result.injected_power assert abs(result.power.sum() / result.absorbed_power - 1) < 1e-9 shine[energy].append(result.shine_through_fraction) shine[energy] = np.array(shine[energy]) assert np.all(np.diff(shine[energy]) < 0) # denser: less shine-throughassert np.all(shine[1.0e6] > shine[0.5e6]) # faster: more shine-through
fig, ax = plt.subplots(figsize=(5.5, 3.5), layout="constrained")for energy in energies: ax.semilogy(densities, 100 * shine[energy], "o-", label=f"{energy / 1e6:.1f} MeV D")ax.set_xlabel("central electron density [m^-3]")ax.set_ylabel("shine-through [%]")ax.legend()plt.show()
Birth profile and wall hits
Section titled “Birth profile and wall hits”The birth profile is the power deposited per shell, divided by the shell volumes for a density. The fast ions are also sampled (births_per_marker), each with its birth point, velocity and pitch v_par / v; they feed the slowing down. With a wall, wall_hits holds where the shine-through leaves the vessel.
surfaces = equilibrium.flux_surfaces(rho_edges)beam = foc.NeutralBeam.gaussian( origin, direction, 1.0e6, 16.5e6, width=0.3, divergence=0.005, markers=400, seed=1)fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(10, 4), layout="constrained")peaks = []for ne0 in (2e19, 1e20): plasma = foc.Plasma.parabolic(ne0, 1e4, ions=ions) result = foc.deposit( beam, equilibrium, plasma, step=0.03, births_per_marker=2, rho_edges=rho_edges, seed=2, ) density = result.power_density(surfaces) / 1e3 peaks.append(np.average(surfaces.rho, weights=result.power)) ax0.plot( surfaces.rho, density, "o-", label=f"n_e0 = {ne0:.0e}, shine-through {100 * result.shine_through_fraction:.1f} %", ) births = np.asarray(result.births.positions) ax1.plot( births[:, 0], births[:, 1], ".", markersize=2, label=f"births, n_e0 = {ne0:.0e}" ) if ne0 == 2e19: hits = np.asarray(result.wall_hits.positions) ax1.plot( hits[:, 0], hits[:, 1], "x", color="k", markersize=4, label="shine-through on the wall", )assert peaks[1] > peaks[0] # denser: deposited further outax0.set_xlabel("rho")ax0.set_ylabel("deposited power density [kW/m^3]")ax0.legend(fontsize=8)phi = np.linspace(0, 2 * np.pi, 400)for radius in ( equilibrium.boundary.vertices[:, 0].min(), equilibrium.boundary.vertices[:, 0].max(),): ax1.plot(radius * np.cos(phi), radius * np.sin(phi), color="0.6", linewidth=0.8)ax1.set_aspect("equal")ax1.set_xlabel("x [m]")ax1.set_ylabel("y [m]")ax1.legend(fontsize=7, loc="lower left")plt.show()
Charge exchange is part of sigma_s: an ion that picks up an electron and leaves as a fast neutral counts as stopped where it was ionised, not followed further. Each marker is a straight line, so the beam is not deflected and no orbit is followed before ionisation.