Neutral Particle Recycling

The native XGC neutral module is a Monte Carlo estimate of the neutral density, temperature, and associated reaction rates for ionization and charge exchange. The contribution from the neutrals to the ion and electron collision operators are thereby calculated.

Input parameters

Click here for default values and descriptions.

Monte Carlo calculation details

The neutrals translate in a 2D cartesian (R,Z) space against the toroidally-averaged background plasma. Sample contributions to density, energy, etc. are tallied every neu_dt seconds. For best accuracy, neu_dt ought to be at least as small as the time it takes for a fast neutral to traverse a mesh cell.

In each step, a random number is drawn and v*neu_dt (where v is the instantaneous neutral speed) is compared to the mean free path of ionization and charge exchange to determine if one of those reactions occur during that step.

When an ionization reaction occurs, the sample trajectory is terminated. When a charge exchange reaction occurs, the neutral attains an energy equal to the local ion temperature. The collision frequency for charge exchange in the native model is independent of neutral velocity (relative to the ion velocity).

Reation rates and sources

XGC models ionization and charge exchange reactions and their impacts on plasma density and energy.

Ionization creates ions at the neutral temperature, and does several things to the electrons. Firstly, it adds a cold electron at a temperature determined by a fit with respect to the electron energy \(E_e\):

\[T_c = E_{ion} / 6\]

for \(E_e>1.5 E_{iz}\) (with \(E_{iz} = 30\) eV; the energy loss per ionization), and

\[T_c = E_e/2 - E_{iz}/3\]

otherwise and is enforced to always be at least 3 eV. Note that an energy-dependent temperature precludes straightforward integration and this needs to be done numerically in-situ.

The source term of “cold” electrons (delta_fe_neut_ion_c2) is:

\[\left[ \delta f_e (\mathbf{v} ) / \delta t \right]_{esource,iz} = n_n \left\langle \sigma v \right\rangle_{iz} \sqrt{T_e} \mu(\mathbf{v}) \int T_c(\mathbf{v'})^{-3/2} f_e(\mathbf{v'}) e^{-E_e(\mathbf{v}) / T_c(\mathbf{v'}) } d^3\mathbf{v'}\]

(Correct for normalizations). This term represents the electron ejected from the atom at a cold temperature \(T_c\). A beam of monoenergetic electrons will produce a Maxwellian of secondary electrons at this temperature (which depends on the energy of the beam). However, there’s an entire distribution of electrons \(f_e\), which needs to be averaged over. This is similar to, but does not exactly correspond, to the treatment of Wersal and Ricci (2015). The sink term of hot electrons (delta_fe_neut_ion_h) is:

\[\left[ \delta f_e (\mathbf{v} ) / \delta t \right]_{sink,iz} = -n_n \left\langle \sigma v \right\rangle_{iz} f_e(\mathbf{v})\]

with a corresponding warm source (delta_fe_neut_ion_h) with energy reduced by \(E_{iz}\)

\[\left[ \delta f_e (\mathbf{v'} ) / \delta t \right]_{eloss,iz} = n_n \left\langle \sigma v \right\rangle_{iz} f_e(\mathbf{v})\]

where \(\mathbf{v'}\) is a velocity with random pitch angle, but an energy equal to \(E_e(\mathbf{v}) - E_{iz}\) (with a minimum of 1 eV unless \(E_e\) itself is below 1 eV. Since \(\mathbf{w}\) is partially random, it will be between velocity space nodes so its contribution needs to be appropriately weighted to neighboring velocity nodes. (isn’t the velocity-space representation finite-volume, so the velocity would instead fall into a unique bin?).

Finally, the source of ions from ionization (delta_fi_neut_ion):

\[\left[ \delta f_i (\mathbf{v} ) / \delta t \right]_{ion,iz} = n_e \left\langle \sigma v \right\rangle_{iz} f_n(\mathbf{v})\]

Charge-exchange has a similar logic, but applied only to the ions. An ion sink at one energy becomes a ion source at the neutral temperature.

The charge-exchange sink term for ions (delta_fi_neut_cx) is:

\[\left[ \delta f_i / \delta t \right]_{sink,cx} = -n_n \left\langle \sigma v \right\rangle_{cx} f_i\]

while the source (delta_fi_neut_cx_M) is:

\[\left[ \delta f_i / \delta t \right]_{source,cx} = \left\langle \sigma v \right\rangle_{cx} f_n\]

this is then rescaled so that its phase-space integral matches that of delta_fi_neut_cx to enforce conservation (why is there no MPI reduction over these sums? at this stage?)

References