Collision models

The distribution $f(p, ξ, r)$ of the electrons obeys a Fokker-Planck equation in the momentum $p$ (units of $m_e c$) and the pitch $ξ = p_\parallel/p$. The equation holds the electric field, the synchrotron loss and, for small-angle collisions, the linearized relativistic test-particle operator of DREAM. Large-angle collisions enter as a separate knock-on source, and a field-particle term restores the momentum the test-particle operator loses. The geometry is a cylinder in a homogeneous magnetic field, so the bounce averages are trivial: each coefficient takes its local value.

The Fokker-Planck operator

fokker_planck_operator fills the face coefficients of an AdvectionDiffusion operator in the flux form

\[\frac{\partial f}{\partial t} = -\frac{1}{\mathcal{V}'} \frac{\partial}{\partial z^k} \left[\mathcal{V}' \left(A^k f - D^{kl} \frac{\partial f}{\partial z^l}\right)\right], \qquad z = (p, ξ),\]

with the phase-space Jacobian $\mathcal{V}'$ (phase_space_jacobian). Each physical process adds to the coefficients:

fokker_planck_operator and set_fokker_planck_coefficients! write the slowing down and the energy diffusion in one sweep over the p faces. Their coefficients are those of add_slowing_down! and add_energy_diffusion! bit for bit. They build their CollisionModel from the faces of the grid. With such a model, they and add_pitch_scattering! read the factors that depend on the momentum and the charge state alone from tables built once per momentum grid. A model built without the faces evaluates these factors directly. These factors are the bremsstrahlung cross section and the bound-electron terms $w_j h_j$ and $w_j g_j$ below. A tabulated factor is the number of the direct evaluation, so the coefficients do not change.

The collision frequencies $ν_s$, $ν_D$ and $ν_\parallel$ come from a CollisionModel of the plasma and the CollisionSettings. The flux through $p = 0$ follows from $\lim_{p \to 0} p^3 ν_s$ (p3_slowing_down_frequency_at_zero). p_zero_boundary = :zero_flux of fokker_planck_operator takes no flux there, as DREAM's runaway grid. The interpolation of $f$ to the faces is the scheme of the operator (see Positivity and flux limiters).

An electric field that is an unknown of the system does not enter the coefficients. The operator then holds the coefficients of a unit field (set_unit_field_coefficients!), and the FokkerPlanckTerm multiplies them by the field of each radius. A plasma that changes in time rebuilds the coefficients through set_fokker_planck_coefficients!, with the scheme and the boundary conditions kept (prescribed_fokker_planck_term). With a plasma read from the unknowns (PlasmaFromUnknowns) the term is nonlinear, and its keyword lagged selects the part of the plasma taken at the start of the step. The default lagged = :ions is a choice of this package (see Choices of this package).

Collision frequencies

The frequencies are those of DREAM's CollisionFrequency classes, with the partial screening of L. Hesslow et al., J. Plasma Phys. 84, 905840605 (2018), doi:10.1017/S0022377818001113. The mode of CollisionSettings sets the target electrons:

• :full: relativistic Maxwellian targets at the cold-electron temperature, with the thermal factors of Pike and Rose and the energy diffusion $ν_\parallel$.
• :superthermal: cold targets, no energy diffusion.
• :ultra_relativistic: cold targets, no energy diffusion, and the limit $v \to c$ of the prefactors.

The thermal factors are those of O. J. Pike and S. J. Rose, Phys. Rev. E 89, 053107 (2014), doi:10.1103/PhysRevE.89.053107.

With $Θ = T/(m_e c^2)$, the cold-electron density $n_c$ and the density $n_j$ of charge state $j$ (charge $Z_{0j}$, nuclear charge $Z_j$), the frequencies of the :full and :superthermal modes are

\[\begin{aligned} ν_s &= 4π r_0^2 c\, \frac{γ^2}{p^3} \Big[\ln Λ_{ee}\, ψ_s(p)\, n_c + \sum_j n_j\, w_j(p)\, h_j(p)\Big] + ν_{s,\mathrm{rad}}, \\ ν_D &= 4π r_0^2 c\, \frac{γ}{p^3} \Big[\ln Λ_{ee}\, ψ_D(p)\, n_c + \ln Λ_{ei} \sum_j n_j Z_{0j}^2 + \sum_j n_j\, w_j(p)\, g_j(p)\Big], \\ ν_\parallel &= Θ γ\, ν_s^{\mathrm{coll}}, \end{aligned}\]

