Shift to the comoving (random) velocity $\mathbf w=\mathbf v-\mathbf u$, where $f=n(m/2\pi k_BT)^{3/2}e^{-mw^2/2k_BT}$ is an isotropic Gaussian. The 0th moment is just the normalized Gaussian,

since $\int e^{-mw^2/2k_BT}\,d^3w=(2\pi k_BT/m)^{3/2}$. The 1st moment splits as $\mathbf v=\mathbf u+\mathbf w$; the $\mathbf w$ part is an odd integrand and drops, so $\frac1n\int\mathbf v f\,d^3v=\mathbf u+\langle\mathbf w\rangle=\mathbf u$. The 2nd central moment uses the per-axis variance $\langle w_i^2\rangle=k_BT/m$,

It is drift-independent because the central moment is measured about $\mathbf u$ — subtracting the bulk motion leaves only the thermal spread, which is precisely what pressure is. Hence the code returns $n=1.0\times10^{19}$, $u_x=1.0\times10^{5}$, and $p=nk_BT=13.8$ Pa.

Integrate the Vlasov equation $\partial_t f+\mathbf v\cdot\nabla_x f+\frac{\mathbf F}{m}\cdot\nabla_v f=0$ over $d^3v$. The first term is $\int\partial_t f\,d^3v=\partial_t n$. In the second, $\mathbf v$ is an independent coordinate, so it passes through the spatial gradient, $\int\mathbf v\cdot\nabla_x f\,d^3v=\nabla\cdot\!\int\mathbf v f\,d^3v=\nabla\cdot(n\mathbf u)$. The force is divergence-free in velocity ($\nabla_v\cdot(\mathbf E+\mathbf v\times\mathbf B)=0$), so $\mathbf F\cdot\nabla_v f=\nabla_v\cdot(\mathbf F f)$ is a pure velocity-space divergence and converts to a surface integral that vanishes as $f\to0$ at infinity,

What survives is the continuity equation,

This is the same local law $\partial_t(\text{density})+\nabla\cdot(\text{flux})=0$ carried by mass in ~CM-22 (density $\rho$) and probability in ~QM-04 (density $|\psi|^2$) — KEY BRIDGE B2; only the conserved density changes. For the rigid bump $\rho=e^{-(x-ut)^2}$ the finite-difference residual is $\approx3\times10^{-3}$, far below the $\approx0.43$ size of $\partial\rho/\partial t$ itself.

Multiply Vlasov by $m\mathbf v$ and integrate. The inertial terms give $\partial_t(mn\mathbf u)+\nabla\cdot(mn\langle\mathbf{vv}\rangle)$ with $\langle\mathbf{vv}\rangle=\mathbf u\mathbf u+\mathsf P/(mn)$, and the Lorentz term gives $qn(\mathbf E+\mathbf u\times\mathbf B)$. Expanding the derivatives and using continuity to cancel the $\mathbf u\,[\partial_t n+\nabla\cdot(n\mathbf u)]$ piece collapses the inertia to the convective form,

The hierarchy never closes: continuity (for $n$) needs $\mathbf u$, momentum (for $\mathbf u$) needs the 2nd moment $\mathsf P$, the pressure equation would need the 3rd-moment heat flux $\mathbf q$, and so on. Truncate polytropically, $p=C\rho^\gamma$ with $\gamma=(N+2)/N$: 1-D plasma-wave compression has $N=1\Rightarrow\gamma=3$, a 3-D monatomic plasma $N=3\Rightarrow\gamma=5/3$. This $\gamma$ fixes the sound speed; with $\gamma=5/3$ and $p=nk_BT$ from P1,

matching sound_speed(5/3, 1e19K_B1.16e6, 1e19*M_P) $\approx1.263\times10^5$ m/s.

With $\gamma=\frac53$, $p=nk_BT\approx160$ Pa and $\rho=nm_p\approx1.67\times10^{-8}$ kg/m³, the two speeds are

