# SPDX-License-Identifier: AGPL-3.0-or-later
# Commercial license available
# © Concepts 1996–2026 Miroslav Šotek. All rights reserved.
# © Code 2020–2026 Miroslav Šotek. All rights reserved.
# ORCID: 0009-0009-3560-0851
# Contact: www.anulum.li | protoscience@anulum.li
# SCPN Fusion Core — Relativistic Runaway Electron Dynamics
"""Runaway electron utility contracts for DREAM-style reduced-order workflows."""
from __future__ import annotations
from dataclasses import dataclass
from typing import Callable
import numpy as np
from numpy.typing import NDArray
FloatArray = NDArray[np.float64]
E_CHARGE = 1.602176634e-19 # C
M_E = 9.1093837e-31 # kg
C_LIGHT = 2.99792458e8 # m/s
EPS_0 = 8.8541878128e-12 # F/m
[docs]
@dataclass
class RunawayParams:
"""Scalar plasma parameters used by reduced-order runaway balances."""
ne_20: float # electron density [10^20 m^-3]
Te_keV: float # electron temperature [keV]
E_par: float # parallel electric field [V/m]
Z_eff: float # effective charge
B0: float # magnetic field [T]
R0: float # major radius [m]
a: float = 2.0 # minor radius [m]
[docs]
@dataclass(frozen=True)
class DreamFluidBalance:
"""DREAM-style fluid runaway density balance for a single plasma state."""
dreicer_source: float
avalanche_source: float
loss_source: float
total_source: float
runaway_fraction: float
growth_time_s: float
[docs]
def dreicer_field(ne_20: float, Te_keV: float, coulomb_log: float = 15.0) -> float:
"""Dreicer field E_D = n_e e^3 lnΛ / (4π ε₀² T_e).
Connor & Hastie, Nucl. Fusion 15, 415 (1975).
"""
n_e = ne_20 * 1e20
Te_J = Te_keV * 1e3 * E_CHARGE
return float(n_e * E_CHARGE**3 * coulomb_log / (4.0 * np.pi * EPS_0**2 * Te_J))
[docs]
def critical_field(ne_20: float, coulomb_log: float = 15.0) -> float:
"""Critical field E_c = n_e e^3 lnΛ / (4π ε₀² m_e c²).
Rosenbluth & Putvinski, Nucl. Fusion 37, 1355 (1997).
E_D / E_c = m_e c² / T_e ≈ 51 at T_e = 10 keV.
"""
n_e = ne_20 * 1e20
return float(n_e * E_CHARGE**3 * coulomb_log / (4.0 * np.pi * EPS_0**2 * M_E * C_LIGHT**2))
[docs]
def coulomb_logarithm_thermal(ne_20: float, Te_keV: float) -> float:
"""Thermal electron-electron Coulomb logarithm.
ln Lambda_T = 14.9 - 0.5 ln(n_e / 10^20 m^-3) + ln(T_e / 1 keV) — the
standard thermal form, identical to DREAM's
``CoulombLogarithm::evaluateLnLambdaT`` (source-verified at DREAM commit
a08edc0d, src/Equations/CoulombLogarithm.cpp).
"""
if not np.isfinite(ne_20) or ne_20 <= 0.0:
raise ValueError("ne_20 must be finite and > 0.")
if not np.isfinite(Te_keV) or Te_keV <= 0.0:
raise ValueError("Te_keV must be finite and > 0.")
return float(14.9 - 0.5 * np.log(ne_20) + np.log(Te_keV))
[docs]
def coulomb_logarithm_relativistic(ne_20: float, Te_keV: float) -> float:
"""Relativistic Coulomb logarithm.
ln Lambda_c = 14.6 + 0.5 ln(T_e[eV] / (n_e / 10^20 m^-3)) — Hesslow et
al., Nucl. Fusion 59, 084004 (2019) (arXiv:1904.00602), identical to
DREAM's ``CoulombLogarithm::evaluateLnLambdaC``.
"""
if not np.isfinite(ne_20) or ne_20 <= 0.0:
raise ValueError("ne_20 must be finite and > 0.")
if not np.isfinite(Te_keV) or Te_keV <= 0.0:
raise ValueError("Te_keV must be finite and > 0.")
return float(14.6 + 0.5 * np.log(Te_keV * 1e3 / ne_20))
# Classical electron radius r_0 = e^2 / (4 pi eps_0 m_e c^2) [m]
R_0_CLASSICAL = E_CHARGE**2 / (4.0 * np.pi * EPS_0 * M_E * C_LIGHT**2)
[docs]
def dreicer_generation_rate(params: RunawayParams, coulomb_log: float | None = None) -> float:
"""Primary (Dreicer) runaway generation rate [m^-3 s^-1].
Corrected Connor & Hastie (1975) rate in the exact form implemented by
DREAM's ``ConnorHastie::RunawayRate`` with corrections enabled
(``DREICER_RATE_CONNOR_HASTIE``; source-verified at DREAM commit
a08edc0d, src/Equations/ConnorHastie.cpp):
gamma_D = (n_e / tau_EE) (E/E_D)^alpha
exp(-lambda/(4 E/E_D) - sqrt(eta (1+Z) / (E/E_D)))
with ``alpha = -3 (1+Z) h / 16``, the h/eta/lambda correction factors
functions of E/E_c, and ``tau_EE = beta_th^3 / (4 pi r_0^2 c n_e
lnLambda_T)`` the thermal collision time. When *coulomb_log* is None
(default), the thermal logarithm lnLambda_T is computed from the plasma
state and the critical-field threshold uses the relativistic
lnLambda_c; passing a value forces both (legacy fixed-lnLambda mode).
Validated against a really-executed DREAM fluid reference
(``validation/reference_data/dream/``): ratio ours/DREAM = 1.007 at the
committed same-case state. The previous implementation used C_D = 0.35,
a malformed exponent, and fixed lnLambda = 15 — it was 5.4x LOW
(BACKLOG 3 finding, closed by this form).
"""
if params.E_par <= 0.0 or params.Te_keV <= 0.0 or params.ne_20 <= 0.0:
return 0.0
ln_lambda_t = (
coulomb_logarithm_thermal(params.ne_20, params.Te_keV)
if coulomb_log is None
else coulomb_log
)
ln_lambda_c = (
coulomb_logarithm_relativistic(params.ne_20, params.Te_keV)
if coulomb_log is None
else coulomb_log
)
E_c = critical_field(params.ne_20, ln_lambda_c)
if params.E_par <= E_c:
return 0.0
n_e = params.ne_20 * 1e20
E_D = dreicer_field(params.ne_20, params.Te_keV, ln_lambda_t)
E_Ec = params.E_par / E_c
E_ED = params.E_par / E_D
if E_ED < 1e-4:
return 0.0
Z = params.Z_eff
h = (E_Ec + 2.0 * (E_Ec - 2.0) * np.sqrt(E_Ec / (E_Ec - 1.0)) - (Z - 7.0) / (Z + 1.0)) / (
3.0 * (E_Ec - 1.0)
)
eta_f = 0.5 * np.pi - np.arcsin(1.0 - 2.0 / E_Ec)
eta = E_Ec**2 / (4.0 * (E_Ec - 1.0)) * eta_f**2
lmbd = 8.0 * E_Ec**2 * (1.0 - 1.0 / (2.0 * E_Ec) - np.sqrt(1.0 - 1.0 / E_Ec))
alpha = -3.0 / 16.0 * (1.0 + Z) * h
beta_th = np.sqrt(2.0 * params.Te_keV * 1e3 * E_CHARGE / (M_E * C_LIGHT**2))
tau_ee = beta_th**3 / (4.0 * np.pi * R_0_CLASSICAL**2 * C_LIGHT * n_e * ln_lambda_t)
exponent = -lmbd / (4.0 * E_ED) - np.sqrt(eta * (1.0 + Z) / E_ED)
if exponent < -500: # pragma: no cover - defensive underflow guard
return 0.0
rate = n_e / tau_ee * E_ED**alpha * np.exp(exponent)
if not np.isfinite(rate): # pragma: no cover - defensive near-threshold guard
return 0.0
return float(max(0.0, rate))
[docs]
def avalanche_growth_rate(
params: RunawayParams, n_RE: float, coulomb_log: float | None = None
) -> float:
"""Avalanche multiplication source [m^-3 s^-1] (``Gamma_ava n_RE``).
Rosenbluth-Putvinski growth rate in the compact form quoted by Hesslow
et al., Nucl. Fusion 59, 084004 (2019) (arXiv:1904.00602, page 2:
``Gamma_0 = e / (lnLambda m_e c sqrt(5 + Z_eff))``), with the
critical-field threshold:
Gamma_ava = e (E_par - E_c) / (m_e c lnLambda_c sqrt(5 + Z_eff))
When *coulomb_log* is None (default) the relativistic Coulomb logarithm
is computed from the plasma state; passing a value forces it (legacy
fixed-lnLambda mode). Compared against a really-executed DREAM fluid
reference (``validation/reference_data/dream/``): ratio ours/DREAM =
0.75 at the committed same-case state — the residual is the compact RP
form versus DREAM's matched-formula effective critical momentum
(Hesslow et al. 2019), not a unit or factor error. The previous
implementation omitted the ``1/(lnLambda sqrt(5+Z_eff))`` factor
entirely and was 32x HIGH (BACKLOG 3 finding, closed by this form).
"""
if n_RE <= 0.0 or params.E_par <= 0.0:
return 0.0
if params.Te_keV <= 0.0 or params.ne_20 <= 0.0:
return 0.0
ln_lambda_c = (
coulomb_logarithm_relativistic(params.ne_20, params.Te_keV)
if coulomb_log is None
else coulomb_log
)
E_c = critical_field(params.ne_20, ln_lambda_c)
if params.E_par <= E_c:
return 0.0
gamma_ava = (
E_CHARGE
* (params.E_par - E_c)
/ (M_E * C_LIGHT * ln_lambda_c * np.sqrt(5.0 + params.Z_eff))
)
return float(max(0.0, n_RE * gamma_ava))
[docs]
def hot_tail_seed(
Te_pre_keV: float,
Te_post_keV: float,
ne_20: float,
quench_time_ms: float,
*,
vc_vte_ref: float = 4.0,
quench_exponent: float = 0.2,
) -> float:
"""Seed RE density from thermal quench hot-tail mechanism [m^-3].
Smith, H.M. et al., Phys. Plasmas 15, 072502 (2008).
Faster TQ → lower v_c/v_te → more seed electrons.
"""
te_pre = float(Te_pre_keV)
te_post = float(Te_post_keV)
ne = float(ne_20)
tau_q = float(quench_time_ms)
vc_ref = float(vc_vte_ref)
q_exp = float(quench_exponent)
if not np.isfinite(te_pre) or te_pre <= 0.0:
raise ValueError("Te_pre_keV must be finite and > 0.")
if not np.isfinite(te_post) or te_post <= 0.0:
raise ValueError("Te_post_keV must be finite and > 0.")
if not np.isfinite(ne) or ne <= 0.0:
raise ValueError("ne_20 must be finite and > 0.")
if not np.isfinite(tau_q) or tau_q <= 0.0:
raise ValueError("quench_time_ms must be finite and > 0.")
if not np.isfinite(vc_ref) or vc_ref <= 0.0:
raise ValueError("vc_vte_ref must be finite and > 0.")
if not np.isfinite(q_exp) or q_exp <= 0.0:
raise ValueError("quench_exponent must be finite and > 0.")
if te_post >= te_pre:
return 0.0
ratio = te_pre / te_post
# v_c/v_te parametric fit to Smith 2008 Fig. 3
v_c_v_te = vc_ref * (tau_q / 1.0) ** q_exp
n_e = ne * 1e20
exponent = -(v_c_v_te**2)
if (
exponent < -500
): # pragma: no cover - defensive (np.exp underflow guard for extreme v_c/v_te)
return 0.0
n_seed = n_e * v_c_v_te**3 * np.exp(exponent) * ratio**1.5
return float(max(0.0, n_seed))
[docs]
def dream_fluid_density_balance(
params: RunawayParams,
n_RE: float,
*,
loss_time_s: float = np.inf,
max_runaway_fraction: float = 1.0,
coulomb_log: float | None = None,
) -> DreamFluidBalance:
"""Evaluate the scalar density balance used by DREAM-style fluid runs.
The contract is ``dn_RE/dt = S_Dreicer + gamma_avalanche n_RE - n_RE/tau_loss``.
It is a fluid benchmark contract, not a kinetic DREAM distribution solver.
"""
n_re = float(n_RE)
tau_loss = float(loss_time_s)
max_fraction = float(max_runaway_fraction)
if not np.isfinite(n_re) or n_re < 0.0:
raise ValueError("n_RE must be finite and >= 0.")
if not np.isfinite(max_fraction) or max_fraction <= 0.0 or max_fraction > 1.0:
raise ValueError("max_runaway_fraction must be finite and in (0, 1].")
if np.isnan(tau_loss) or tau_loss <= 0.0:
raise ValueError("loss_time_s must be positive or infinite.")
n_e = params.ne_20 * 1e20
if not np.isfinite(n_e) or n_e <= 0.0:
raise ValueError("electron density must be finite and > 0.")
if n_re > n_e * max_fraction:
raise ValueError("n_RE exceeds the configured runaway density cap.")
dreicer_source = dreicer_generation_rate(params, coulomb_log)
avalanche_source = avalanche_growth_rate(params, n_re, coulomb_log)
loss_source = 0.0 if np.isinf(tau_loss) else n_re / tau_loss
total_source = dreicer_source + avalanche_source - loss_source
growth_time_s = n_re / total_source if total_source > 0.0 and n_re > 0.0 else np.inf
return DreamFluidBalance(
dreicer_source=float(dreicer_source),
avalanche_source=float(avalanche_source),
loss_source=float(loss_source),
total_source=float(total_source),
runaway_fraction=float(n_re / n_e),
growth_time_s=float(growth_time_s),
)
[docs]
class RunawayEvolution:
"""Time-domain reduced-order solver for runaway electron density evolution."""
def __init__(self, params: RunawayParams):
self.params = params
[docs]
def balance(
self,
n_RE: float,
E_par: float,
*,
loss_time_s: float = np.inf,
max_runaway_fraction: float = 1.0,
) -> DreamFluidBalance:
"""Return the instantaneous DREAM-style fluid density balance."""
self.params.E_par = float(E_par)
return dream_fluid_density_balance(
self.params,
n_RE,
loss_time_s=loss_time_s,
max_runaway_fraction=max_runaway_fraction,
)
[docs]
def step(
self,
dt: float,
n_RE: float,
E_par: float,
*,
loss_time_s: float = np.inf,
max_runaway_fraction: float = 1.0,
) -> float:
"""Advance one explicit time step of the reduced-order runaway balance."""
dt_s = float(dt)
n_re = float(n_RE)
e_par = float(E_par)
if not np.isfinite(dt_s) or dt_s <= 0.0:
raise ValueError("dt must be finite and > 0.")
if not np.isfinite(n_re) or n_re < 0.0:
raise ValueError("n_RE must be finite and >= 0.")
if not np.isfinite(e_par):
raise ValueError("E_par must be finite.")
balance = self.balance(
n_re,
e_par,
loss_time_s=loss_time_s,
max_runaway_fraction=max_runaway_fraction,
)
n_cap = self.params.ne_20 * 1e20 * float(max_runaway_fraction)
return float(min(n_cap, max(n_re + balance.total_source * dt_s, 0.0)))
[docs]
def evolve(
self,
n_RE_0: float,
E_par_profile: Callable[[float], float],
t_span: tuple[float, float],
dt: float,
) -> tuple[FloatArray, FloatArray]:
"""Integrate runaway density over a profile-driven E∥(t) history."""
if not isinstance(t_span, tuple) or len(t_span) != 2:
raise ValueError("t_span must be a (t_start, t_end) tuple.")
t_start, t_end = float(t_span[0]), float(t_span[1])
dt_s = float(dt)
n0 = float(n_RE_0)
if not np.isfinite(dt_s) or dt_s <= 0.0:
raise ValueError("dt must be finite and > 0.")
if not np.isfinite(t_start) or not np.isfinite(t_end) or t_end <= t_start:
raise ValueError("t_span must be finite with t_end > t_start.")
if not np.isfinite(n0) or n0 < 0.0:
raise ValueError("n_RE_0 must be finite and >= 0.")
n_steps = int(np.ceil((t_end - t_start) / dt_s))
if (
n_steps < 1
): # pragma: no cover - unreachable (t_end > t_start and dt > 0 already validated)
raise ValueError("integration requires at least one time step.")
t_arr = np.linspace(t_start, t_end, n_steps + 1)
n_RE_arr = np.zeros(n_steps + 1)
n_RE_arr[0] = n0
for i in range(n_steps):
E_par = E_par_profile(t_arr[i])
if not np.isfinite(E_par):
raise ValueError("E_par_profile must return finite values.")
n_RE_arr[i + 1] = self.step(dt_s, n_RE_arr[i], float(E_par))
return t_arr, n_RE_arr
[docs]
def current_fraction(self, n_RE: float, I_p_MA: float) -> float:
"""Fraction of total plasma current carried by REs (v ≈ c)."""
n_re = float(n_RE)
i_p_ma = float(I_p_MA)
if not np.isfinite(n_re) or n_re < 0.0:
raise ValueError("n_RE must be finite and >= 0.")
if not np.isfinite(i_p_ma):
raise ValueError("I_p_MA must be finite.")
if i_p_ma <= 0.0:
return 0.0
j_RE = E_CHARGE * n_re * C_LIGHT
I_RE = j_RE * np.pi * self.params.a**2
return float(min(1.0, I_RE / (i_p_ma * 1e6)))
[docs]
class RunawayMitigationAssessment:
"""Collection of reduced-order mitigation formulas and guard checks."""
[docs]
@staticmethod
def required_density_for_suppression(
E_par: float, Z_eff: float, coulomb_log: float = 15.0
) -> float:
"""Density [10^20 m^-3] needed to make E_c > E_par."""
e_par = float(E_par)
z_eff = float(Z_eff)
ln_lambda = float(coulomb_log)
if not np.isfinite(e_par) or e_par <= 0.0:
return 0.0
if not np.isfinite(z_eff) or z_eff <= 0.0:
raise ValueError("Z_eff must be finite and > 0.")
if not np.isfinite(ln_lambda) or ln_lambda <= 0.0:
raise ValueError("coulomb_log must be finite and > 0.")
# E_c = n_e * e^3 lnL / (4 pi eps0^2 m_e c^2)
# Effective collisional drag scales with impurity content.
collision_multiplier = 0.5 * (1.0 + z_eff)
coeff = collision_multiplier * E_CHARGE**3 * ln_lambda
# n_e = E_par * (4 pi eps0^2 m_e c^2) / (e^3 lnL)
n_e_m3 = e_par * (4.0 * np.pi * EPS_0**2 * M_E * C_LIGHT**2) / coeff
return float(n_e_m3 / 1e20)
[docs]
@staticmethod
def maximum_re_energy(
B0: float,
R0: float,
*,
ne_20: float = 1.0,
Z_eff: float = 1.5,
pitch_angle_rad: float = 0.35,
) -> float:
"""Maximum RE energy [MeV] from orbit and synchrotron limits."""
b = float(B0)
r = float(R0)
ne = float(ne_20)
zeff = float(Z_eff)
pitch = float(pitch_angle_rad)
if not np.isfinite(b) or b <= 0.0:
raise ValueError("B0 must be finite and > 0.")
if not np.isfinite(r) or r <= 0.0:
raise ValueError("R0 must be finite and > 0.")
if not np.isfinite(ne) or ne <= 0.0:
raise ValueError("ne_20 must be finite and > 0.")
if not np.isfinite(zeff) or zeff <= 0.0:
raise ValueError("Z_eff must be finite and > 0.")
if not np.isfinite(pitch) or pitch <= 0.0 or pitch >= (0.5 * np.pi):
raise ValueError("pitch_angle_rad must be finite and in (0, pi/2).")
# Drift-orbit limit: E_orbit[eV] ~ B[T] R[m] c[m/s].
orbit_limit_mev = (b * r * C_LIGHT) * 1.0e-6
# Synchrotron-loss constrained practical limit around ITER conditions.
sin_pitch = np.sin(pitch)
synch_limit_mev = 28.0
synch_limit_mev *= (b / 5.3) ** (1.0 / 3.0)
synch_limit_mev *= (r / 6.2) ** (2.0 / 3.0)
synch_limit_mev *= (ne / 1.0) ** (-1.0 / 6.0)
synch_limit_mev *= (zeff / 1.5) ** (-1.0 / 6.0)
synch_limit_mev *= (max(sin_pitch, 1e-3) / np.sin(0.35)) ** (-1.0 / 3.0)
return float(max(0.1, min(orbit_limit_mev, synch_limit_mev)))
[docs]
@staticmethod
def wall_heat_load(
n_RE: float,
E_max_MeV: float,
A_wet: float,
volume: float = 800.0,
*,
mean_energy_MeV: float | None = None,
) -> float:
"""Deposited runaway-electron beam energy density [MJ/m^2]."""
n_re = float(n_RE)
e_max = float(E_max_MeV)
area = float(A_wet)
vol = float(volume)
if not np.isfinite(n_re) or n_re < 0.0:
raise ValueError("n_RE must be finite and >= 0.")
if not np.isfinite(e_max) or e_max <= 0.0:
raise ValueError("E_max_MeV must be finite and > 0.")
if not np.isfinite(vol) or vol <= 0.0:
raise ValueError("volume must be finite and > 0.")
if mean_energy_MeV is None:
e_mean = 0.5 * e_max
else:
e_mean = float(mean_energy_MeV)
if not np.isfinite(e_mean) or e_mean <= 0.0 or e_mean > e_max:
raise ValueError("mean_energy_MeV must be finite and within (0, E_max_MeV]")
E_avg_J = e_mean * 1e6 * E_CHARGE
W_total_J = n_re * vol * E_avg_J
if area <= 0.0:
return float("inf")
load_MJ_m2 = (W_total_J / 1e6) / area
return float(load_MJ_m2)