2.1. Uniform loss with planar symmetry

The simplest case that nevertheless exhibits key features of the analysis to follow is ballistic transport from a source at the boundary through a one-dimensional domain with a uniform loss rate. Let $f(x,v)$ be the phase space distribution of neutrals produced with a velocity distribution $S(v)$ and eliminated at a loss rate $\gamma$ (in units of $\mathrm{s}^{-1}$ ). The kinetic equation describing this case is

(2.1) \begin{equation} v \frac {\partial f}{\partial x} + \gamma f = S(v) \delta (x), \end{equation}

with $\delta (x)$ representing the Dirac distribution. With an integrating factor, (2.1) admits an analytic solution

(2.2) \begin{equation} f(x,v) = \frac {S(v)}{v} e^{-\gamma x/v}. \end{equation}

While the distribution function can be written in closed form, its density moment $n_n \equiv \int f \,\mathrm{d}v$ cannot. This is of chief concern since the local ionisation rate density $R_{iz}$ (and other important quantities of interest) are directly proportional thereto: $R_{iz} = n_n \gamma _{iz}$ . In order to perform such an integration, one must choose a velocity dependence of the source. Throughout this work, the choice is an isotropic Maxwellian

(2.3) \begin{equation} S(v) =\frac {2 S_0}{\sqrt {\pi } v_{tn}} e^{-v^2/v_{tn}^2}, \end{equation}

where $v_{tn}$ is the thermal speed associated with the source such that $T_n = m_n v_{tn}^2 / 2$ (with neutral particle mass $m_n$ ), $S_0$ is the rate of particles produced per unit time per unit area, the other two velocity coordinates are already averaged upon and $S(v)$ is only defined for $v\gt 0$ . The Maxwellian source of (2.3) is not necessarily physically motivated. It would only be expected in the fluid limit where neutrals equilibrate by frequent charge exchange with ions locally. Instead, one expects a non-Maxwellian population of neutrals, with different populations interacting with ions of different temperatures before crossing the separatrix. This requires further analysis that depends on a detailed description of the scrape-off layer (SOL). Since an analytic form of the source needs to be chosen to estimate key quantities such as ionisation rate, this work will proceed with the Maxwellian source. In this case, the density moment is

(2.4) \begin{equation} n_n =\frac {2 S_0}{\sqrt {\pi } v_{tn}} \int \limits _0^\infty v^{-1} \exp \left (-\frac {v^2}{v_{tn}^2} - \frac {\gamma x}{v} \right ) \,\mathrm{d}v. \end{equation}

To approximate integrals of this type, this work will use the saddle point method, taking advantage of the relatively large loss rate $\gamma \gg v_{tn}/x$ . In (2.4), the integrand has a saddle point at $v_* = v_{tn} (\gamma x / 2 v_t )^{1/3}$ , and the integral in the asymptotic limit is

(2.5) \begin{equation} n_n(x) \approx \frac {2 S_0}{\sqrt {3} v_{tn}} \left ( \frac {\gamma x}{2 v_{tn}} \right )^{-1/3} \exp \left [ - 3 \left ( \frac {\gamma x}{2 v_{tn}} \right )^{2/3} \right ]\!. \end{equation}

This is equivalent to the asymptotic limit of the singular eigenfunction solution of Connor (Reference Connor1977).

To test the applicability of (2.5), several cases are prepared with DEGAS2 using electron-impact ionisation as the only reaction included. For the purposes of estimating neutral density, this reaction is purely a sink precisely as it acts in (2.1), and DEGAS2 directly solves along trajectories by sampling the source distribution and dividing the domain into many sub-volumes to estimate neutral tallies. In these cases, the neutral thermal velocity corresponds to a temperature $T_n =$ 100 eV. The electron density and temperature were chosen so that the three cases correspond to, within three significant figures, $\hat {\gamma } = (3, 30, 100)$ , where $\hat {\gamma } \equiv \gamma / v_{tn}$ is an inverse ‘mean free path’. This is a known to be an important quantity for fuelling (Fujii et al. Reference Fujii, Hasuo, Goto and Lore2025) and will remain dimensional with units of inverse length to stress that it does not directly scale with device size (Mordijck Reference Mordijck2020). Table 1 lists the properties of these test cases. Note that, by scanning in $\hat {\gamma }$ , not only are different plasma conditions accounted for, but also different source thermal velocities.

