PIC modeling of negative ion sources for fusion

This work represents the first attempt to model the full-size ITER negative ion source prototype including expansion, extraction and part of the acceleration regions keeping the resolution fine enough to resolve every single aperture of the extraction grid. The model consists of a 2.5-dimensional Particle-in-Cell/Monte Carlo Collision representation of the plane perpendicular to the filter field lines. Both the magnetic filter and electron deflection fields have been included. A negative ion current density of j H − = 500 A m − 2 produced by neutral conversion from the plasma grid is used as fixed parameter, while negative ions produced by electron dissociative attachment of vibrationally excited molecules and by ionic conversion on plasma grid are self-consistently simulated. Results show the non-ambipolar character of the transport in the expansion region driven by electron magnetic drifts in the plane perpendicular to the filter field. It induces a top-bottom asymmetry detected up to the extraction grid which in turn leads to a tilted positive ion flow hitting the plasma grid and a tilted negative ion flow emitted from the plasma grid. As a consequence, the plasma structure is not uniform around the single aperture: the meniscus assumes a form of asymmetric lobe and a deeper potential well is detected from one side of the aperture relative to the other side. Therefore, the surface-produced contribution to the negative ion extraction is not equally distributed between both the sides around the aperture but it come mainly from the lower side of the grid giving an asymmetrical current distribution in the single beamlet.


Introduction
Neutral beam injection (NBI) system based on neutralization of negative ions is essential for sustaining plasma in controlled thermonuclear fusion reactors [1]. Negative hydrogen and deuterium ions are often generated in a RF-inductively coupled plasma (RF-ICP) discharge [2] where both, volumetric by dissociative electron attachment (DEA) and surface by atomic/ionic conversion, mechanisms [3] are involved. This hybrid (volume and surface negative ion production) source mainly consists of three parts: firstly, the cylindrical driver, where RF coils are used to generate and sustain the plasma. Secondly, the rectangular expansion region, where the plasma expands into the actual source body and finally the extraction region characterised by a three-grid system: the plasma grid (PG) in contact with the source and caesiated for negative ion production, the extraction grid (EG) equipped with permanent magnets in order to separate the co-extracted electrons from the negative ion beam and the grounded grid (GG). A key role in such a discharge is played by the so called magnetic filter: a magnetic field perpendicular to the plasma flow reduces the electron temperature and density towards the multi-aperture extraction system [4,5]. This reduction has a twice beneficial effect: it decreases the co-extracted electron current from the source and it limits the negative ion electron detachment by fast electrons. Although it was agreed that the electron temperature drop through the filter is due to the enhanced electron residence time and therefore to the increased collisional energy losses of electrons there [7], the particle and energy transport through the magnetic filter is still not clear. Recent numerical models have shown the importance of the chamber walls for the electron transport across the filter and probably non-classical turbulent phenomena may be involved too [7,8]. In experiments [5], ratios of coextracted electron current to extracted negative ion current around or below 1 can be reached, strongly depending on the status of the Caesium conditioning. Typically, in deuterium much higher values of this ratio are measured than in hydrogen [5,6]. In addition, a strong inhomogeneity (top-bottom asymmetry in the plane perpendicular to the filter filed lines) is present along the plasma grid [4], as effect of electron magnetic drifts. This probably creates inhomogeneous plasma entrance conditions along the apertures (thus affecting the quality of the extracted beam [9]) and it makes the problem of negative ion extraction not reducible to the study of the single aperture in the plasma grid.
Due to the difficulties of experimental diagnostic techniques (invasive and/or non local), numerical models are becoming more and more requested. It goes without saying that an appropriate study of the negative ion source requires a self-consistent kinetic approach. The self-consistency is a quality necessary to reveal the electric field structures developing across the filter field and the fine potential map around every single extraction aperture. The kinetic description is a must for a low pressure discharge and whenever the production of reactive species (hot atoms and negative ions) is involved. In addition, such a study must necessarily couple the expansion and extraction regions, including as well part of the acceleration grid in order to reproduce the effects of the source on the extracted beam.