Their ratio is $v_A/c_s\approx55$, so $c_s\ll v_A$: magnetic tension is far stiffer than the gas here — the same verdict as $\beta\ll1$ in P6. Because $v_A\propto\rho^{-1/2}$, quadrupling the density multiplies $v_A$ by $1/\sqrt4=\tfrac12$. This matches sound_speed(...) $\approx1.26\times10^5$ m/s and alfven_speed(1.0, 1e19M_P) $\approx6.90\times10^6$ m/s (ratio $\approx55$), with alfven_speed(1.0, 41e19M_P) exactly half of alfven_speed(1.0, 1e19M_P).

For $\mathbf k\perp\mathbf B$ a compression squeezes the gas and the field together, so gas pressure and magnetic pressure restore the wave in parallel and their squared speeds add, $v_{\text{fast}}=\sqrt{c_s^2+v_A^2}$. Since both terms are non-negative,

with equality only when the other speed vanishes: $B\to0\Rightarrow v_A\to0\Rightarrow v_{\text{fast}}\to c_s$ (pure sound), and $p\to0\Rightarrow c_s\to0\Rightarrow v_{\text{fast}}\to v_A$ (pure compressional Alfvén). Numerically

only a hair above $v_A$ because $c_s^2/v_A^2\approx(1/55)^2\approx3\times10^{-4}$. This matches fast_magnetosonic_speed(1.263e5, 6.898e6) $\approx6.90\times10^6$ m/s $=\sqrt{c_s^2+v_A^2}$.

The magnetic pressure is

while the gas pressure is $p=nk_BT\approx160$ Pa. Their ratio is the plasma beta,

so the plasma is magnetically dominated — the field pressure beats the gas pressure by a factor $\sim1/\beta\approx2500$. With $p$ fixed and $P_B\propto B^2$, the beta scales as $\beta\propto1/B^2$, so doubling $B$ drops it to $\beta/4\approx1.0\times10^{-4}$. (A KrF discharge/e-beam plasma, by contrast, is essentially unmagnetized, $\beta\gg1$, so its momentum balance is set by pressure and collisions, not $\mathbf J\times\mathbf B$.) This matches magnetic_pressure(1.0) $\approx3.98\times10^5$ Pa and plasma_beta(1e19, 1.16e6, 1.0) $\approx4.0\times10^{-4}\ll1$, with $\beta\approx1.0\times10^{-4}$ after doubling $B$.

Pull a slab of electrons aside by $x$; the exposed surface charge $\sigma=nex$ makes a parallel-plate field $E=\sigma/\varepsilon_0=nex/\varepsilon_0$, so every electron feels a linear restoring force,

Screening is the static cousin: linearizing Boltzmann electrons $n=n_0e^{e\phi/k_BT}\approx n_0(1+e\phi/k_BT)$ in Poisson's equation gives $\nabla^2\phi=\phi/\lambda_D^2$ with

For $n=10^{18}\,\mathrm{m^{-3}}$ and $T=10^4\,$K this gives $\omega_p=5.64\times10^{10}\,$rad/s and $\lambda_D=6.90\times10^{-6}\,$m, equal to $v_{th}/\omega_p$ with $v_{th}=3.89\times10^5\,$m/s — matching plasma_frequency and debye_length.

Linearize the warm-electron fluid about a uniform background: continuity $\partial_t n_1+n_0\partial_x u=0$, momentum $m n_0\partial_t u=-en_0E-\partial_x p_1$, and Poisson $\varepsilon_0\partial_x E=-en_1$, closed by the 1-D adiabatic law $p_1=3k_BT\,n_1$ ($\gamma=3$ for one compressional degree of freedom). Eliminating $u,E,p_1$ for a plane wave $e^{i(kx-\omega t)}$ collapses them to a single wave equation,

using $v_{th}=\omega_p\lambda_D$. As $k\to0$ the pressure term dies and $\omega\to\omega_p$ (the cold cutoff); since $3k^2\lambda_D^2\ge0$ the wave always rides above $\omega_p$. Numerically $k\lambda_D=0.01$ gives $\omega/\omega_p=\sqrt{1+3(0.01)^2}\approx1.0001$ and $k\lambda_D=0.5$ gives $\sqrt{1+3(0.25)}=\sqrt{1.75}\approx1.323$, matching bohm_gross.