Table 1. List of the parameters used for the uniform-plasma cases considered here. Labels U1–U3 are the planar symmetry cases, and X1–X3 are the X-point cases. Also shown are the normalised rates for ionisation and charge exchange, respectively.

Figure 2. Predicted neutral density through a uniform-plasma domain with several normalised ionisation loss rates $ \hat {\gamma } = n_e \langle \sigma _{iz} v \rangle / v_{tn}$ and identical source rates. On the left (a) is a linear scale at low $x$ , and on the right (b) is a logarithmic scale. The crosses are results from DEGAS2, the solid lines represent the numerical evaluation of the integral (2.2) and the dashed line is the analytic closed-form approximation (2.5). Cases correspond to those labelled U1, U2 and U3 in table 1.

Figure 2 shows these three cases comparing the DEGAS2 results with the model of (2.5). Also shown are the results from direct numerical calculation of the integral in (2.4). It is worth noting that all these cases have the same neutral source: a rate of $10^{22}$ particles per second per $\mathrm{m}^2$ and a source temperature of 100 eV. Agreement between the model and simulation is generally good, with some exceptions near $x=0$ . Since this is the foundation upon which more general models are built, it is worth examining the discrepancy. DEGAS2 utilises ‘suppressed ionisation’ whereby it directly uses the analytic solution of $f$ in (2.2), so it should, in principle, exactly match the direct numeric integration of (2.4). The reason any disagreement between DEGAS2 and numerical integration is seen at all in figure 2 is because DEGAS2 calculates a volume-averaged tally, whereas the function (2.2) is integrated at the particular values of $x$ . Since $n_n$ has a strong second derivative near $x=0$ , the volume-averaged density in a cell is slightly more than the value evaluated at the midpoint of the cell. A similar issue was encountered in the benchmarks of Bernard et al. (Reference Bernard, Halpern, Francisquez, Mandell, Juno, Hammett, Hakim, Wilkie and Guterl2022), requiring DEGAS2 to use higher resolution to capture higher-order variability within spatial cells. Additional disagreement is seen between the closed-form saddle point solution (2.5) and both the DEGAS2 results and direct integration of the kinetic equation. This disagreement is due to the breakdown of the saddle point approximation itself, wherein $\hat {\gamma }$ is not reliably considered asymptotically large, especially near $x=0$ . Even so, a disagreement of only $\sim 25\,\%$ when the ‘large’ asymptotic parameter is much smaller than unity (approximately 0.03 for the $\hat {\gamma } = 3$ case at $x \approx 0.01$ ) is remarkably fortunate. This might be an indication that the analytic approximation of (2.4), and those like it that follow, can be derived by other means.

2.2. Pedestal loss with planar symmetry

The solution to (2.1) is readily found for any $\gamma (x)$ that admits a closed-form indefinite integral, but this work shall restrict itself to a particular functional form

(2.6) \begin{equation} \gamma (x) = \gamma _0 + \frac {\gamma _\infty - \gamma _0}{2} \left [ 1 + \tanh \left (\frac {x-x_0}{\varDelta } \right ) \right ]\!, \end{equation}

where $x_0$ and $\varDelta$ are parameters defining the shape of $\gamma (x)$ , while $\gamma _0$ and $\gamma _\infty$ are the left and right asymptotes of $\gamma$ , corresponding to the SOL and core asymptotes, respectively. This form is inspired by the plasma profiles typically fit in the pedestal region of tokamaks (Groebner & Osborne Reference Groebner and Osborne1998). The interpretation of $\varDelta$ and $x_0$ are the ‘pedestal’ half-width and the location of largest gradient, respectively. The loss profile shapes that will be used are listed in table 2. Note that the shapes and asymptotes of the internally fit loss-rate profile do not necessarily match that of the plasma itself, since the former is a tabulated function of the latter.