slowing_down_frequency, deflection_frequency and parallel_diffusion_frequency evaluate them. With screening = :none, $n_c$ includes the bound electrons and $Z_{0j}^2$ becomes $Z_j^2$. $ψ_s$ and $ψ_D$ are the thermal factors of the :full mode and one otherwise. $h_j$ and $g_j$ are the bound-electron terms of partial screening, $w_j$ their screening weight, and $ν_{s,\mathrm{rad}}$ the bremsstrahlung stopping. $ν_s^{\mathrm{coll}}$ is $ν_s$ without $ν_{s,\mathrm{rad}}$, unless bremsstrahlung_diffusion = :dream. With screened_diffusion = :zero the energy diffusion takes the free-electron part alone.

Coulomb logarithms

coulomb_logarithm gives the logarithm of electron-electron and electron-ion collisions in three models:

• :thermal: $\ln Λ_T = 14.9 - \tfrac12 \ln\bigl(n/(10^{20}\,\mathrm{m^{-3}})\bigr) + \ln(T/\mathrm{keV})$ (thermal_coulomb_logarithm).
• :relativistic: the momentum-independent $\ln Λ_c = 14.6 + \tfrac12 \ln\bigl((T/\mathrm{eV})/(n/(10^{20}\,\mathrm{m^{-3}}))\bigr)$ (relativistic_coulomb_logarithm).
• :energy_dependent: $\ln Λ = \ln Λ_T + \tfrac{1}{k} \ln(1 + x^k)$ with $k = 5$, $x = \sqrt{2(γ - 1)}/p_T$ for electrons and $x = 2p/p_T$ for ions, $p_T = \sqrt{2Θ}$.

These are eqs. 2.7 to 2.10 of Hesslow et al. (2018), as in DREAM's CoulombLogarithm.cpp. $n$ is the free-electron density.

Screening

The screening of CollisionSettings sets how the bound electrons of partially ionized ions act on a fast electron:

• :complete: the bound electrons are invisible, and an ion scatters with its net charge $Z_{0j}$.
• :none: the bound electrons are targets like free ones, and an ion scatters with its nuclear charge $Z_j$.
• :partial: the model of Hesslow et al. (2018), between the two.

Partial screening adds the Bethe stopping on the bound electrons, $h_j$, with the mean excitation energy $I_j$ of the charge state (mean_excitation_energy). It adds the deflection $g_j$ by the screened nucleus, with the effective ion size of the charge state (ion_size). The sizes are the DFT values of Table 1 of Hesslow et al. where tabulated and their eq. 2.28 otherwise. This is ion_size_table = :hesslow (the default); ion_size_table = :dream takes DREAM's lookup of that table. $h_j$ takes DREAM's form $\ln(1 + (h\,e^{-β^2})^k)/k$ of eq. 2.31 of Hesslow et al. (2018), with $h = p\sqrt{γ - 1}/I_j$ and $I_j$ in units of $m_e c^2$. This form stays positive where the paper's $\ln(1 + h^k)/k - β^2$ turns negative below $h ≈ 0.23$ ($p ≈ 0.012$ for Ar⁺). The two agree to 3e-9 for $h ≥ 16$, for Ar⁺ at $p ≥ 0.1$.

Fig. 3(b) of Hesslow et al. (2018) shows the deflection frequency with partial screening over a deflection frequency without screening. The curve of the figure and the text of the paper (p. 12) differ in that denominator. Read from its curve, the frequency without screening keeps only the free electrons among the targets of the electron-electron collisions. The text counts the bound electrons among these targets, as screening = :none and DREAM's NON_SCREENED do.

test/test_independent_references.jl divides the package's frequency with partial screening by the denominator of the curve. It compares that ratio with the readings of Fig. 3(b), for Ar⁺ at 10 eV with $n_e = n_\mathrm{Ar} = 10^{20}$ m⁻³. The readings are at four momenta, from $p = 0.1$ to $p = 100$. Each difference is relative to the reading at the same momentum.

