Integration of kinetic ions in a three‐dimensional Monte‐Carlo neutral transport code

Charged particles of species different from the major constituents of a plasma might not reach local thermodynamic equilibrium during their lifespan. Such low collisional ions render the often‐used fluid description, assuming Maxwellian velocity distribution, as insufficient. Hence, kinetic treatment of such minorities becomes inevitable. The three‐dimensional kinetic neutral particle Monte‐Carlo code EIRENE [D. Reiter et al., Fusion Science and Technology, 47 (2), 172 (2005)] possesses a model for ion transport reduced in physics, consisting of field‐line tracing and energy relaxation of particles. The trajectory integration is now extended to account for first‐order particle drifts, anomalous cross‐field diffusion denoting for turbulence effects, and the consideration of the magnetic mirror force. The technical implementation, as well as the verification of the enhanced kinetic ion transport is being reported.


INTRODUCTION
Nuclear fusion via magnetic confinement is a research topic requiring experiments very expensive in construction and execution. The vast amount of non-linearities in such hot, dense, and turbulent plasmas makes interpretation of experimental data difficult. Numerical simulation possesses the power to validate physical models by comparison to measurements. Thus, modelling opens the possibility to both, predicting future fusion campaigns and assessing physical insight of experiments.
Particles in the scrape-off layer might hit wall or target plates and neutralize by binding a surface electron. Therefore, for a complete physical picture in the plasma edge, a combined neutral-plasma code package is required. The European standard code for kinetic neutrals is EIRENE, [1,2] mainly for three arguments. First, it is coupled to extensive atomic and molecular databases, enabling a multitude of different plasma chemical reactions and physical collisions, respectively. Second, the code is applied to realistic geometries, as, for example, the ITER divertor design has been conceived using EIRENE. [3] Third, it is coupled to all major European fluid plasma edge codes, [4][5][6][7] allowing for such sought-after combined neutral-plasma simulations in the edge of nuclear fusion devices.
While in tokamak studies there is often assumed toroidal symmetry, one uses only two generalized coordinates in such simulations. However, intrinsic three-dimensional effects such as local gas pumps or puffs, the addition of external error field perturbation, and the consideration of the sophisticated stellarator geometries require a full three-dimensional treatment. The edge Monte-Carlo plasma fluid code EMC3 [6] does meet these requirements. Additionally, its algorithmic structure via Markov chains renders the code a well aligned coupling partner to EIRENE: While in the latter, one follows single particles undergoing chemical and physical transitions depending on the local background conditions, the former does track plasma parcels, that is, packages of energy and matter, from which aforementioned background gets calculated.
For charged particles that are not of the main plasma species, the assumption of local thermodynamic equilibrium might be wrong due to a short lifetime life compared with the collision frequency coll , where coll life ≪ 1. This renders a fluid description, in which the ion velocity distribution is assumed to be Maxwellian, as insufficient. It has been shown that, for example, the lower charge state impurities of carbon C + -C 3+ [8] and the minorities of helium He + [9] have to be treated kinetically.
EIRENE holds the possibility of describing ion transport, albeit with reduced physical fidelity, namely field-line tracing and energy relaxation only. Early approaches of including more effects for application in the two-dimensional domain [10] have been extended to full 3D geometry. [11,12] In this work, we introduce the technical details of the introduction of first-order drift effects, anomalous cross-field diffusion, and magnetic mirror force into the kinetic ion transport of EIRENE in the EMC3-EIRENE environment.
This paper is organized as follows. In Section 2, we introduce the underlying formulae, while Section 3 hosts the explanation of the technical implementations into the full unstructured EMC3 grid. In Section 4, we evaluate and verify the introduced algorithms and wrap up, and we draw conclusions for a forthcoming publication in Section 5.