Critical revision of existing numerical models
State of the art concerning the self-consistent kinetic simulations of ITER-relevant negative ion source is represented by two different kinds of model.
On one hand, there are three-dimensional models of the region around an isolated single aperture of the extraction grid (see figure 1 where the cross sectional view of computational domain, boundary conditions and injection criteria are shown). One of the most important problems of these simulations is related to the arbitrariness of injection conditions at the source plane (left boundary in figure 1) and to the fact that the single aperture is considered as isolated from the others [10][11][12]. A large number of free parameters are used as input data for these models, among which: (i) initial plasma composition, density and temperature; (ii) size of initial plasma region and its distance from the plasma grid PG (shaded area in figure 1); (iii) system of re-injection of particles absorbed on the grid or extracted when two different negative ions are present (from volume and from surface).
All these parameters are arbitrary and artificial and often them are chosen on the basis of results considered acceptable and meaningful. At the same time, experimental evidences [13] and expansion region models [14] show that 5 cm from PG, plasma is not Maxwellian, the axial flux is not ambipolar (locally ¹ + -j j z z , , , where z is the direction normal to the grid) and the electron current has a strong transversal (to the grid) component [15]. The electron transport across the magnetic filter field is anomalous and non-local. Such conditions influence the physics of the plasma boundary layer and are impossible to be reproduced in a local zoomed model, where periodic boundary conditions are imposed on the lateral sides of the computational box (top and bottom lines in figure 1). Moreover, the reliability of such a models has been recently [16,17] questioned in relation to the accuracy of Poisson equation solution and to the particle re-injection method. Nevertheless, the same model proposed in [16,17] keeps some arbitrariness due to the size of initial plasma and to the periodicity on the upper and lower sides of the computational domain. In addition, the alternative suggested re-injection method is not suitable for electronegative plasmas (presence of negative ions produced in the bulk).
On the other hand, two-and three-dimensional full-size source simulations [8,14,18,19] have shown the importance of vertical plasma inhomogeneity along the extraction grid (here indicated as y-direction in figure 2) due to the presence of the magnetic filter but without the possibility of solving the transport and extraction of negative ions through the single aperture: the cell size used for the mesh is much larger than the single aperture diameter and the extraction field is not included making the sheath and meniscus around the grid hole not realistic. This is due to the fact that a scaled model with a larger vacuum permittivity (about 10 4 times the real value) is used to cover the large volume of the expansion region. In fact, such a trick increases the Debye length and reduces the plasma frequency allowing the use of larger cell sizes and time steps. Moreover, the weakness of these models is related to the fact that the density-scaling (increasing the vacuum permittivity is equivalent to reduce the plasma density) is not appropriate when turbulence and Coulomb collisions affect the electron crossfield transport.
For this reason, in the present work, a model of the full-size ITER-relevant negative ion source has been developed keeping the grid cell small enough to resolve in detail the extraction dynamics of every apertures till the EG grid. The details of the model will be presented in the next section, while results are discussed in section 4.