Both values are measured. The test asserts the first ratio against the readings within 3 %. The second is recorded only in a comment of the test.

With screening_window = :excitation_energy (the default) the bound-electron parts are weighted by $w_j$, which is a cubic smooth step in the kinetic energy $E_k = (γ - 1) m_e c^2$:

\[w_j = s^2 (3 - 2s), \qquad s = \operatorname{clamp}\!\left(\frac{E_k - a I_j}{(b - a) I_j},\, 0,\, 1\right),\]

with $a$ and $b$ the settings screening_window_start and screening_window_end. Below the window the frequencies are those of complete screening. The default factors are a choice of this package (see Choices of this package). screening_window = :dream takes $w_j = 1$ at all momenta, as DREAM does.

With bound_friction = :ionization_split the stopping on the bound electrons splits into two parts. The ionization part is the energy of the ionizations that $f$ causes. The excitation part is the remainder, floored at zero. The ionization part uses the cross sections of the kinetic ionization (ionization_cross_section), so that $f$ pays the energy of the ionizations that the KineticIonizationTerm counts (see Atomic data and ionization).

Bremsstrahlung

With bremsstrahlung = :stopping_power the radiative energy loss on the nuclei of all charge states enters $ν_s$ as a drag:

\[ν_{s,\mathrm{rad}} = \sum_j n_j Z_j^2\, α r_0^2 c\, \frac{γ}{p}\, φ_\mathrm{rad}(p),\]

with the unscreened cross section 4BN of H. W. Koch and J. W. Motz, Reviews of Modern Physics 31, 920 (1959), doi:10.1103/RevModPhys.31.920 (bremsstrahlung_cross_section). There is no electron-electron bremsstrahlung, as in DREAM. The stopping is on by default, where DREAM leaves it out, a choice of this package (see Choices of this package). bremsstrahlung_constant selects the constant term of the formula: $4/3$ (Koch and Motz) with :koch_motz (the default), or $1$, DREAM's integer value, with :dream.

bremsstrahlung_diffusion decides whether the energy diffusion of the :full mode takes the radiative stopping as well:

• :none (the default): the stopping is a pure drag, as in the model of M. Hoppe et al., Comput. Phys. Commun. 268, 108098 (2021), doi:10.1016/j.cpc.2021.108098, eqs. 8, 9, 14 and 15.
• :dream: DREAM's code, whose $ν_\parallel$ rescales the whole $ν_s$.

With :none a Maxwellian at the bulk temperature loses its bremsstrahlung power, which tends to Gould's formula at low temperature (GOULD_BREMSSTRAHLUNG_PREFACTOR). The formula is that of R. J. Gould, The Astrophysical Journal 238, 1026 (1980), doi:10.1086/158068. With :dream a Maxwellian at the bulk temperature radiates nothing, because the energy diffusion rescaled with the whole $ν_s$ balances the radiative drag on it. The default :none is a choice of this package (see Choices of this package).

Synchrotron loss

add_synchrotron! adds the radiation reaction of the gyration in the magnetic field $B$ of each radius, with

\[C = \frac{e^4}{6π ε_0 m_e^3 c^3}\]

(SYNCHROTRON_RATE_PER_FIELD_SQUARED). It is on by default in fokker_planck_operator, set_fokker_planck_coefficients! and the FokkerPlanckTerm, where DREAM leaves it out. The default is chosen for accuracy. synchrotron = false leaves the loss out.

Knock-on collisions

The small-angle operator holds the collisions of $f$ with the thermal electrons. Close collisions of a fast primary with an electron at rest produce secondaries at large momentum transfer, the knock-ons. A KnockOnOperator adds their source, the Boltzmann gain term of the cross section of C. Møller, Annalen der Physik 406, 531 (1932), doi:10.1002/andp.19324060506, above a cutoff momentum $p_\mathrm{cut}$. It is ported from DREAM's MollerBoltzmannOperator. The target density is the total electron density $n_\mathrm{tot}$, free and bound. The operator is the source of the knock-ons alone. It has no loss term for the primaries.

The source in a knock-on cell is $n_\mathrm{tot}$ times a sum over the primary cells of three factors: an energy kernel, a pitch kernel and the number of primaries in the cell (the formula is in KnockOnOperator).