Table 2. List of the parameters determining the loss-rate profiles tested for the spatially varying planar cases. Here, $\varDelta _e$ and $x_{0e}$ correspond to the pedestal shape properties of the electron population, where $\varDelta$ and $x_0$ are the parameters for the fitted shapes of the ionisation loss-rate profile $\gamma (x)$ .

Under the functional form of (2.6), the solution to (2.1) is

(2.7) \begin{equation} f(x,v) = \frac {S(v)}{v} \exp \left ( \frac {x_0 - x}{v}\frac {\gamma _0 + \gamma _\infty }{2} \right ) \left \{ \frac { \cosh \left [{(x - x_0)}/{\varDelta } \right ]}{\cosh ({x_0}/{ \varDelta })} \right \}^{{\varDelta }/{v}\left ( \gamma _\infty - \gamma _0 \right )}. \end{equation}

The density moment under the same approximation as the previous section is

(2.8) \begin{equation} n_n(x) \approx \frac {2 S_0}{\sqrt {3 \pi }} \left ( \frac {u(x)}{2} \right )^{-1/3} \exp \left [- 3 \left ( \frac {u(x)}{2} \right )^{2/3}\right ]\!, \end{equation}

where

(2.9) \begin{equation} u(x) = \varDelta \frac {\gamma _\infty - \gamma _0}{2}\ln \left\{ \frac {\cosh \left [{(x -x_0)}/{\varDelta } \right ]}{ \cosh \left ({x_0}/{\varDelta } \right )}\right \} + \frac {\gamma _0+\gamma _\infty }{2}\left (x-x_0\right )\!. \end{equation}

With this model, the neutral density throughout a medium with a spatially varying loss rate can be estimated. Again, DEGAS2 is utilised for testing. Four artificial pedestals are constructed with the properties listed in table 2. In DEGAS2, electron density and temperature profiles are constructed similarly to (2.6). Then, the ionisation rate throughout the domain is fit to (2.6) to find the listed properties that determine $\hat {\gamma }(x)$ . The results from these simulations are shown in figure 3. Similar agreement is seen as for the uniform case: direct integration agrees better with the DEGAS2 than the model, especially at relatively low $\hat {\gamma }$ . Significant disagreement is seen for the ‘P1’ case at small $x$ since the ‘large’ asymptotic parameter is as low as $5\times 10^{-3}$ . Agreement is recovered deeper into the domain as $\hat {\gamma } x$ becomes larger.

Figure 3. Comparison between DEGAS2 (solid line) using constructed pedestal-like profiles versus the analytic model of (2.8) (dashed line) and direct numerical integration of (2.7) (crosses). The four colours correspond to the four pedestal shapes listed in table 2.

3.1. Uniform lossy plasma

In coordinates where $r$ is the distance to the X-point and $v$ is the corresponding radial velocity, the kinetic equation is

(3.1) \begin{equation} v \frac {\partial f}{\partial r} + \gamma f = \frac {S(v)}{2 \pi R} \frac {\delta (r)}{r}, \end{equation}

which admits the solution

(3.2) \begin{equation} f(r,v) = \frac {S(v)}{2 \pi R r v_{tn} } \frac {1}{v} e^{- \gamma r/v}. \end{equation}

With the now-familiar procedure of treating the source as Maxwellian and applying the saddle point method, the analytic model for neutral density is

(3.3) \begin{equation} n(r) \approx \frac {S_0}{\sqrt {12 \pi ^3} R v_{tn} r} \exp \left [ - 3 \left ( \frac {\gamma r}{2 v_{tn}} \right )^{2/3} \right ]\!. \end{equation}

In this geometry, the neutral density takes a slightly different form owing to the cancellation of the factor of $v$ due to the velocity integration measure $\int [\ldots ]\,v\mathrm{d}v$ . This change is not inconsequential, as the solution more clearly departs from purely exponential behaviour (Fujii et al. Reference Fujii, Hasuo, Goto and Lore2025); see figure 4 for sample cases with three different pairs of uniform electron density and temperature (as listed in table 1).

