Collisions

This page derives the Monte Carlo Coulomb collision operator applied with swcoll = .true.: coefficients, discrete kicks, substep bounds and random streams. On this page only, \(p = v/v_0\). Usage: Running SIMPLE.

Scope and splitting

A test alpha particle scatters on a Maxwellian background of electrons and two ion species with densities \(n_a\) and temperatures \(T_a\) from the &profiles namelist. The operator always uses the alpha mass \(m_\alpha = 4m_p\) and charge \(Z_\alpha = 2\) (src/collis_alphas.f90:107,126): n_e does not enter, and a run with n_d \(\neq 4\) stops in init_collisions with v0_coll != v0. Collisions cannot be combined with classification.

Each orbit microstep \(\Delta\tau\) is followed by a collision step of the same length at fixed position (collide in src/simple_main.f90) that changes \(p\) and \(\lambda\); update_momentum re-encodes \(J_\perp\), \(\bar v_\parallel\), \(\bar p_\varphi\) (Guiding-centre model). A failed collision step rejects the microstep as a numerical failure (Integrators).

Coefficients

For species \(a \in \{1, 2, e\}\) let \(v_{T,a} = \sqrt{2T_a/m_a}\), \(x_a = p\,v_0/v_{T,a}\), \(n_e = Z_1n_1 + Z_2n_2\), \(Z_e = 1\), \(E_\alpha = m_\alpha v_0^2/2 = 3.5\,\text{MeV}/\texttt{facE\_al}\), and

\[\nu_a = \frac{2\pi Z_\alpha^2 Z_a^2 e^4 n_a \ln\Lambda_a}{m_\alpha^2 v_0^3\, v_{T,a}} \quad [\text{cm}^{-1}].\]