On the slow ion timescale the light electrons stay Boltzmann-distributed, $n_{e1}/n_0=e\phi/k_BT_e$, while cold ions carry the inertia. Ion momentum $m_i(-i\omega)u=-ike\phi$ and continuity $-i\omega n_{i1}+ikn_0u=0$ give $n_{i1}=n_0ek^2\phi/(m_i\omega^2)$; substituting both densities into Poisson $-\varepsilon_0k^2\phi=-e(n_{i1}-n_{e1})$ yields

Solving for $\omega$ and using $\omega_{pi}\lambda_D=\sqrt{k_BT_e/m_i}=c_s$,

For $k\lambda_D\ll1$ this is non-dispersive sound $\omega\simeq kc_s$; the phase speed $\omega/k=c_s/\sqrt{1+k^2\lambda_D^2}$ falls as $k$ grows; and for $k\lambda_D\gg1$ it saturates at $\omega\to c_s/\lambda_D=\omega_{pi}$. That ceiling is plasma_frequency(n, e, m_i)$\approx1.32\times10^9\,$rad/s, the large-$k\lambda_D$ limit that ion_acoustic approaches.

A particle at the phase velocity $v=\omega/k$ sees a nearly static field and trades energy resonantly — those just slower than the wave gain, those just faster lose — so the net rate tracks the slope of $f_0$ there (notes §5), $\gamma=\tfrac{\pi}{2}(\omega_p^3/k^2)\,\partial_v f_0|_{\omega/k}$. For a Maxwellian $f_0\propto e^{-v^2/2v_{th}^2}$ the slope $\partial_v f_0=-(v/v_{th}^2)f_0$ is negative, so $\gamma<0$ — loss. Setting $v=\omega/k$ with the Bohm–Gross $\omega^2=\omega_p^2(1+3k^2\lambda_D^2)$ puts the exponent at $-\tfrac{v^2}{2v_{th}^2}=-\tfrac{1}{2(k\lambda_D)^2}-\tfrac32$, and the prefactor collapses via $\tfrac{\pi}{2}/\sqrt{2\pi}=\sqrt{\pi/8}$ and $\omega_p^4/(k^3v_{th}^3)=\omega_p(k\lambda_D)^{-3}$ to

It is exponentially weak for $k\lambda_D\ll1$ because the resonance $v\approx v_{th}/k\lambda_D$ then sits far on the tail where almost no particles live. This returns $\gamma/\omega_p\approx-0.020,-0.096,-0.151$ at $k\lambda_D=0.3,0.4,0.5$ and $|\gamma|/\omega_p<10^{-4}$ at $k\lambda_D=0.15$, matching landau_damping_rate.

Multiply through by $(\omega-kv_0)^2(\omega+kv_0)^2=(\omega^2-a)^2$ with $a=(kv_0)^2$, using $(\omega+kv_0)^2+(\omega-kv_0)^2=2(\omega^2+a)$:

This is a quadratic in $\omega^2$ whose discriminant simplifies to $\omega_p^2(8a+\omega_p^2)$:

The lower ($-$) root is negative iff $(2a+\omega_p^2)^2<\omega_p^2(8a+\omega_p^2)$, which reduces to $4a^2<4a\omega_p^2$, i.e. $a<\omega_p^2$ — exactly $kv_0<\omega_p$. There $\omega^2<0$ makes $\omega=i\gamma$ a purely growing mode. So two_stream_growth_rate is $>0$ at $kv_0=0.5\,\omega_p$ (inside the band) and clips to $0.0$ at $kv_0=1.5\,\omega_p$ (outside).

Inside the band the growing root of P5 gives $\gamma^2=-\omega^2=\tfrac12\big[\omega_p\sqrt{8a+\omega_p^2}-(2a+\omega_p^2)\big]$, $a=(kv_0)^2$. Maximize over $a$:

so $kv_0=\sqrt{3/8}\,\omega_p\approx0.612\,\omega_p$. Substituting back,

The peak sits at $k\sim\omega_p/v_0$, i.e. wavelength $\lambda\sim v_0/\omega_p$ — the beam's "Debye-like" length (how far a beam particle drifts in one plasma period): shorter modes have $kv_0>\omega_p$ and fall outside the band, longer ones couple weakly. Scanning two_stream_growth_rate reproduces $\mathrm{Im}\,\omega/\omega_p\approx0.3536=1/\sqrt8$ at $kv_0/\omega_p\approx0.612$.

