Knock-on collisions

max_knockon_momentum(p1)

Largest knock-on momentum from an incident electron of momentum p1, in units of mₑc: the momentum at half the incident kinetic energy, γ_max = (γ₁ + 1)/2. Above it, the outgoing electron with momentum p is the primary, not the knock-on.

Ported from DREAM KnockOnUtilities.hpp (Kinematics::MaximumKnockOnMomentum).

moller_cross_section(p, p1)

Møller differential cross section Σ = c dσ/dγ [m³/s] for producing a knock-on electron of momentum p from an incident electron of momentum p1 on a stationary target, both in units of mₑc. The rate coefficient per target electron is (p1/γ1) Σ dγ.

Ported from DREAM KnockOnUtilities.hpp (Kinematics::EvaluateMollerDifferentialCrossSection).

moller_flux(p, p1)

Antiderivative S [m³/s] of the Møller cross section in the knock-on momentum p (in units of mₑc), ∂S/∂p = (p/γ) Σ, i.e. ∂S/∂γ = Σ. Integrating the knock-on source over a momentum cell reduces to differences of S.

Ported from DREAM KnockOnUtilities.hpp (Kinematics::EvaluateMollerFlux).

moller_total_rate(p1, p_cutoff)

Rate coefficient [m³/s] for an incident electron of momentum p1 to produce a knock-on electron with momentum above p_cutoff, both in units of mₑc, (p1/γ1) ∫ Σ dγ = (p1/γ1) [S(p_max) - S(p_cutoff)] with p_max = max_knockon_momentum(p1); zero when p_cutoff ≥ p_max. Multiplied by the target density it gives the production rate per incident electron. Corresponds to DREAM's EvaluateMollerFluxIntegratedOverKnockonGrid on an unbounded knock-on grid.

momentum_at_xi_star(xi_s, p1)

Knock-on momentum p (in units of mₑc) with xi_star(p, p1) = xi_s, the inverse of xi_star in its first argument (ξ* grows monotonically with p): from $\xi_*^2 = (\gamma_1 + 1)(\gamma - 1) / ((\gamma_1 - 1)(\gamma + 1))$,

\[\gamma - 1 = \frac{2 \xi_*^2 (\gamma_1 - 1)}{\gamma_1 + 1 - \xi_*^2 (\gamma_1 - 1)}, \qquad p = \sqrt{(\gamma - 1)(\gamma + 1)}.\]

Used to split the receiving momentum interval of a pair at the ξ* values of the pitch table (KnockOnKernels).

xi_star(p, p1)

Cosine ξ* of the angle between the incident electron (momentum p1, colliding with a stationary electron) and the knock-on electron (momentum p), both in units of mₑc, ξ* = √[(γ₁ + 1)(γ - 1) / ((γ₁ - 1)(γ + 1))] (notation of Embréus et al. 2018).

Ported from DREAM KnockOnUtilities.hpp (Kinematics::EvaluateXiStar).

MollerEnergyKernel{T}

Momentum part of the knock-on operator between a primary momentum grid (incident electrons, cells k) and a knock-on momentum grid (cells i):

\[\mathtt{rates}[i, k] = \beta_1(p_{1,k}) \, \frac{S(p_U) - S(p_L)}{\Delta p_i},\]

the Møller production rate coefficient [m³/s] averaged over knock-on cell i (S_ik of the DREAM notes), with [p_L, p_U] the part of the cell inside [p_cutoff, max_knockon_momentum(p1_k)] and the primary evaluated at its cell centre p1_k. S is moller_flux.

• total_rates[k] = Σ_i Δp_i rates[i, k]: production rate coefficient onto the knock-on grid.
• rates_above_grid[k]: production above the top face of the knock-on grid (knock-ons that leave the grid, e.g. into a fluid runaway population).

Ported from DREAM MollerEnergyKernel.cpp and KnockOnUtilities.cpp (EvaluateMollerFluxMatrixElementOnGrid, EvaluateMollerFluxIntegratedOverKnockonGrid).

MollerEnergyKernel(knockon, primary, p_cutoff)

Build the kernel for knock-on momentum grid knockon and primary momentum grid primary (both MomentumGrid; only the p grids are used) with cutoff momentum p_cutoff > 0 in units of mₑc.

MollerPitchKernel{T}