Figure 4. Comparing DEGAS2 simulation results with the closed-form model of (3.3) for three different uniform loss rates with a linear source. For visual clarity, different symbols: crosses and dots, are used to denote DEGAS2 results in the left-hand (linear scale) and right-hand (logarithmic scale) plots, respectively. Neutral density is plotted versus distance to the X-point for each volume element in the mesh. The three cases are listed as X1, X2 and X3 in table 1.

Note that the DEGAS2 results exhibit significant spread, with (3.3) capturing well the upper bound. The reason for this spread is a combination of: poloidal curvature of the lemniscate, toroidal curvature and finite extent of the numerical source. These effects are beyond the scope of the analytic model, but are worth noting and some aspects are discussed further in § 4.

Considering the spatial dependence of the neutral density from simulation, one can be easily imagine attributing the departure from exponential behaviour to charge exchange, but this reaction is not included in these cases. This is addressed in the following section.

3.2. Inclusion of charge exchange

In general, charge exchange involves the inclusion of a collision integral in the kinetic equation. However, it is often useful to approximate it as a relaxation operator (Wersal & Ricci Reference Wersal and Ricci2015; Bernard et al. Reference Bernard, Halpern, Francisquez, Mandell, Juno, Hammett, Hakim, Wilkie and Guterl2022) of the form $\gamma _{cx} \left[ (n/n_i) f_i - f \right]$ , where $f_i$ is a Maxwellian distribution of ions with density $n_i$ and thermal velocity $v_{ti}$ . Effectively, this operator acts a sink of neutrals and replaces them with those sampled from the ion distribution $f_i$ . It is exact only if ion velocities are decoupled from those of neutrals: either through a particular (unphysical) form of the cross-section or by ion velocities being asymptotically large compared with neutrals (Wilkie et al. Reference Wilkie, Laggner, Hager, Rosenthal, Ku, Churchill, Horvath, Chang and Bortolon2024). Therefore, this operator does not act on velocity space rigorously. One can imagine the approximation $v_{tn} \ll v_{ti}$ being valid, but it will not be formally applied here. The reason for not doing so is that neutrals penetrating the separatrix are likely to have already undergone a charge exchange reaction with a similar-temperature ion population outside the separatrix. If the neutral energy is very low compared with ions, they are not likely to have reached the separatrix in the first place.

With the BGK form of the charge exchange operator, the kinetic equation (with a Maxwellian source) reads

(3.4) \begin{equation} v \frac {\partial f}{\partial r} + \left (\gamma _{iz} + \gamma _{cx} \right ) f = \frac {S_0}{2 \pi ^{2} R v_{tn}^2} e^{-v^2/v_{tn}^2} \frac {\delta (r)}{r} + \gamma _{cx} n(r) \frac {e^{-v^2/v_{ti}^2}}{\pi v_{ti}^2}. \end{equation}

The ion population is taken to be spatially uniform. Equation (3.4) introduces an additional complication: a source which depends on the undetermined neutral density. This remains tractable though by solving with a similar integrating factor

(3.5) \begin{align} f(r,v) &= \frac { S_0 }{2 \pi ^2 R v_{tn}^2 } \frac {1}{rv}\exp \left [ - \frac {v^2}{v_{tn}^2} - \frac {(\gamma _{iz} + \gamma _{cx})r}{v} \right ] \nonumber \\ &\quad + \frac {\gamma _{cx}}{\pi v_{ti}^2 v} \int \limits _0^r n(r') \exp \left [ -\frac {v^2}{v_{ti}^2} - \frac {\left ( \gamma _{iz} + \gamma _{cx} \right ) }{v} \left (r - r' \right )\right ] \, \mathrm{d}r'. \end{align}

Upon integrating over velocity (with the asymptotic solution of (3.3)), the result is

(3.6) \begin{equation} n(r) = n_0(r) + \frac {\gamma _{cx}}{\pi v_{ti}^2} \int \limits _0^r n(r') \int \exp \left [ -\frac {v^2}{v_{ti}^2} - \frac {\left ( \gamma _{iz} + \gamma _{cx} \right ) }{v} \left (r - r' \right )\right ] \, \mathrm{d}v\mathrm{d}r', \end{equation}

where