PHYSICAL MODEL
The equation of motion for kinetic ions in EIRENE is of Fokker-Planck type where f is the particle distribution function, / t is the partial derivative with respect to time t, ∇ is the Nabla-operator, v is the particle velocity, and D ⊥ is a perpendicular diffusion coefficient in the units m 2 /s denoting for turbulence effects. A solution of Equation (1) is obtained by a Monte-Carlo approach, following the guiding centre orbits Δx gc of particles given by where v gc is the guiding centre velocity, Δt is an incremental time duration, and̂is a normalized random vector |̂|= 1 perpendicular to the magnetic field̂⋅ B = 0, where B is the magnetic flux density in vacuum ∇ × B = 0. In the field-line tracing approach, the guiding centre velocity is given by the component parallel to the magnetic field v gc = v ‖ and the diffusion is set to zero D ⊥ = 0. Note that in case of nonvanishing D ⊥ , the diffusion term in Equation (2) introduces position kicks, which do not change the particle's kinetic energy.
A more complex physical model produces deviations from this field-line tracing approach. We introduce where E is the electric field, B is the absolute value of B, v ‖ and v ⊥ are the absolute velocities in parallel ‖ and perpendicular ⊥ direction to B, respectively, q j and m j are the charge and mass, respectively, of a particle of species j. Equations (3) and (4) are the standard expressions [13,14] for first-order electromagnetic particle drifts v drift , consisting of E × B-, ∇B-, and curvature-drift, and the magnetic mirror force v mirror being in charge of the conservation of the magnetic moment = m j v ⊥ 2 /2B. Note that the total guiding centre velocity is now given by Equations (2) and (5) imply that drifts and magnetic mirror force do change the guiding centre velocity v gc , while the cross-field diffusion causes discrete jumps in random directionŝ.
Implementation of Equation (3) is straightforward, while Equation (4) is an ordinary differential equation in time, which must be integrated. We implemented both, a Euler-forward and a fourth-order Runge-Kutta scheme, the latter giving results more accurate (Section 4).
Note that drift motion is purely transversal to the parallel motion v drift ⋅v ‖ = 0, while the mirror force acts purely parallel v mirror × v ‖ = 0. In the vicinity of a target, a gyro angle gets sampled and the integration is performed using the fully resolved particle motion, including its gyration.

TRAJECTORY INTEGRATION IN UNSTRUCTURED GRIDS
In EMC3-EIRENE, grid information is stored on the EMC3 side, where cells are called bricks. While EIRENE uses Cartesian coordinates (x, y, z), EMC3 uses cylindrical coordinates (r, , z), and, additionally for charged particles, local cell coordinates (l, m). A computational brick in EMC3 consists of radial, toroidal, and poloidal surfaces, which are again split into two typically non-parallel triangles each. Toroidal surfaces share a fixed azimuthal angle and two toroidally consecutive brick nodes lie on one single magnetic field line, as can be seen in blue in Figure 1. Trajectory integration in such unstructured grids is different for neutral and charged particles, as for the former it is straightforward. For the latter, if transverse effects are disregarded, the next intersecting surface can only be toroidal. If an ion does not perform any collisional events during its time of flight in one brick, its local coordinates (l, m) on one toroidal face get mapped to the next one. For EIRENE, which intrinsically is unable to host curved trajectories, the two intersection points produce a virtual and straight path v virt which is followed, resulting in a small mistake by underestimating the travelled distance. If the particle does perform a collision, above's procedure stays the same, but for the fact that a virtual toroidal surface is set up in the brick to which the local coordinates (l, m) get mapped. By the addition of drifts, diffusion, and the mirror force, this trajectory integration is now changed: The above-mentioned virtual velocity v virt is altered in its components transverse to the magnetic field by drifts and diffusion, resulting in a new virtual velocity v virt ′ . As now it is also possible that the trajectory intersects a radial or poloidal surface, the location on the initial toroidal surface (l, m) can no longer be simply mapped. With v virt ′ , one enters the EMC3 integration scheme for neutral particles, which returns the correct next intersection point and surface. After the projection to the next surface is done, the change to the parallel velocity v ‖ (and, because of energy conservation, to the perpendicular velocity v ⊥ ) due to the mirror force is calculated. This change might result in a sign-flip, causing the particle to turn. Obviously, this introduces another error, as the sign-flip might occur before the particle is projected to the next intersection. However, this error is on the same order of accuracy as the whole trajectory integration in unstructured bricks and, thus, negligible.
Equations (3) and (4) require the gradient of the magnetic field strength ∇B at an arbitrary position inside a brick. An optimized scheme, harnessing the B-field accuracy of EMC3 bricks, has been implemented. There is one absolute value B stored at each brick node, the direction B/B implicitly given in the grid construction. For interpolation, a bilinear mapping is applied, namely a calculation of B at the projected position (l, m) at each toroidal face in combination with a linear interpolation between the two angles t at one surface and t + 1 at the consecutively next one. This way, the gradient of the magnetic field strength can be obtained much faster than in the previously used Cartesian stencil [11] approach, which is important if the code enhancements introduced in this paper shall be optimized and, thus, be relevant to the user community of EMC3-EIRENE.
Note that Feng et al. are currently working on the self-consistent calculation of the electric field. [15] Therefore, the here presented enhancements can only host either a hard-coded or externally read-in electric field, which is necessary to calculate the first term on the right-hand-side of Equation (3), denoting the E × B-drift contribution. As this term is not self-consistently calculated, it is artifically omitted in the following investigations. Furthermore, we remark that parallel variations of the equilibrium potential might be considerable in the tokamak edge region. Hence, the electric field will exhibit a parallel compontent E ‖ which must be considered in the parallel motion however, it is neglected for the time being as no self-consistent calculation is possible.