Knock-on pitch kernel δ[j, l](ξ*) between primary pitch cells l and knock-on pitch cells j, tabulated on a uniform grid of ξ* values and interpolated linearly to the ξ* of each momentum pair (i, k). Every tabulated column satisfies Σ_j Δξ_j δ[j, l] = 1, and so does every interpolated one.

• xi_star_values: tabulated ξ*, one per table slice.
• table[j, l, m]: δ of slice m.
• lower_slice[i, k], upper_weight[i, k]: the single-ξ* interpolation of pair (i, k), which uses slices lower_slice and lower_slice + 1 with weights 1 - upper_weight and upper_weight. With xi_star_mode = :interval the operator instead spreads the pair over every slice its ξ* range covers (KnockOnKernels); these two fields then describe the midpoint interpolation, which serves as a reference.
• xi_star_mode: how ξ* of a pair enters (:interval, :midpoint or :centre).
• no_contributing_pairs: no (i, k) pair produces knock-ons (cutoff above the grid).

In cylindrical geometry the kernel does not depend on radius; one table serves all radii (DREAM stores one per radius).

MollerPitchKernel(knockon, primary, p_cutoff; table_size=80, dream_compatible=false,
                  xi_star_mode=(dream_compatible ? :centre : :interval),
                  primary_pitch=(dream_compatible ? :centre : :split),
                  primary_pitch_points=4)

• The part of knock-on cell i that receives knock-ons from primary cell k is [max(p_{i-1/2}, p_cutoff), min(p_{i+1/2}, p_max(p1_k))], with the momenta in units of mₑc, the interval the energy kernel integrates over. ξ* varies over it, fastest at low p.
• xi_star_mode selects how that is treated. :interval (the default) spreads each pair over the table slices its ξ* range covers, with the weights of KnockOnKernels; the table then covers the ξ* of the interval ends of all contributing pairs, so that no pair is clamped and the production stays exact. :midpoint evaluates ξ* once, at the midpoint of the interval, and :centre once at the knock-on cell centre (DREAM). :midpoint keeps ξ* < 1. :centre does not on a knock-on grid coarser than the primary grid, where a cell centre can lie above p1: a pair with ξ* > 1 in the table range is refused with an ArgumentError, since no pitch distribution exists for it. The table of the two single-ξ* modes covers the ξ* of all contributing pairs.
• dream_compatible = true reproduces DREAM (MollerDeltaAngleKernel.cpp): ξ* at the knock-on cell centre, and a table range from the pairs with γ(p_{i-1/2}) < γ(p1_k) - 1 (a total energy compared with a kinetic energy). For primaries below γ₁ ≈ 3 this misses contributing pairs, whose ξ* is then clamped to the table edge; with a knock-on grid coarser than the primary grid, a cell centre above p1 gives ξ* > 1, which the range criterion excludes for most grids and clamps.
• The knock-on pitch grid must span [-1, 1]; the primary pitch grid may cover part of it.
• primary_pitch: how the primary pitch enters. :centre evaluates it at the cell centre (DREAM), which breaks reciprocity and isotropy at O(1) in the max norm. :split averages it over the cell, split at the kink points (pitch_kernel_column_split!), with primary_pitch_points Gauss points per piece: isotropic primaries then give isotropic knock-ons to quadrature accuracy. :gauss averages with primary_pitch_points plain Gauss points over the whole cell (slow convergence; kept for comparison).

pitch_cell_fraction(lower_face, upper_face, xi1, xi_s)

Fraction of knock-ons with pitch in [lower_face, upper_face] produced by a primary with pitch xi1 (ξ₁) at scattering-angle cosine xi_s (ξ*), averaged over the gyrophase (DREAM D, AnalyticDelta::EvaluateXiIntegral). The argument of asin is clamped to [-1, 1], so no separate reachability test is needed. When the distribution degenerates to a point (b² ≤ 5ε), the point is assigned to the half-open cell (lower_face, upper_face], as in DREAM.

pitch_kernel_column!(column, xi_faces, xi_s; primary_cell, nodes=(0.0,), weights=(2.0,))

Fill column[j] = δ_j for knock-on pitch cells with faces xi_faces, normalized so that Σ_j Δξ_j δ_j = 1, and return the norm before normalization. The primary pitch is averaged over primary_cell = (xi1_lower, xi1_upper) with Gauss nodes and weights on [-1, 1]; with primary_cell = (xi1, xi1) and the default single node the primary sits at xi1 (DREAM). A primary at xi1 = -1 whose knock-on pitch is exactly -1 is assigned to the first cell.