(3.7) \begin{equation} n_0(r) = \frac {S_0}{\sqrt {12 \pi ^3} R v_{tn} r} \exp \left \{ - 3 \left [ \frac { \left (\gamma _{iz} + \gamma _{cx} \right ) r}{ 2 v_{tn}} \right ]^{2/3} \right \}\!, \end{equation}

is a slightly modified version of (3.3) with $\gamma = \gamma _{iz} + \gamma _{cx}$ , treating the charge exchange rate as an additional loss. The form of (3.6) is a Fredholm integral equation of the second kind. Formally, this is solved with a Neumann series (Spanier & Gelbard Reference Spanier and Gelbard2008), and this is the formal foundation upon which the Monte Carlo method is based. This is effectively a series of successive application of the integral kernel to an iterated $n(r')$ , adding additional dimensionality to the integral at each iteration. A single such iteration is applied in this work, corresponding to neutrals undergoing only one charge exchange (‘first generation’). After applying the asymptotic expansion on the velocity integral and substituting $n_0(r')$ in the integrand

(3.8) \begin{align} n(r) &\approx n_0(r) + \frac {\gamma _{cx} S_0}{\sqrt {12 \pi ^3} R v_{tn}^2 } \frac {T_n}{T_i} \left ( 1+ \frac {T_n}{T_i} \right )^{-1/6} \\ &\quad \times \left [ \frac { \left ( \gamma _{iz} + \gamma _{cx} \right ) r}{2 v_{tn}} \right ]^{-1/3} \exp \left \{ - 3 \left [ \frac { \left (\gamma _{iz} + \gamma _{cx} \right ) r}{ 2 v_{tn}} \right ]^{2/3} \left (1 + \frac {T_n}{T_i} \right )\right \}\!. \nonumber \end{align}

Once again, the model (3.8) is compared with results from DEGAS2 in figure 5 along with several other iterations for interpretation. These three cases are the same as in figure 4 (cases X1–X3 in table 1), but now with the charge exchange reaction included with ion populations at temperatures 300, 500 and 1000 eV, respectively. Several features of these cases are worthy of comment. Firstly, both the simulation results and model depart significantly from the case without charge exchange, so the effect is overall, as expected, non-trivial. Additionally, until the neutral density drops by nearly two orders of magnitude, both the model and the simulation are nearly indistinguishable from $n_0$ (3.7), which is labelled as the ‘Pure loss model’. This indicates that the effect of charge exchange near the X-point is primarily as a loss mechanism. Finally, for the most extreme case X3, there remains an order-unity departure of the simulation results from the model(s) at $r \gtrsim 1\, \mathrm{cm}$ . There are three possible reasons for this discrepancy in terms of the physics included in DEGAS2, but not the analytic model: full collision operator and cross-section for charge exchange rather than the BGK form; neutrals are permitted to undergo more than one charge exchange interaction in the simulation; and speed of neutrals with respect to the at-rest ion population may not be well captured by $v_{tn}$ .

Figure 5. DEGAS2 results for neutral density throughout the mesh plotted versus distance of the centroid each element to the X-point (pale dots), compared with the analytic model of (3.8) for the three different uniform-plasma cases. Also shown (as dotted lines) are the cases without charge exchange and (dashed lines) the prediction of (3.7) with charge exchange treated purely as a loss. Plasma parameters correspond to those listed as X1, X2 and X3 in table 1.

The fact that charge exchange behaves as a loss near the X-point is evident from direct comparison of the simulation results in figure 6. There, it is evident that the neutral penetration is decreased when charge exchange is included. This is due to the geometry of the problem. After a charge exchange event, the resulting neutral is sampled from the ion Maxwellian distribution. Most – approximately 75 % – of the post-charge exchange neutrals will therefore be on a trajectory that leaves the domain rather than penetrating deeper toward the core.

Figure 6. Examples (corresponding to case X1) of neutral density from DEGAS2 when ion charge exchange is ignored (left) or included (right).

The lossy behaviour of charge exchange observed here is not universal. In fact, the energy of the source neutrals is driven largely by the charge exchange reaction outside the separatrix, manifesting as a larger $v_{tn}$ and a lower effective loss rate $\hat {\gamma } = \gamma / v_{tn}$ .