Description and assumptions of numerical model
The model (whose computational domain is reported in blue in figure 2(a) consists of a 2.5-dimensional Particle-in-Cell/Monte Carlo Collision (PIC-MCC) technique [20] simulating the region going from the driver exit plane (located at = = z z 16 cm 0 and between 17 cm < < y 41 cm) to the entrance of the EG (at = z z EG ) of the ITER ion source prototype BATMAN [2].
The simulation starts from scratch with no particles in the computational domain. The driver is not selfconsistently simulated, but a prescribed production rate of electrons, H + and + H 2 ions is imposed at the driver exit plane (red area named 'Injection region' in figure 2(a)). Following [21], the total ionization rate is set to = --· w 3 10 m s ion 23 3 1 , whose 56% concerns + H 2 via direct ionization of H 2 , while the remaining 44% concerns the H + ion production by two channels: dissociative ionization of H 2 and direct ionization of H. To ensure that the injected particles keep the prescribed fixed temperature (set to 12 eV for electrons and to 1 eV for ions), a continuous Maxwellianization of particles present in the injection area is processed. A fixed atomic and molecular density, = -· n 0.8 10 m  The z (perpendicular to PG) and y coordinates are simulated while uniformity is considered along the direction x of the magnetic filter field (B x is directed inward the yz plane). This choice is justified by the interest of studying the E×B and diamagnetic Ṕ Belectron drifts [25] lying in the (yz) plane: (here k B and e are the Boltzmann constant and the elementary charge, respectively). The use of the nomenclature 2.5-dimensional (2.5D) comes from the fact that although the self-consistent electric field is not solved along the x-direction, the particles are tracked along it and particle-wall interaction with a thin-sheath approximation is considered at the x-boundaries as follows: if an electron reaches the lateral wall with an energy larger than the local electric potential, , it is lost and a secondary electron is applied; negative ions reaching the walls are lost if their energy is larger than the local potential, otherwise they are reflected back; finally, all positive ions reaching the walls are lost. Particles produced by ionization, as well as those injected into the injection region, are all created anywhere between   x L 0 x (where L x =32 cm) and all the two-dimensional maps in the plane (yz) are averaged along the x-direction.
Based on the 2.5D representation, the driver exit is not circular but a rectangle of area is the driver diameter and = L 32 cm x is the extension along x of the expansion region), while the circular apertures in the PG become slits of length x h =L x and aperture = D 8 mm h (see figure 2(b)). In the present model, ten flat apertures are considered. PG thickness between each aperture is The distance between PG and EG is 3.6 mm. The choice of extraction grid parameters in the model has been done in order to: (i) preserve the same production area A prod  530 cm 2 (black region in figure 2(b)) of the so-called Large Extraction Grid (LAG) used in BATMAN [27]; this means that the total current of negative ion produced on the surface -I H s , is equal to the experiment and therefore, a correct current on PG = -- , is guaranteed for the continuity equation; (ii) keep the size of the apertures and the gap between apertures similar to that of LAG; this guarantees a true representation of the ion dynamics and extraction.
Nevertheless, due to 2.5D representation, this choice leads to a larger extraction area (red region in figure 2(b)) by a factor 4 (A extr =256 cm 2 against 63 cm 2 of LAG [27]). For this reason, in order to compare the model with experiment, we have decided to present results in term of extracted current density instead of total extracted current.
The filter field assumes a bell-shaped z-profile with a peak value of = B 7 x,max mT at = z 3 cm max from PG location (so called standard configuration): (no dependences along x and y are included and the possible mirroring effects due to the filter field curvature are neglected). The electron deflection field (y and z components) in each aperture (reported in figure 1) is assumed to have a cusp-like behaviour: 2 cm and y is referred to the center of the aperture. The B y component sign is alternated from one aperture to the other.
Negative ions are self-consistently produced by DEA to vibrationally excited molecules and by ionic conversion on caesiated PG. The conversion yield is taken from [28]: where E is the positive ion impact energy per atomic mass, Here, R N = 0.8 and R E =0.7 are the particle and energy reflection coefficients, respectively, while h 0 is the electron transfer probability and E th is the threshold energy. In addition, a fixed negative ion flux emitted from the PG surface In order to avoid that the initial condition influences steady state results, the negative ion production by neutral conversion is activated only when electron extraction takes place, after about 10 μs from the beginning of the simulation.
The bias plate and PG (in green and black, respectively, in figure 2) bias are set to a potential f f = = 0 BP PG (usually, to reduce the co-extraced electron current, PG is set to plasma potential condition at about f = 23 V PG [6]), while the EG (in red in figure 2) potential is f = 1 kV (s is the sheath thickness). Fitting the experimental values [27] of the extracted electron current density j e as a function of the extraction potential for the LAG system, gives an exponent of a = 1.2 in place of the classical value 1.5. All the other side walls are grounded and Poisson equation is solved using Portable, Extensible Toolkit for Scientific Computation (PETSc) software [29].
Different bulk collisions are implemented using Test Particle and Direct Simulation Monte Carlo Collision techniques [30]. The list of relevant collisions included are reported in [31].
Differently from previous full size source models [8,14,18,19], a more realistic value of vacuum permittivity has been used (e e ¢ = 25 (the index j runs over the different charged species) detected inside the source. The site corresponding to l D,min is located inside the potential well attached to the PG surface (see figure 9 where the behaviour in the last 5 mm ending to the PG surface is reported) and it is mainly due to the negative ion contribution to equation (7). The code uses a hybrid OpenMP/MPI parallelism paradigm and it spans about 2 μs/day using 8 Intel®Xeon E5-2650Lv2 10-Core CPU.

Results and discussion
In this section, results of the 2.5D PIC-MCC model at steady state are presented. In order to reduce the numerical noise, all the global integral quantities (such as current, energy, etc.), one-dimensional profiles and two-dimensional maps have been obtained averaging over the last 1000 PIC cycles. We first analyse the electron transport across the magnetic filter in the expansion region and then the electron and negative ion extraction through the multiaperutre grid.

Non-ambipolar transport across the magnetic filter field
In figure 3 electron G e and molecular positive ion (the dominant ionic species) G + H 2 fluxes with streamlines are displayed in the plane (y, z).
The electron flow exhibits a quite complex distribution. It is first evident a strong electron flow directed towards the bottom surface at the entrance of the magnetic filter region = z 28 cm. It is the effect of electron density and temperature gradients (see also figure 5) that induces a diamagnetic drift in the negative y-direction (see equation (1) with a negative pressure gradient). This electron flow interferes with the electron flow coming from the driver bending it in anticlockwise direction in the upper part and in clockwise direction in the lower part. As a result of this interaction, an electron backflow towards the driver sidewall and a short-circuit close to the bottom surface are established. Only some electrons (the low-energy one after many collisions) are able to escape and cross the filter region towards the PG. Due to the strong filter action, the streamlines show a transition from classical to anomalous regime at = z 31 cm, i.e. about 8 cm from the extraction grid. From here to the extraction grid a more chaotic transport characterised by the presence of eddies and magnetic mirroring along the electron deflection magnetic field lines is established. It is also worth noting that, as a result of the strong electron flow along y, the positive ions flow close to PG has a marked vertical deviation towards the upper part (see figure 3(b)). Enlarged views of the extraction region and effects of this asymmetry will be discussed in the next paragraph. Comparing the electron and positive ion flow, the most important macroscopic effect created by the magnetic filter field is the non ambipolar character of the transport inside the expansion chamber. Although the total diffusion must be ambipolar (it is not the fluxes, but the divergence of the fluxes that must be set equal to each other,  G =  G · · e i ), the one-dimensional part of the losses need not be ambipolar: ions diffuse almost uniformly inside the source, while electrons diffuse primary along x and towards the expansion region backplate at = z z 0 . Since G G  e z i z , , , in a simple 1D(z) case it would be expect an electric field E z to be established which accelerates electrons and retards ions. However, in 3D this electric field is short-circuited by an imbalance of the fluxes along the direction of the magnetic filter field x. Here, a self-consistent solution of the field along x is missed (uniformity is assumed) but results show that electron flux along x is much larger than the ionic counterpart G G A. The almost totality of positive ion current is equally shared between z 0 , z PG and at the walls normal to x-direction, with the + H 2 contributing with one third more respect to H + ion. In figure 4, electric potential, electron temperature and density maps are reported, respectively, while in figure 5 axial profiles of the three different quantities have been extracted from these 2D maps at three different y locations, = y 19 cm (bottom ), = y 29 cm (center +) and = y 39 cm (top ), visualized in figure 4(a). Following the central line, the electric potential decreases almost linearly in the expansion region from the driver exit (where a value of 36 V is detected) to about 25 V in the filter region. Here the typical oscillations of a plasmawall transition in an oblique magnetic field [32] are present. A similar behaviour is also observed in the electron temperature and density maps: after a peak of 10 eV and -· 5 10 m 17 3 at = z 16 cm, respectively, both decrease to 3 eV and -· 3 10 m 17 3 at the entrance of the extraction region. This clearly shows the working principle of magnetic filter field [3]: making the electron transit time larger and the flux almost entirely parallel to the magnetic field lines in the filter region, electrons reduce their energy by collisional events and move faster along transversal direction. These effects reduce both the electron temperature and density downstream of the magnetic filter region. All the quantities show a marked top-bottom asymmetry due to the aforementioned magnetic drifts. The asymmetry becomes more pronounced at = z 30 cm in the filter field region: the electric potential takes values 5 V larger than those of the bottom part. Note that the electron temperature is lower where the electron density is larger (top region) and vice-versa, in agreement with experimental measurements [4] and results from different models [8]. This tends to make the electron pressure more symmetric. It is also interesting to see in the electric potential map, fine microstructures on the top part in the magnetic filter region. The electric potential is modulated: the amplitude of the oscillations reaches 10%, while the wavelength is about 1 cm. Such structures have already been detected by 1D(z) model [7] where they were identified as electron-ion hybrid (EIH) mode due to the presence of a magnetic field gradient. Similar structures were also found in Hall-effect thrusters (HETs) [33] (where the E × B drift is closed along the azimuthal direction of a cylindrical coaxial channel; no walls are present along the direction perpendicular to the magnetic field and to electron flow). In HETs, it is commonly accepted to identify as the origin of azimuthal modulation, the electron cyclotron drift instability (ECDI). Here, a dedicated study is necessary to estimate the contribution of this modulation to the electron transport along z compared to the that induced by the ambipolar field E y .

Extraction through the multiaperture grid
Despite the inhomogeneous plasma condition induced by the magnetic filter field, the extracted current show a reasonable uniform distribution along the different apertures ( figure 6). Electron current shows a deviation of about 30% between the bottom and top apertures. In particular, except the two extreme apertures, all the others   figure 1), the H − flow direction is opposite to that of positive ions (shown in figure 3(b)), a fact that has recently found experimental evidences [15]. As a consequence, the  is presented as asymmetrical lobes with a penetration length inside the source increasing from the bottom to the top aperture. Even the potential well (result of the unbalanced negative charge produced on the surface) attached to the PG surface is not uniform between two adjacent apertures, but it reaches a maximum depth of 1.5 eV from one side, the upper one, of the aperture and it is absent on the other side, the lower one (see also figure 10(a) where the y-profile of the electric potential one cell from PG is reported for the same central aperture). The location of absence of the potential well corresponds to the place where a larger positive ion flow arrives [34]. Therefore, surface-produced negative ions have preferential trajectories always directed towards the bottom part of the source (see both figures 7(c) and 8(c)) and the extraction occurs where the negative ion streamlines intercept the meniscus. The remaining part of the surfaced-emitted H − flux is lost inside the source.
In order to demonstrate the good spatial resolution of the simulation, the z-profiles of electric potential and plasma density in the last 5 mm along a line in the middle of the y-coordinate (ending at the PG surface) is reported in figure 9. As it can be seen the plasma bulk potential takes a value of 25 V, while electron density at the entrance is   The asymmetry detected around every single aperture has an impact on the quality of the single beamlet extracted. In figure 10(b) the y-profile of the negative ion current density for the central aperture is reported as computed at EG plane. It shows the presence of an important tail in the lower-part of the negative ion current density profile due to the fact that the majority of negative ions are extracted from the lower-side of the aperture. Here, the first half of PG surface has no potential well allowing surface-produced negative ions to move into bulk up to the meniscus line where they reverse the motion towards the aperture. On the contrary, the contribution coming from the upper-side of the aperture is strongly reduced by the presence of a potential well laying for the first half of PG surface. The area of the tail increases going from the lower to the upper aperture.
Among the different surface contributions (the volume one is negligible) to the negative ion production, about 54% comes from the positive ion conversion, while the remaining 46% from atomic conversion. This is due to the fact that 12% of the total atomic flux coming from positive ions reaching PG (G + G + + 2 H H 2 ) is converted into negative ion flux. Such a high conversion yield (in comparison the conversion yield of atoms is only 2%) is due to the energy distribution function of ions hitting PG: a large fraction of H + and + H 2 has a kinetic energy larger than the conversion threshold (for + H 2 the threshold is twice in order to allow the conversion of each atoms of the molecular ion). A surface-produced negative ion has to overcome first the potential well and then reverse its motion towards the aperture. This makes the extraction probability reaching a value of 10% in agreement with a previous Monte Carlo transport model [35]. The most important mechanism leading to motion reversal towards the apertures is the electrostatic penetration of the EG field inside the source region, i.e. the meniscus interception of the surface-produced negative ion flow which occurs up to a distance of 5 mm from the aperture (see the bending of negative ion streamlines in figure 8(c)). However, it should be noted that the EG field penetration depth is larger than in the real case [36] due to 2.5D representation (a penetration through a slit is deeper than a penetration through a circular aperture) and to the fact that the extraction area is 4 times bigger. On the other side, collisions play only a minor role to the motion reversal: when extracted only 30% of the surface-produced negative ions have suffered at least one collision with a neutral particle.

Conclusions
A first attempt to kinetically and self-consistently model the whole expansion region of the ITER-prototype negative ion source has been done, keeping the mesh size small enough to also resolve the extraction and part of the acceleration regions of the multiaperture grid. Results show the strong non-ambipolar character of the transport inside the source: diamagnetic drift contributes to the creation of a top-bottom asymmetry which also affects the positive ion flow hitting the plasma grid. As a consequence it induces a top-bottom asymmetry around every single aperture: the potential well is not uniformly distributed at the lower and upper sides of the aperture creating a preferential direction of surface-produced negative ion flow from the PG surface. This has important effects on the quality of the beam optics. Future investigations will be devoted to study the effects of the  0 -scaling, of magnetic filter field, PG bias voltage and aperture size and shape. Figure 10. y-profiles of (a) electric potential f one cell from PG surface and (b) negative ion current densityj H as computed at EG plane, for the central aperture.