Simulations of negative hydrogen ion sources

The development and the optimisation of negative hydrogen/deuterium ion sources goes hand in hand with modelling. In this paper a brief introduction on the physics and types of different sources, and on the Kinetic and Fluid theories for plasma description is made. Examples of some recent models are considered whereas the main emphasis is on the model behind the concept and design of a matrix source of negative hydrogen ions. At the Institute for Nuclear Research and Nuclear Energy of the Bulgarian Academy of Sciences a new cyclotron center is under construction which opens new opportunities for research. One of them is the development of plasma sources for additional proton beam acceleration. We have applied the modelling technique implemented in the aforementioned model of the matrix source to a microwave plasma source exemplifying a plasma filled array of cavities made of a dielectric material with high permittivity. Preliminary results for the distribution of the plasma parameters and the φ component of the electric field in the plasma are obtained.


Introduction
The negative hydrogen/deuterium (H − /D − ) ion sources are an integral part of devices as proton accelerators [1] and magnetic confinement fusion reactors of the tokamak type (ITER -currently under construction at Cadarache France). In the field of fusion applications the Negative Ion (NI) sources are a vital component of the systems for additional heating -Neutral Beam Injection (NBI) -of tokamak fusion plasmas, as they complement the ohmic one. Development of these systems for the ITER tokamak is an extensive task allocated in fusion research institutes worldwide [2]. The H − ion sources found their first application in the Van der Graaf accelerators, which dates back to the 50's of the 20-th century, and later on in the cyclotron accelerators. The currently used and developed H − ion sources in proton accelerators are of the same type as the ones in the systems for additional heating of thermonuclear plasma. Regarding the discharge type the sources are: DC (direct current, based on glow discharge or arc discharge), Radio Frequency (RF) and microwave discharges.
There are two mechanisms for production of H − ions -volume (in the plasma volume) and surface production (on the discharge walls). The H − ions are created in the plasma volume through dissociative attachment of electrons to vibrationally excited molecules (e + H 2 (ν) → H + H − ). This is a two step process: excitations of the vibrational levels (ν) of the ground electron state through collisions with fast electrons (energy above 10 eV), followed by dissociative attachment of slow electrons (energy 1 eV) to vibrationally excited molecules. The effectiveness of H − surface production (H, H + + e wall → H) depends on the surface material. Since a low work function material is required, usually Cs or Ba are used.
The two main types of H − ion sources are DC and RF. In both of them the volume production of H − could be optimised by employing a magnetic filter (MF) (Figure 1 and 2 magnetic cusps with fields that effectively separate the volume of the source in one region with fast electrons and another region with slow electrons, i.e. it reduces the electron temperature. NI sources with MF are called tandem sources. The plasma in the DC source is ignited by a hot filament (Figure 1), whereas the typical RF sources are inductive discharges and the plasma is generated by the RF field introduced in the driver region by RF coil (Figure 2). The other two regions in the RF NI sources are, respectively, expansion (the plasma generated in the driver region diffuses into it) and extraction regions. The NI are extracted via an extraction system which consists typically of 3 multi-aperture grids: the plasma grid (PG), the extraction grid (EG), the grounded grid (GG). The NI are extracted from the plasma by applying a positive voltage (5-10 kV) between the PG and the EG. Subsequently, electrons are accompanying the NI extraction, hence they need to be deflected. For this purpose secondary a magnetic field, created by permanent magnets embedded in the EG (Figure 2), deflects the electrons out of the beam on the grid components.  Depending on their application, NI sources vary in design and beam parameters. For example, the beam requirements for the full size ITER neutral beam injection experiment Source for Production of Ions of Deuterium Extracted from Radio-frequency plasma (SPIDER) are: 40 A beam current; 3600 s pulse length; 1 MeV beam energy [2]. NI source for the spallation neutron source project at the Oak Ridge National Laboratory has the following characteristics: RF source with surface production; 40 mA beam current; 1 ms pulse length; 2.5 MeV beam energy [3].
In the research and development of H − ion sources the key issues are summarised as follows: improvement of the H − ions production efficiency under low pressure conditions; development of large volume H − ion sources with spatial uniformity; optimisation of the H − ions extraction from the source; high energy acceleration and good beam optics. The rapid development of computer resources has boosted the computational methods in plasma physics. In the next part the two theories for plasma description (respectively, the kinetic and fluid theories) shall be presented and, respectively, simulation models in which the Particle In Cell (PIC; for collisional models the Monte Carlo Collision (MCC) scheme is applied [4]; we shall use the PIC-MCC notation for such models) or Finite Element Method (FEM) are implemented. Both of them are well established tools in the field of plasma physics simulations, more information could be found in the literature [4][5][6]. These models will show part of the progress that has been made in the field of simulations of H − ion sources. The main topic will be the H − matrix ion source which concept [7] is built upon self-consistent model coupling gas discharge part within the fluid plasma theory with electrodynamics. The modelling equations are numerically solved by applying FEM. In the current study an example of the aforementioned modelling technique will be presented as we applied it in the development of a two-dimensional, axially symmetric, self-consistent model of a new plasma source. The source exemplifies a plasma filled linear array of cavities made of material with high dielectric permittivity. The discharge is in argon. The gas pressure is 10 Torr. We have obtained preliminary results for the distribution of the electron density and temperature; phi (ϕ) component of the electric field. The plasma is generated by a transverse electric (TE) wave with a frequency of 2.4 GHz, and 160 W deposited power. The application of this source is in the field of a plasma source for additional acceleration of a proton beam. The development of such a plasma source could be a part of the research program of the cyclotron center which is under construction at the Institute for Nuclear Research and Nuclear Energy of the Bulgarian Academy of Sciences [8].

Plasma description theories, H − ions source modelling
There are two theories for plasma description, respectively -kinetic and fluid. The kinetic theory is based on obtaining the distribution function f (r, v, t), which gives us the distribution of the particles in the six dimensional phase space (position r and velocity v space of the particles) at a given moment of time t. The distribution function has to satisfy the Boltzmann equation: Where F is the force acting on the particles and ∂/∂v -is the velocity space gradient. The right hand side term of the equation is the collisional term. In a sufficiently hot plasma the collisions can be neglected and furthermore the force is usually entirely electromagnetic (Lorentz force: F = q(E + v × B); q -electrical charge of the particles; E-electric field; Bmagnetic field), then the Boltzmann equation takes it's special form -the Vlassov equation.
Having obtained the distribution function by solving the Boltzmann or respectively the Vlassov equation we can obtain the macroscopic parameters. Typically we need the distributions of the following (macroscopic) parameters: the density of charged particles of certain species n α (r, t); their velocity v α (r, t) and their temperature T α (r, t) (α-electrons or ions). The kinetic theory is applied in the cases when non-Maxwellian plasmas and non-local effects are studied like the electron plasma oscillations and the Landau damping (particle-wave interactions).
Fortunately in many cases the plasma species have Maxwellian velocity ditribution. Thus a more crude model (than the kinetic theory) for description of plasma phenomena is applicablethe fluid theory for plasma description. In the fluid description typically the first three moment equations are included, respectively: the continuity equation (2); the fluid equation of motion (momentum transport equation) (3); the electron energy balance equation (4): Where, in these equations n α , m α , T α and p α are, respectively, the density, mass, temperature and partial pressure for given type of species. u α is the fluid velocity, µ αβ = m α m β /(m α + m β ) is the reduced mass for elastic collisions between particle species α and β, Z α = ±1 respectively for ions and electrons, Γ e is the electron flux, J e = −(5/2)n e D e ∇T e + (5/2)n e T e u e -electron energy flux, consisting of convective and heat fluxes (D e -electron diffusion coefficient); Q ext = (1/2)(Reσ)|E| 2 is the deposited power by the external RF field, the δn α /δt term is the collisional generation and loss of charged particles through the processes of ionization and recombination. On the right hand side of (4) despite the collisional losses P coll the losses for generation of constant electric field are also noted (−eΓ e .E DC , where E DC is the constant electric field inside the discharge).
Since the two theories for plasma description are being discussed we have to note that the plasma is a dynamic system and respectively its charged particles are in constant motion. The charged particle fluxes induce currents inside the plasma volume, which interact with each other. Those interactions on their turn cause modifications to the macroscopic plasma properties. Thus in order to have a self-consistent decription of the plasma we need to include in our set of equations in a proper manner the Maxwell equations. Depending on the geometry and the manner of power deposition in the discharge the Maxwell equations could be reduced, respectively, to: Poisson's equation (simulation of DC discharges; extraction region of the NI sources [9,10]); equations for the azimuthal components of the electric or magnetic field (simulation of cylindrical RF inductive discharges), other alternative is the vector potential equation of the RF magnetic field; electromagnetic wave equations (in the case of surface microwave heated discharges).
The NI sources are sophisticated devices, thus a number of issues should be taken into account in the research and development process. Here we shall consider issues related to the extraction region and the RF power deposition in the NI sources. RF power deposition models are the tools which give us insight on the physics of the different states and operating regimes of the discharge, and thus help us to optimize their design. The physics of the NI extraction region requires detailed modelling by itself, since many processes and parameters should be taken into account in order to insure the beam quality.
An example of electromagnetic PIC-MCC models for studying the RF power deposition are the models employed for the plasma modelling of CERN Linac4 H − and The J-Parc Linac sources [11,12]. The Linac4 H − ion source model was further developed [13,14] giving insight on the effect of design and operational parameters on the plasma dynamics in the source, and finally the plasma heating in the Linac4 H − source in the nominal operational conditions was simulated for the first time, taking into account the neutral populations as well as the external static magnetic fields [15]. The transition between the two stable operating regimes (capacitive and inductive [16]) of the inductive discharges as well as their ignition and stable state are other main issues when RF heating of inductive discharges is studied. Some recent works employing PIC models [12,14,17,18] consider these issues.
The complexity of the physics involved in the extraction process is a result of: sophisticated shape of the 3D electromagnetic field; kinetic reactions and plasma surface processes. The key parameter is the meniscus (Figure 2) -the 3D frontier between the plasma and the beam -as its shape and position determine to a great extent the beam optics. Concerning the beam quality one of the main issues are the co-extracted electrons which induce thermal loads on the electron dumping system, which contributes to increased beam emittance. Formation of beam halo (part of the beam with poor optics) and top-bottom assymetry [19] (caused by the magnetic field filter) are other issues which should be taken into account. Recent works take into account and/or study these effects, as some of the models for simulation of the extraction region are capable to simulate the extraction systems of ITER like NBI in three dimensions and even a first attempt for simulation of the full-size ITER negative ion source is made [9,10,[20][21][22].
The aforementioned models are based on kinetic plasma description, but as it was noted a great deal of the plasma phenomena could be described accurately enough by the fluid theory. Here the concept for matrix source of volume produced H − ions [7], based on fluid models, shall The small discharge radius is in the basis of the matrix source, theoretically motivated by results obtained within two-dimensional fluid models of discharges inductively driven by a cylindrical coil [23,24]. The small radius (2-3 cm) ensures high density of the negative ions in the region with the highest value of the plasma potential [23] as a result of the negative ion flux in the DC electric field in the discharge. Increasing the radius decreases the electric field and respectively the ion flux, and thus the accumulation effect of NI on the discharge axis is diminished [24]. The measured [25] higher density of the NI in the driver region (a smallradius quartz tube) of a two chamber, tandem small-size plasma source, inductively driven by a cylindrical coil confirms the theoretical results. The source design -a matrix of gas-discharge tubes -has orientated the work to planar-coil driving of the matrix [24]. Thus, the suggested construction of the source is for a matrix of small-radius inductively driven discharges with: a single aperture extraction from each discharge; RF inductive driving by a planar coil (a large area planar coil, common for the whole matrix). Finding RF driving which would ensure both efficient RF power transfer to the plasma and equal plasma parameters (including their spatial distribution) in the separate discharge tubes is an issue which comes to the fore in the development of matrix plasma sources. For the purpose a three-dimensional model coupling the gas discharge part within the fluid plasma theory with electrodynamics was developed. Applying the model several tests [26] of different coil designs and source constructions were made. Finally a zigzag coil with an Ω-shaped conductor on the bottom of each discharge tube proved to be an effective coil for RF driving of a source completed as a matrix of small-radius discharges [27].

Microwave excitation of linear array of plasma filled cavities
In this section will be presented the obtained by us preliminary modelling results, for macroscopic plasma parameters and the ϕ component of the electric field, of microwave heated plasma source. The chosen gas as discharge medium is argon due to its simple chemistry in discharges at low pressure. The intended application of the plasma source is in the field of additional acceleration of cyclotron proton beam. The simulated source represents a plasma filled linear array of cavities made of a diellectric material with high permittivity (ε = 25) (Figures 3 and 4). In the current study it is shown that this structure is capable of sustaining a surface microwave discharge with high plasma density, i.e. surface microwaves heating the plasma in the cavities. For the purpose we have implemented the same modelling techniques, as in the study of RF heating of the matrix source [26,27].  Since the geometry of the source is intrinsically axially-symmetric (Figure 3), we reduced the modelling domain to a two-dimensional domain with axial symmetry (Figure 4). On the plasma domain a wall boundary condition (for the fluid equations) is applied to account for the loss of electrons due to their net and random flux to the walls and, respectively, ground (V = 0) for the Poisson's equation. The electric field wave equation has the following boundary conditions: input port -microwave deposited power of 160 W (f = 2.4 GHz TE-wave) fed through the top side of the dielectric; output port on the bottom side; perfect electric conductor -lateral side of the dielectric and on the top and botom sides of the plasma domain. The chemical mechanism in the plasma volume is simplified as it consists only of 3 species and 7 reactions: elastic collisions; direct ionzation via collision with electrons (15.8 eV ); all electronically excited levels lumped into a single excited level (11.5 eV ); excited argon atoms are consumed via superelastic collision with electrons (−11.5 eV ); stepwise (4.24 eV ) and Penning ionization; quenching with neutral argon atoms. Two surface reactions (particle-wall interactions) are taken into account: excited atoms convert to ground state argon atoms and of the argon ions convert to neutral argon atoms. The cross-sections and rate coefficients of the argon volume and surface reactions are taken from the available data within the Comsol Multiphysics plasma interface.
The gas discharge part of the implemented model solves electron continuity and energy balance fluid equations, thus the results obtained for the macroscopic plasma properties are for the distribution of the electron density and temperature Figure 5  electron density obtained is high ∼ 10 19 m −3 , thus the electron plasma frequency is also highf e = (1/2π)(e 2 n e /m e ε 0 ) 1/2 = 63 GHz (where e is the electron charge; n e -the electron density; m e -the electron mass; ε 0 -the vacuum permittivity). The wave frequency is much lower than the electron plasma frequency (f << f e ), thus the wave does not propagate in the plasma, it is absorbed and its energy is dissipated in the region near the dielectric wall. As a result the electron temperature has a maximum in this region ( Figure 5 (b)). The electrons gain energy in the region with maximum electron temperature and carry it in the plasma volume where it is dissipated through collisions, therefore the electron density has its maximum in the plasma volume ( Figure 5 (a)). The well known formula from the literature [28] E max = (m e cω p )/e (cspeed of light (in vacuum); ω p = 2πf e -plasma frequency) gives us the maximum electric field in plasma, i.e. the maximum electric field which could be waked inside the plasma. Since it is a function of the electron density (ω p ∝ n 1/2 e ) the value of electric field in the plasma is high E max ∼ 680 MV/m. The main conclusion of these results is that this structure of the plasma source and the chosen manner of plasma heating ensure high enough electron density which is a requirement for plasma acceleration sources. Next step in the development of the model is to change the discharge gas to hydrogen.