With the Chandrasekhar function \(G(x) = [\operatorname{erf}x - x\operatorname{erf}'x]/(2x^2)\), onseff evaluates (series for \(x < 0.01\), asymptotes for \(x > 6\))

\[d_p = \frac{2G(x)}{x}, \qquad d_\lambda = \frac{\operatorname{erf}x - G(x)}{x}, \qquad d_F = \frac1x\frac{d(x^2 d_p)}{dx} = \frac{4e^{-x^2}}{\sqrt\pi} - d_p ,\]

and coleff sums the species:

\[D_p = \sum_a \nu_a d_p(x_a), \qquad D_\lambda = \frac{1}{p^2}\sum_a \nu_a d_\lambda(x_a), \qquad F_p = \sum_a \nu_a\Bigl[\frac{d_F(x_a)}{p} - 2p\frac{E_\alpha}{T_a}d_p(x_a)\Bigr] = \frac{1}{p^2}\frac{\partial(p^2 D_p)}{\partial p} - 2p\sum_a\frac{E_\alpha}{T_a}\nu_a d_p(x_a).\]

\(F_p\) and \(D_p\) are the Itô drift and diffusion of \(p\), and \(D_\lambda\) the pitch diffusion, of the Fokker–Planck equation for \(f(p,\lambda)\) per unit velocity volume

\[\frac{\partial f}{\partial\tau} = \sum_a \frac{1}{p^2}\frac{\partial}{\partial p} \Bigl[p^2\nu_a d_p(x_a)\Bigl(\frac{\partial f}{\partial p} + 2p\frac{E_\alpha}{T_a}f\Bigr)\Bigr] + D_\lambda\frac{\partial}{\partial\lambda}\Bigl[(1-\lambda^2)\frac{\partial f}{\partial\lambda}\Bigr],\]

whose speed flux vanishes for each species on the Maxwellian \(f\propto\exp(-m_\alpha v^2/2T_a)\).

The Coulomb logarithms (compute_collision_coeffs; \(n\) in cm^-3, \(T\), \(E_\alpha\) in eV, \(A_i\) = am1, am2) are the NRL electron–ion form and the NRL mixed ion–ion form with \(m_\alpha = 4m_p\), \(T_\alpha \to E_\alpha\) and only the background-ion density under the root:

\[\ln\Lambda_e = 24 - \ln\frac{\sqrt{n_e}}{T_e}, \qquad \ln\Lambda_i = 23 - \ln\Bigl[\frac{Z_\alpha Z_i(4 + A_i)}{4T_i + A_iE_\alpha} \sqrt{\frac{n_iZ_i^2}{T_i}}\Bigr].\]

init_collision_profiles tabulates \(\nu_a\), \(v_0/v_{T,a}\) and \(E_\alpha/T_a\) on 101 equidistant points in \([0,1]\), with floors \(T_a \ge 1\) eV and \(n_a \ge 10^6\) m^-3 (both ion densities zero is vacuum). Each step interpolates linearly at the clamped normalised toroidal flux \(s\) of the marker, which profile_s (src/profile_radius.f90) computes from the tracing state: \(z_1\) itself for integ_coords -1, 0, 1 and 2; for 3 and 4 the radius of the source field after integ_to_ref; for 5 the reference radius; for 6 the enclosed toroidal flux of the SPECTRE chart \(\rho_g\), normalised to tflux(Mvol). A rho_tor chartmap radius gives \(s = \rho^2\), and vmec_extended/vmec_to_wall charts give \(s = (\rho/\rho_\mathrm{lcfs})^2\). Radii without toroidal-flux meaning (geometric or poloidal chartmaps) are accepted only for radially uniform &profiles.

Discrete step

stost (mode 1) advances \((p,\lambda)\) over \(\Delta\tau\) in substeps \(\delta\), with coefficients at the current \(p\):

\[\lambda' = \lambda\cos\alpha + \sqrt{1-\lambda^2}\,\sin\alpha\cos\psi, \quad \alpha = \sqrt{4D_\lambda\delta},\ \psi\sim\mathcal U[0,2\pi); \qquad p' = p + F_p\delta + \sqrt{2D_p\delta}\,\xi,\quad \xi\sim\mathcal U[-\sqrt3,\sqrt3].\]

The pitch rotation has \(\langle\Delta\lambda\rangle = -2\lambda D_\lambda\delta\) and \(\langle\Delta\lambda^2\rangle = 2(1-\lambda^2)D_\lambda\delta\) to leading order, the generator of the Lorentz operator, and keeps \(|\lambda| \le 1\) exactly; \(\xi\) has zero mean and unit variance. Before the random draws,

\[\delta = \min\Bigl(\text{remaining},\ \frac{f^2}{D_\lambda},\ \frac{f\,p}{|F_p|},\ \frac{f^2p^2}{6D_p}\Bigr), \qquad f = 0.1,\]

which bounds the rotation angle by \(2f\) and the drift and the largest noise amplitude \(\sqrt{6D_p\delta}\) by \(fp\). stost accepts \(0 < f \le 0.25\) (max_substep_fraction); simple.x uses 0.1. A non-finite or negative coefficient, \(p \le 0\), a substep that no longer advances time, or more than 100,000 substeps restores the initial state and returns an error. In vacuum neither the state nor the generator changes.

Monoenergetic Lorentz mode

collision_model = 'lorentz_nustar' (with swcoll = .true.) replaces the Coulomb operator by pitch-angle scattering at fixed \(p\). The collision frequency is \(\nu = \nu_* |\iota| v_0 / R_0\), with nu_star_standard \(= \nu_*\), lorentz_iota \(= \iota\) and lorentz_major_radius_cm \(= R_0\); all three must be set (\(\nu_* \ge 0\), \(R_0 > 0\), \(\iota \ne 0\)). With path length \(\tau = v_0 t\), one call to stost_lorentz applies the Euler kick \(\Delta\lambda = \sqrt{2 D (1-\lambda^2)\Delta\tau}\,\xi - 2\lambda D \Delta\tau\), \(D = \nu/(2 v_0)\), \(\xi = \pm 1\), reflected at \(|\lambda| = 1\). It has no substep control: keep \(\nu \Delta t\) small. p_abs is unchanged. The particle energy can be set directly with particle_energy_eV (default -1: use \(3.5\,\mathrm{MeV}/\) facE_al); a nondefault facE_al that disagrees with an explicit energy is rejected. Mode and energy are written to orbits.nc and results.nc as global attributes.

Random streams

With deterministic = .true., reset_seed_if_deterministic (src/params.f90) reseeds the generator before each marker from ran_seed and \(k\) = particle_index_offset + source row:

\[s_0 = (\texttt{ran\_seed} + 104729\,k) \bmod (2^{31}-2) + 1, \qquad s_i = (48271\,s_{i-1} + 7919\,i) \bmod (2^{31}-2) + 1 .\]

A marker thus draws the same kicks for any thread count or contiguous shard; distinct seeds do not prove statistical independence.

References

• J. D. Huba, NRL Plasma Formulary, NRL/PU/6790–19-640, Naval Research Laboratory (2019)