• The energy kernel $S_{ik}$ (MollerEnergyKernel) is the Møller production rate from primary momentum cell $k$, averaged over knock-on cell $i$. It integrates the cross section (moller_cross_section, moller_flux) from the cutoff to half the kinetic energy of the primary (max_knockon_momentum).
• The pitch kernel $δ_{jl}$ (MollerPitchKernel) is the distribution in knock-on pitch of the knock-ons from primary pitch cell $l$, normalized so that $\sum_j Δξ_j δ_{jl} = 1$, from the gyrophase-averaged fraction of each cell (pitch_cell_fraction). It is taken at the scattering cosine $ξ^*$ that the kinematics fixes (xi_star).

The pitch kernel is tabulated in the scattering cosine and interpolated per momentum pair (KnockOnKernels). By default the primary pitch is averaged over its cell, and the scattering cosine is integrated over the momentum interval that receives the knock-ons. With dream_compatible = true both are taken at the cell centres, as DREAM does. The primary and the knock-on grids may differ, with the same radial grid, and the momentum grid may differ between radii (apply_at_radius!).

A KnockOnTerm places the source in the equation of a kinetic unknown. Its treatment takes three values. :explicit_primary (the default) takes the primary at the start of the step, as DREAM's Boltzmann operator does. :explicit takes the whole term, target density included, at the previous time level. :implicit takes the whole term at the new time level, matrix-free, since it couples all cells of a radius. The operator's own matrix-free Jacobian is knockon_linear_operator.

The cutoff

The cutoff separates the knock-ons counted by the source from the collisions of the small-angle operator. It may differ between radii and change in time: set_cutoff! takes any positive cutoff per radius, and the kernels (MollerEnergyKernel, MollerPitchKernel) take the part of a cell above a cutoff inside it. A KnockOnCutoffUpdate sets it after every accepted step from the plasma of each radius, on a face of the momentum grid, so that no cell is cut. Its default placement (cutoff_target) is the logarithmic midpoint in energy between the temperature $T$ and the kinetic energy $E_\mathrm{crit}$ at the critical momentum:

\[E_\mathrm{cut} = \sqrt{T\, E_\mathrm{crit}}.\]

Where no critical momentum exists (the electric field below the critical field), or its kinetic energy is not above the temperature, the cutoff energy is no_ceiling_factor times the temperature. A face is kept while the target stays between its neighbouring faces (choose_cutoff_face), so the kernels are rebuilt rarely. The placement and no_ceiling_factor are choices of this package (see Choices of this package). The threshold is total_critical_momentum by default. Its critical field takes the Coulomb logarithm named by the keyword coulomb_logarithm of KnockOnCutoffUpdate, as DREAM's critical fields do. The default :thermal gives the thermal logarithm, and the other models the relativistic one. A run whose collisions take another model passes it, and :relativistic is allowed in any run. DREAM keeps one cutoff for all radii and times.

With a hot-tail and a runaway grid the pairings of the two grids share one cutoff, through the keyword linked of KnockOnCutoffUpdate, so that every knock-on above it is counted once.

Fluid avalanche

A fluid runaway density above a grid creates knock-ons as well. RosenbluthPutvinskiSourceTerm is DREAM's kinetic Rosenbluth-Putvinski source, the knock-ons of a primary at the speed of light (rosenbluth_putvinski_rate). It is the source of M. N. Rosenbluth and S. V. Putvinski, Nucl. Fusion 37, 1355 (1997), doi:10.1088/0029-5515/37/10/I03. ExternalAvalancheTerm is the avalanche growth rate of the fluid runaways that DREAM uses above the grid. The rate is that of L. Hesslow et al., Nucl. Fusion 59, 084004 (2019), doi:10.1088/1741-4326/ab26c2, with the fluid_avalanche_growth_rate of the effective_critical_field and the avalanche_critical_momentum. The fluid growth rate takes the collision settings of DREAM's RunawayFluid (runaway_fluid_settings), not those of the kinetic equation.

Field-particle term