pitch_kernel_column_split!(column, xi_faces, xi_s; primary_cell, nodes, weights)

Like pitch_kernel_column! with the primary averaged over primary_cell = (xi1_lower, xi1_upper), but exact to high order: as a function of the primary pitch ξ₁, the fraction D_j(ξ₁) reaching knock-on cell j has square-root kinks where an edge a ± b of the gyro-averaged distribution crosses a face t of the knock-on grid, i.e. at ξ₁ = cos(arccos t ± θ*) with cos θ* = ξ* (and jumps at the faces when ξ* = 1). The cell is split at these points, and each piece is integrated with the Gauss nodes and weights (on [-1, 1]) in the variable s of ξ₁ = mid + h sin(πs/2), which turns the square-root behaviour at the piece ends into a smooth one.

pitch_kernel_matrix!(D, kernel::MollerPitchKernel, i, k)
pitch_kernel_matrix!(D, kernels::KnockOnKernels, i, k)

Pitch kernel of the momentum pair (i, k): from a MollerPitchKernel the kernel interpolated at the pair's single ξ, `D[j, l] = δ[j, l](ξ(pi, p1k)); from [KnockOnKernels](@ref) the effective kernel the operator applies, the table slices of the pair weighted by its entries and normalized by the pair's production, which forxistarmode = :intervalis the kernel integrated over ξ*(p) across the receiving interval and otherwise the same interpolation.rates[i, k]times this matrix is what [apply!`](@ref) contributes for the pair.

KnockOnKernels{T}

Energy and pitch kernels for one cutoff momentum, together with the production each momentum pair (i, k) sends to each slice of the ξ* table. The entries of pair (i, k) are pair_start[i, k]:pair_stop[i, k] (empty when the pair produces no knock-ons); entry t puts the weight entry_weight[t] [m³/s] on table slice entry_slice[t], and the weights of a pair sum to energy.rates[i, k], so the production is exactly the energy kernel's.

With xi_star_mode = :interval the weights come from the Møller flux integrated over the receiving interval against the hat functions φ_m of the table,

\[\mathtt{entry\_weight}[t] = \frac{\beta_1}{\Delta p_i} \int \beta(p) \, \Sigma(p, p_{1,k}) \, \varphi_{m_t}\bigl(\xi_*(p, p_{1,k})\bigr) \, \mathrm{d}p,\]

which is the pitch kernel integrated over ξ(p) instead of evaluated at one ξ. With :midpoint or :centre a pair has the two entries of the linear interpolation at its single ξ* (one when the interpolation weight is 0 or 1).

For the :blas evaluation, the primary momenta whose pairs use slice m are slice_primaries[m], their products occupy columns slice_offsets[m] + 1 : slice_offsets[m] + length(slice_primaries[m]) of the operator's slice_products, and entry_column[t] is the column of entry t.

KnockOnOperator{T}

Knock-on source on the grid knockon produced by primaries on the grid primary (the same grid for in-grid knock-on) colliding with stationary targets of density n_tot [m⁻³] (the total electron density):

\[\begin{aligned} C_{ijr} &= \frac{n_{\mathrm{tot},r}}{\mathcal{V}'_{K,ijr}} \sum_k S_{ik} \sum_l W_{klr} \, \delta_{jl}\bigl(\xi_*(p_i, p_{1,k})\bigr), \\ W_{klr} &= \Delta p_k \, \Delta \xi_l \, \mathcal{V}'_{P,klr} \, f_{klr}. \end{aligned}\]

C has the units of $\partial f / \partial t$. Its density moment is the total production rate, with total_rates from the energy kernel:

\[\frac{1}{V'} \sum_{ij} \mathcal{V}'_{ij} \Delta p_i \Delta \xi_j C_{ij} = n_\mathrm{tot} \sum_{kl} \frac{\mathcal{V}'_{kl} \Delta p_k \Delta \xi_l f_{kl}}{V'} \, \mathtt{total\_rates}[k].\]

Only the ratio 𝒱'_P / 𝒱'_K at equal radius enters, so any constant normalization of 𝒱' cancels (DREAM's cylindrical Vprime is 4π² times Runes2P's).

The cutoff momentum p_cutoff[r] (in units of mₑc) may differ between radii and change in time (set_cutoff!); each radius uses the kernels of its cutoff, and radii with equal cutoffs share them. Both grids must have the same radial grid. The operator is linear in f and in n_tot; DREAM evaluates it explicitly, from f at the previous time level, with one cutoff for all radii and times.

KnockOnOperator(knockon, primary, p_cutoff; method=:blas, table_size=80,
                dream_compatible=false, primary_pitch, primary_pitch_points=4)

Build the energy and pitch kernels for knock-on grid knockon and primary grid primary (both KineticGrid) with cutoff momentum p_cutoff in units of mₑc (a number, or one value per radius). method selects the evaluation (:blas or :loop, see the source file); the other keyword arguments are passed to MollerPitchKernel. dream_compatible = true, table_size = 80 (primary at the cell centre) with one cutoff reproduces DREAM's MollerBoltzmannOperator in cylindrical geometry.

apply!(df, op::KnockOnOperator, f, total_electron_density; scale=1)

Add scale times the knock-on source produced by the primary distribution f[k, l, r] to df[i, j, r] on the knock-on grid; total_electron_density[r] is the target density [m⁻³]. f and df are 3D arrays, so both grids must have a shared momentum grid; for per-radius grids use apply_at_radius!. Allocation free; not safe to call concurrently on the same operator (shared scratch).

apply_at_radius!(df, op::KnockOnOperator, f, total_electron_density, ir; scale=1)

Knock-on source at radius index ir alone: f is the primary distribution at that radius ((p cells, ξ cells) of the primary momentum grid there), df the rate of change of the knock-on distribution there, total_electron_density the target density [m⁻³] there. Works with a different momentum grid at each radius.

energy_kernel(op::KnockOnOperator, ir=1)

The MollerEnergyKernel used at radius index ir.

knockon_kernels(op::KnockOnOperator, ir=1)

KnockOnKernels used at radius index ir: the energy and pitch kernels together with the table slices and weights of every momentum pair.

knockon_linear_operator(op::KnockOnOperator, total_electron_density; scale=1)

Matrix-free linear operator scale ∂C/∂f of size (cell_count(knockon), cell_count(primary)), acting on vec(f) of the primary grid and returning vec of the knock-on grid: the implicit (Jacobian) form of the knock-on term, for use in Krylov solvers and in sums with other LinearOperators or matrices. C is linear in f, so the Jacobian-vector product is apply! itself. total_electron_density [m⁻³] is referenced, not copied: updating it in place (e.g. between Newton iterations) changes the operator. Both grids need a shared momentum grid, as apply! does; with a momentum grid per radius the products go through apply_at_radius!, as in KnockOnTerm.

The derivative with respect to the total electron density is C(f)/n_tot at each radius (DREAM's only Jacobian block for this term, which it treats explicitly in f). Treating f implicitly is a Runes2P departure; the explicit form stays available through apply!.

The product uses the operator's scratch arrays and is not safe to call concurrently. No transpose product is provided.

pitch_kernel(op::KnockOnOperator, ir=1)

The MollerPitchKernel used at radius index ir.

pitch_kernel_matrix!(D, kernels::KnockOnKernels, i, k)

Effective pitch kernel of the momentum pair (i, k), see the MollerPitchKernel method.

set_cutoff!(op::KnockOnOperator, p_cutoff)

Set the cutoff momentum in units of mₑc (a number for all radii, or one value per radius) and rebuild the kernels of the radii whose cutoff changed. Kernels are cached by the momentum grids of the radius and the cutoff value, so radii with the same grids and cutoff share one set, and a combination already in use costs nothing. Returns the number of kernel sets built. When a cutoff is invalid or a kernel cannot be built, the operator is unchanged. A rebuild costs a few operator applications; snapping the cutoff to momentum cell faces (choose_cutoff) keeps rebuilds rare when the cutoff follows the plasma.

A knockon_linear_operator of op sees the new kernels, and a KnockOnTerm of op recomputes its lagged source at its next preparation.

resolved_pitch_options(options::NamedTuple)

The kernel options of a KnockOnOperator with the defaults filled in: the keywords of MollerPitchKernel and xi_star_points of KnockOnKernels. The operator keeps the keywords it was given, and the record of a run needs the values used (describe).

choose_cutoff(p_faces, lower, upper; current_cutoff=nothing)

Cutoff momentum on a momentum-grid face (p_faces) inside the window [lower, upper), returned as (p_cutoff, inside_window), all momenta in units of mₑc. A current_cutoff that is a face of the grid and still inside the window is kept, so the kernels are rebuilt only when the window moves past it; a value that is not a face (another grid's cutoff) is replaced. Otherwise the face closest (in log p) to √(lower·upper) is taken, with lower raised to the smallest positive face when it lies below it (the smallest face ≥ lower when upper = Inf), which leaves room for the bounds to move in either direction. If no positive face lies inside the window (window closed or narrower than a cell), the largest positive face below upper is returned (keeping knock-ons above the critical momentum, at the cost of the lower bound), or the smallest positive face if none is below upper, with inside_window = false.

Snapping to faces also keeps the cutoff off the interior of a cell, so the energy kernel has no partial cell at the cutoff.

choose_cutoff_face(p_faces, target, lower, upper; current_cutoff=nothing)

The face of the momentum grid for a cutoff with the target target, inside the open window (lower, upper), returned as (p_cutoff, inside_window), all momenta in units of mₑc: the face nearest to the target in the logarithm of the energy. A current_cutoff that is a face inside the window is kept while the target lies between its two neighbouring faces, so a target that hovers between two faces does not rebuild the kernels at every time point; the cutoff is then at most one cell from the target. Without a face inside the window, the face nearest to the target is taken and inside_window is false.

critical_field(density, coulomb_logarithm)

Connor-Hastie critical field E_c = n e³ lnΛ / (4π ε₀² mₑ c²) [V/m] (Connor and Hastie, Nucl. Fusion 15, 415 (1975), doi:10.1088/0029-5515/15/3/007) for an electron density n [m⁻³] and a Coulomb logarithm lnΛ. In the limit of complete screening n is the density of the free electrons, and without screening the free plus the bound electrons (Hesslow et al., Nucl. Fusion 59, 084004 (2019), doi:10.1088/1741-4326/ab26c2; DREAM's Ec_free and Ec_tot, src/Equations/RunawayFluid.cpp:208-209).

critical_momentum(electric_field, critical_electric_field)

Critical momentum 1/√(E/E_c - 1) in units of mₑc, for the field electric_field (E) and the critical field critical_electric_field (E_c) (relativistic test particle, no pitch scattering or radiation), Inf for E ≤ E_c. The two fields are in the same unit, V/m in the package, and only their ratio enters. DREAM's effective critical momentum (fluid.pCrit) includes screening and pitch-angle scattering and is larger; a caller with a better estimate passes it to cutoff_window directly.

cutoff_target(temperature, critical_p; no_ceiling_factor=10)

Knock-on cutoff by the following rule (a choice of this package; see the docs page Choices of this package), returned as (p_cutoff, energy, policy) with p_cutoff in units of mₑc and energy [eV], for the critical momentum critical_p in units of mₑc:

• with a runaway threshold (critical_p finite), the logarithmic midpoint in energy between the thermal energy, which is the temperature [eV], and the kinetic energy at the critical momentum, E = √(T E_crit) (policy = :log_midpoint);
• without one (critical_p = Inf, the field below the critical field, where no avalanche is physical), no_ceiling_factor times the temperature (:no_ceiling);
• where the threshold is not above the thermal energy, which an estimate of the critical momentum can give in a hot plasma, the lower bound counts, since the operator heats the bulk when the cutoff comes near the thermal momentum: no_ceiling_factor times the temperature again (:closed).

no_ceiling_factor = 10 is a choice of this package (see the docs page Choices of this package).

cutoff_window(thermal_p, critical_p; thermal_factor=4)

Admissible range (lower, upper) = (thermal_factor · thermal_p, critical_p) of the knock-on cutoff momentum, all in units of mₑc: above the thermal momentum by a safety factor, below the critical momentum. The window is empty when lower ≥ upper. The factor thermal_factor = 4 is a choice of this package (see the docs page Choices of this package).

kinetic_energy(p)

Kinetic energy (γ - 1) mₑc² in eV of an electron with the momentum p in units of mₑc, in a form that keeps its digits for small p.

momentum_at_energy(energy)

Momentum in units of mₑc of an electron with the kinetic energy energy in eV, the inverse of kinetic_energy.

thermal_momentum(temperature)

Thermal momentum √(2T/mₑc²) in units of mₑc for the temperature temperature (T) in eV.