On the Landau contour the resonant denominator splits by the Plemelj rule, $1/(v-\omega/k)\to\mathrm{P.V.}\,1/(v-\omega/k)-i\pi\,\delta(v-\omega/k)$, so the dielectric function acquires an imaginary part fixed by the residue at the resonance:

Writing $\omega=\omega_r+i\gamma$ and solving $\varepsilon=0$ to first order gives $\gamma=-\operatorname{Im}\varepsilon/(\partial_\omega\!\operatorname{Re}\varepsilon)\propto \partial_v f_0|_{\omega/k}$: the growth rate inherits the sign of the slope at $v=\omega/k$. A single-hump Maxwellian falls for all $v>0$, so the resonance lands on a negative slope and the wave damps (P4); two counter-streaming beams instead carve a valley at $v=0$ with a rising slope between them, and a mode with $\omega/k$ in that gap meets $\partial_v f_0>0$ and grows (P5–P6). Same residue, opposite sign — the two-stream $f_0$ is just the Maxwellian's falling slope flipped upward where the resonant particles sit.

With $r_1=k_{\rm pump}\,g(t)\,n_{\rm Kr}$ (pump), $r_2=k_h\,n_{\rm Kr^}n_{\rm F_2}$ (harpoon), $r_3=A\,n_{\rm KrF^}$ (radiative) and $r_4=k_q\,n_{\rm KrF^*}n_{\rm F_2}$ (quench), the ledger $\dot n_i=$ production $-$ loss is

$\mathrm{F_2}$ is the limiting reagent: it starts at $2\times10^{23}\ \mathrm{m^{-3}}$ (vs $10^{24}$ for Kr), is consumed by harpooning and never refilled. KrF\* has a fast total loss rate $A+k_q n_{\rm F_2}\sim10^{8}\ \mathrm{s^{-1}}$ (an effective lifetime of tens of ns), much faster than the pump envelope, so it sits in quasi-steady state:

tracking the instantaneous formation rate $r_2\propto g(t)\,n_{\rm Kr^}$. While the Gaussian pump fills Kr\, $r_2$ climbs; once the pump passes and Kr\ drains, $r_2$ collapses and KrF\ follows it down. So KrF\ rises to a peak then decays — simulate_krf() gives peak_KrFs $\approx1.45\times10^{22}\ \mathrm{m^{-3}}$ at peak_time $\approx54$ ns, with the final KrF\ far below 10 % of the peak, and excimer_rhs is literally production minus loss.

Average the flux $\sigma v$ over the Maxwellian $f(\mathbf v)=(m/2\pi kT)^{3/2}e^{-mv^2/2kT}$; for an isotropic cross-section $d^3v=4\pi v^2\,dv$, so

Substituting $\varepsilon=\tfrac12 mv^2$ (so $v^3\,dv=2\varepsilon\,d\varepsilon/m^2$) collapses the prefactor to $\sqrt{8/\pi m}\,(kT)^{-3/2}$:

For the step $\sigma=\sigma_0$ above $E_{\rm th}$ the integral runs from $E_{\rm th}$, and $\int_{E_{\rm th}}^\infty\varepsilon\,e^{-\varepsilon/kT}\,d\varepsilon=(kT)^2(1+E_{\rm th}/kT)e^{-E_{\rm th}/kT}$, so with $\langle v\rangle=\sqrt{8kT/\pi m}$

This closed form reproduces the numeric energy integral: rate_coefficient_maxwellian(T, 3e-20, 12.0) equals rate_coefficient_step_closed_form(...) to $\sim10^{-8}$ (well under $10^{-5}$) relative, both rising with $T$.

Take the log of the closed form and differentiate with respect to $\beta\equiv1/kT$. Since $k=\sigma_0\sqrt{8/\pi m}\,(kT)^{1/2}(1+E_{\rm th}/kT)e^{-E_{\rm th}/kT}$,

The apparent activation energy is