The test-particle operator treats the thermal electrons as a fixed Maxwellian. It removes the parallel momentum that the electrons give to each other, which lowers the conductivity. A FieldParticleTerm returns that momentum to the bulk with a model field-particle operator (FieldParticleOperator):

\[C_\mathrm{fp}(f) = \sum_{a,b} g_a A_{ab} (\ell_b \cdot f), \qquad a, b = 0, \dots, N - 1,\]

per radius, with $N$ = modes. $\ell_b \cdot f$ is the rate at which the discrete electron-electron operator (electron_electron_operator) changes the parallel momentum weighted by the Laguerre polynomial $L_b^{(3/2)}(x^2)$, $x = v/v_T$. The polynomials are those of the expansion of S. I. Braginskii, Rev. Plasma Phys. 1, 205 (1965). $g_a \propto v_\parallel L_a^{(3/2)}(x^2) f_M$ is the shape of mode $a$ on the Maxwellian $f_M$. $A$ fixes the matrix elements to those of the Landau operator between the first modes Laguerre modes. The improved Sugama operator does the same (H. Sugama et al., Phys. Plasmas 26, 102108 (2019), doi:10.1063/1.5115440). Row 0 of $A$ makes the parallel momentum conserved to round-off. With modes = 1 the term is the momentum-restoring model of S. P. Hirshman and D. J. Sigmar, Phys. Fluids 19, 1532 (1976), doi:10.1063/1.861356, and the higher modes add the heat-flux corrections.

The term needs the :full mode, where the bulk is part of $f$. Its matrix elements are those of the nonrelativistic Landau operator, so it holds for a nonrelativistic bulk. It couples all cells of a radius and is applied matrix-free, with a Jacobian of rank modes per radius (low_rank_jacobian) that the linear solvers take by a Woodbury update. set_field_particle_coefficients! rebuilds it for a plasma that changes in time.

Ion energy exchange

The thermal exchange between the electrons and the ion species is the IonEnergyExchangeTerm (see Energy and current balance). Its electron-ion logarithm is the thermal one, as in DREAM. Among the ions, ion_coulomb_logarithm = :nrl (the default) takes the mixed ion-ion logarithm of the NRL Plasma Formulary (A. S. Richardson, NRL/PU/6770–19-652, 2019). It takes the temperature of each species and a screening sum over all ion species and the electrons. It is cut at zero so that the exchange keeps its sign. The cut guards against cold iterates of the solver and cold prescribed input. The docstring of IonEnergyExchangeTerm gives examples of the logarithm for deuterium and argon, in equilibrium and with Ar¹⁶⁺ left cold. :dream takes DREAM's form, built from the electron temperature and the same for every pair.

Comparisons

• The collision frequencies against DREAM: test/test_collisions_dream.jl.
• Partial screening against eqs. 2.22, 2.25 and 2.31 and Figs. 3 and 4 of Hesslow et al. (2018), Ar⁺ at 10 eV: test/test_independent_references.jl.
• One step of the small-angle kinetic operator against DREAM: test/test_kinetic_operator_dream.jl.
• The knock-on operator in the DREAM-compatible configuration against DREAM's MollerBoltzmannOperator: test/test_knockon_operator.jl.
• The default knock-on operator against a quadrature reference, and its physics limits: test/test_knockon_reference.jl, test/test_knockon_physics.jl.
• An avalanche run against DREAM's Boltzmann knock-on mode: test/test_avalanche_dream.jl.
• The Rosenbluth-Putvinski source and the fluid growth rate against DREAM: test/test_avalanche_rp_dream.jl.
• The growth rate of the Rosenbluth-Putvinski source against that of the knock-on operator in one avalanche run, whose primaries have a mean momentum of $34\, m_e c$: test/test_runaway_physics.jl. The source gives 1.04 times the growth rate of the knock-on operator at the step of the test (1.03 extrapolated to a zero step).
• The Dreicer rate with the self-consistent bulk against that with a fixed bulk: test/test_populations.jl (below).
• The field-particle term, test/test_field_particle.jl: momentum conservation, the matrix elements of the test-particle operator against Braginskii's, and the field-particle ratios against the like-species sums. The same file checks the Lorentz limit and the Spitzer-Härm conductivity. The Spitzer-Härm conductivity is that of L. Spitzer and R. Härm, Phys. Rev. 89, 977 (1953), doi:10.1103/PhysRev.89.977. The file also checks the conductivity at $Z = 1$ from 1 to 10 keV against that of B. J. Braams and C. F. F. Karney, Phys. Fluids B 1, 1355 (1989), doi:10.1063/1.858966. It is 1.1 % above theirs at 10 keV, a departure that grows with the temperature and that the nonrelativistic field-particle matrix elements would explain.
• The field of the first step against Braginskii: test/test_first_step_field.jl (below).