VERIFICATION
Verification of new code development is an indispensable aspect of proper enhancement of the included physics. For the case study presented here, we use an envisaged ITER scenario from the baseline 2008 database, Ψ N = 0.83, case 2297. [16] This is a deuterium plasma case in which we artificially seed nitrogen for verifying our enhancements, using a fourth-order Runge-Kutta integration scheme for solving Equation (4). Evaluating multiple single particle trajectories tracing along the magnetic field while collisions were turned off, we analysed the the transverse displacement L ⊥ in the projected R-z plane with respect to the travelled distance in longitudinal direction L ‖ . From this transverse displacement L ⊥ , which is purely due to numerically erroneous integration, we can deduce a numerical diffusion coefficient D num ∼ L ⊥ 2 /2 t, where t = L ‖ /v gc . For t = 1 s we found a parallel length of L ‖ ≈ 200 km and a perpendicular displacement of L ⊥ ≈ 3.5 cm, resulting in a numerical diffusion coefficient of D num ≈ 10 −5 m 2 /s. Compared with typically used quanta for transverse diffusion coefficients of D ⊥ ≈ 1 m 2 /s, the numerical diffusion is negligibly small D num ≪ D ⊥ , and, hence, the integration scheme introduced in Section 3 is positively evaluated.
A more advanced verification scheme uses the differentiation of particle orbits into passing and trapped ones. For such banana orbits, as the latter are called, there exists analytic theory valid in spherical plasmas, which are axisymmetric and where B is given analytically. By analysing such orbits, one can evaluate the interplay between ∇Band curvature-drift, as well as mirror force, which are necessary for observing banana orbits. [13] The critical value for deciding whether a particle is trapped or passing is the pitch angle , which is given by the square root of the ratio of energies = √ E ‖ ∕E ⊥ that is stored in the parallel motion and in the gyration, respectively. One finds, for spherical plasmas and an axisymmetric magnetic field, where crit is the critical pitch angle above which an orbit is supposed to be passing, r is the distance of a particle orbit from the centre of the plasma, R 0 is the major radius, d banana is the diameter of a banana of particle j, and B is the magnetic field value in polar direction, when (r, , ) is the set of spherical coordinates. For ITER we find R 0 = 6.2 m and insert non-colliding singly charged nitrogen ions m j = 14 u, where u ≈ 1.66 × 10 −27 kg is the atomic mass unit, q j = 1 e, where e ≈ 1.6 × 10 −19 C is the elementary charge with v ‖ ≈ 1.9 × 10 5 m/s at the position r ≈ 1.9 m, where we find B ≈ 1.2 T. These values give us a critical pitch angle of crit ≈ 0.25. As applying this theory Equations (6) and (7) to ITER geometry is already overextending the model boundaries, we choose two pitch angles 1, 2 , which are well above and below, respectively, the critical one. F I G U R E 2 Trapped ion orbit with the pitch angle below critical 1 < crit F I G U R E 3 Passing ion orbit with the pitch angle above critical 2 > crit Figures 2 and 3 show the two-dimensional projection of two particle orbits, one with 1 = 0.21 (in Figure 2), and one with 2 = 0.3 (in Figure 3). The injection point is marked with a red star. We do find indeed the pitch angle

SUMMARY AND OUTLOOK
We introduced some very fundamental physical enhancements of the kinetic ion transport part of EIRENE, namely by adding first-order drift effects, cross-field diffusion, and magnetic mirror force. These additions, which are relevant for thoroughly investigating the full three-dimensional influence of impurities on actual fusion devices, have been verified by analysing passing and trapped particle orbits, as well as checking on the introduction of numerical diffusion by our integration scheme. A full three-dimensional ITER case study, where the influence of nitrogen impurity seeding on the deuterium background plasma is being incorporated using this enhanced EIRENE kinetic ion transport, is subject of a forthcoming publication.
For the fluid code community, it might be interesting to model ion species kinetically and track their collisionality. This way, justification for the physically always less accurate description assuming local thermodynamic equilibrium might be found.
Application of the here-presented enhancements to heavier impurity species such as tungsten requires a revision of the Fokker-Planck collision operator on which such Coulomb collisions of heavier constituents rely. [17] ACKNOWLEDGMENTS The author thanks Michael Rack for proofreading and constant support. The author thanks Marco Wischmeier for fruitful discussions. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. This work was supported by the EUROfusion-Theory and Advanced Simulation Coordination (E-TASC) and has received funding from the Euratom research and training programme 2019-2020 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. The authors gratefully acknowledge the computing time granted through JARA-HPC on the supercomputer JURECA at Forschungszentrum Jülich.