For $E_{\rm th}\gg kT$ the middle term $\to kT$, leaving

just below threshold — so $k\propto e^{-E_{\rm th}/kT}$ is Arrhenius with $E_a\simeq E_{\rm th}$. Numerically, for $E_{\rm th}=15$ eV the slope of $\ln k$ between 12 000 and 14 000 K gives $E_a\approx14.5$ eV $=0.97\,E_{\rm th}$ (inside the 0.85–1.0 band), while arrhenius(T, A, Ea) returns its input $E_a$ exactly because $\ln\!\big(Ae^{-E_a/kT}\big)$ is exactly linear in $1/kT$.

Detailed balance equates the ionization rate $A\to A^++e$ with the recombination rate $A^++e\to A$; in equilibrium the populations obey the Saha relation

where $(2\pi m_e kT/h^2)^{3/2}$ is the electron quantum concentration (inverse thermal-de-Broglie volume) and the factor $2$ counts the two electron spin states. Put $n_e=n_i$, $n_0=n_{\rm tot}-n_i$, and let $x=n_i/n_{\rm tot}$, $R=S/n_{\rm tot}$:

As $T\to0$, $S\to0$ so $x\to0$ (neutral); as $T\to\infty$, $S\to\infty$ so $x\to1$ (fully ionized), rising monotonically between. Hence saha_ionization_fraction(T, 1e24, 14.0) climbs from $\approx0.02$ at $10^4$ K to $\approx0.99$ at $3\times10^4$ K.

The upper $B/C$ state is ionic, $\mathrm{Kr^+F^-}$, Coulomb-bound with a well a few eV deep; the lower state correlates to ground-state $\mathrm{Kr}+\mathrm{F}$ and is repulsive (no chemical bond). A vertical $B\to X$ photon therefore lands on the repulsive wall and the lower-state pair flies apart in $\sim$ps, so its population is pinned near zero. The inversion is then automatic — the instant any KrF\* exists,

with no lower-level depletion scheme needed (an effectively empty-lower-level laser). The band-centre photon energy is

So photon_energy_eV(248) $\approx4.999$ eV, and since every radiative event adds one bound–free photon, simulate_krf() accumulates a monotonically increasing photons count (the gain/Einstein-coefficient side is ~QO-03).

Add the rate equations with the right stoichiometric weights. Every Kr nucleus sits in exactly one of {Kr, Kr\, KrF\}, so (using the code's channels, with quench returning Kr + F)

$\mathrm{F_2}$ carries two F atoms, F and KrF\* one each, so

Both vanish because each reaction merely relabels atoms. What is not conserved is the raw particle count: the dissociative emission $\mathrm{KrF^*}\to\mathrm{Kr}+\mathrm{F}+h\nu$ splits one heavy particle into two, so $\sum_i n_i$ grows ($\approx+4\%$ here). Accordingly simulate_krf() keeps kr_nuclei and f_nuclei constant to $<10^{-6}$ relative across the pulse while the total particle count grows.

The energy-space Maxwellian is $F(\varepsilon)=\frac{2}{\sqrt\pi}(kT)^{-3/2}\sqrt\varepsilon\,e^{-\varepsilon/kT}$. With $\int_0^\infty\varepsilon^{s-1}e^{-\varepsilon/kT}\,d\varepsilon=\Gamma(s)(kT)^s$ and $\Gamma(\tfrac32)=\tfrac{\sqrt\pi}{2}$, $\Gamma(\tfrac52)=\tfrac{3\sqrt\pi}{4}$,

A discharge EEDF is generally not Maxwellian: the reduced field $E/N$ feeds energy into the electrons while inelastic thresholds (excitation, ionization) drain the high-energy tail, so $F(\varepsilon)$ develops a depleted/enhanced tail no single $T$ can fit — it must be solved from the electron kinetic (Boltzmann) equation (~PK-01). That is why maxwell_energy_pdf integrates to 1 with $\langle\varepsilon\rangle=\tfrac32 kT$ here, whereas phantom11 replaces the Maxwellian with LoKI Monte-Carlo tables read at $E/N=$ CHOW_EN_TD_DEFAULT $=51$ Td.