The Dreicer rate with a self-consistent bulk

test/test_populations.jl compares the Dreicer rate under two models of the thermal bulk. MaxwellianSourceTerm drives the density below $p = 0.12$ back to its initial value. It restores at each step what the previous step lost, so that the bulk stays near its initial density. This model is the fixed bulk, which the published Dreicer rates assume.

ParticleSourceTerm holds the density on the grid at the free-electron density less the runaway density. The flux through the upper momentum boundary of the grid feeds the runaway density. With this self-consistent bulk, as in DREAM, the electrons that ran away are missing from the bulk.

The lost fraction is the share of the electrons that have run away by the end of the self-consistent run, 0.78 %. The depleted run is a run with MaxwellianSourceTerm whose start and target are scaled by 1 − the lost fraction.

The depleted run lacks the lost electrons from the start. The self-consistent bulk loses them over the run. So 0.78 % bounds the effect of the depletion from above. At least 3.0 % of the drop therefore comes from the difference of the two source terms.

The source terms differ in the density they hold, that of the whole grid or that below $p = 0.12$. They also differ in their shapes, a relativistic or a nonrelativistic Maxwellian. The test does not trace the share of each.

What the test asserts:

• The self-consistent rate lies between 0.93 and 0.99 of the rate with the fixed bulk.
• The lost fraction lies between 0.005 and 0.02.
• The depleted rate equals (1 − the lost fraction) times the rate with the fixed bulk, to a relative 1e-9.
• The self-consistent rate lies at least 0.02 below the depleted rate, both relative to the rate with the fixed bulk.

The field of the first step

test/test_first_step_field.jl takes one short step of a system that holds the current of a shifted Maxwellian. A constraint on the current moment holds the current, with the field as its unknown. The field of the step balances the friction of the collisions. The test compares it with $E_B = m_e j/(n e^2 τ_e)$ of S. I. Braginskii, Rev. Plasma Phys. 1, 205 (1965), eqs. (2.5e), (2.7) and (4.6). Here $τ_e$ is Braginskii's electron collision time, his eq. (2.5e).

At 120 eV, $τ_e$ is 6.89e-7 s and $E_B$ is 1.030 V/m. The extrapolated ion rows compare with $E_B (1 + 13Θ/8)$, with $Θ = T/(m_e c^2)$. The header of the test derives that correction to first order in $Θ$. In the nonrelativistic limit, the test-particle ratio is $1 + 1/(\sqrt{2} Z)$, derived in the test from the slowing-down frequencies.

With FieldParticleTerm, the electron-electron collisions conserve the parallel momentum. The heat-flux modes of the term carry no current, so the field stays at that of the ions alone, the reference of its rows. That argument is nonrelativistic.

Each deviation is relative and signed: (value − reference) / reference. The bound of the test is on the same quantity.

• For the extrapolated ion rows, the test gives deviation and bound as absolute values of $E/E_B - (1 + 13Θ/8)$, which agree with the relative ones to the digits shown.
• The step biases the extrapolated field by about −3e-6, the level of the deviations of the extrapolated ion rows at 120 and 30 eV. So the bound of these rows is set by the step, not by the operator.
• Shorter steps than 1e-8 $τ_e$ are no better, since the round-off of the solve grows as 1/dt.
• At 120 and 30 eV, the deviations of the test-particle ratio are of the order of $Θ$, with an uncertainty of about 3e-6 from the step and the extrapolation.
• The :linear weights are exact on the first-order shift, which is linear in $ξ$. On the exact shift at 10 MA/m², the current with the :linear weights differs from the momentum that FieldParticleTerm conserves. The field then lies 2.2e-3 above that of the ions alone with 4 pitch cells, and 1.9e-4 above with 16.