Experimental benchmark of the NINJA code for application to the Linac4 H− ion source plasma

For a dedicated performance optimization of negative hydrogen ion sources applied at particle accelerators, a detailed assessment of the plasma processes is required. Due to the compact design of these sources, diagnostic access is typically limited to optical emission spectroscopy yielding only line-of-sight integrated results. In order to allow for a spatially resolved investigation, the electromagnetic particle-in-cell Monte Carlo collision code NINJA has been developed for the Linac4 ion source at CERN. This code considers the RF field generated by the ICP coil as well as the external static magnetic fields and calculates self-consistently the resulting discharge properties. NINJA is benchmarked at the diagnostically well accessible lab experiment CHARLIE (Concept studies for Helicon Assisted RF Low pressure Ion sourcEs) at varying RF power and gas pressure. A good general agreement is observed between experiment and simulation although the simulated electron density trends for varying pressure and power as well as the absolute electron temperature values deviate slightly from the measured ones. This can be explained by the assumption of strong inductive coupling in NINJA, whereas the CHARLIE discharges show the characteristics of loosely coupled plasmas. For the Linac4 plasma, this assumption is valid. Accordingly, both the absolute values of the accessible plasma parameters and their trends for varying RF power agree well in measurement and simulation. At varying RF power, the H− current extracted from the Linac4 source peaks at 40 kW. For volume operation, this is perfectly reflected by assessing the processes in front of the extraction aperture based on the simulation results where the highest H− density is obtained for the same power level. In surface operation, the production of negative hydrogen ions at the converter surface can only be considered by specialized beam formation codes, which require plasma parameters as input. It has been demonstrated that this input can be provided reliably by the NINJA code.


Introduction
Negative hydrogen ion sources are applied at many particle accelerators over the world for generating highenergy proton beams [1][2][3][4]. Designated experiments and future accelerators impose high demands on the ion source performance which are only partly met by existing H − sources. One example is the ion source of the Linac4 accelerator currently under commissioning at CERN. It is part of an injector chain upgrade of the Large Hadron Collider (LHC) aiming at improving beam brightness and luminosity of the LHC. In routine operation, the ion source of the Linac4 has to deliver an H − current of 45mA within an emittance of 0.25 πmmmrad in pulses of 500μs at a repetition rate of up to 2Hz [5]. However, some particular experimental cases demand a much higher H − current of 80mA, a value that has not been achieved at such low emittance by any existing ion source yet. In order to reach all desired parameters, a dedicated optimization of the ion source is inevitable.
In the ion sources, a low-pressure hydrogen discharge is used for creating negative hydrogen ions. The formation of H − in the plasma can proceed in two different ways [6]: with the first method, the so-called volume process, H − ions are created in the plasma volume via dissociative attachment of electrons to vibrationally excited H 2 molecules. The second method, the so-called surface process, relies on the conversion of hydrogen atoms or ions to H − on a surface with low work function. In the latter case, caesium is evaporated into the ion source for establishing a low work function conversion surface. If the surface process is utilized at ion sources, a higher H − current and a reduced co-extracted electron current is achieved in general compared to pure volume operation. However, due to the high reactivity of caesium the work function of the conversion surface and therefore the source performance degrades with time as caesium is removed from the surface by the plasma and it forms compounds with the background gas. On the other hand, sources operated in volume mode (i. e. without evaporating caesium) have the disadvantage of a lower H − and a higher co-extracted electron current but have the advantage of a higher temporal stability of the source performance and reduced maintenance requirements. The Linac4 H − ion source is typically operated in surface mode reaching a negative hydrogen ion current between 40 and 50mA [7]. In pure volume operation, the achievable H − current is only around 30mA [8] which is below the nominal source parameters.
The creation rate of H − is directly linked to the plasma parameters for both methods: in the volume mode, a high vibrational excitation of ground state H 2 molecules is beneficial as the reaction rate of dissociative attachment increases by orders of magnitude if the molecule is vibrationally excited [9]. The vibrational population of H 2 is in turn influenced by the electron density and temperature. For the surface process, the creation rate of H − is determined by the flux of atomic and ionic hydrogen particles onto the conversion surface [6]. Therefore a higher H − production rate can be expected when the dissociation and ionization degree of the discharge increases (assuming a constant work function of the surface). On the other hand, also the processes leading to a destruction of the negative hydrogen ions like stripping of the excess electron via electron or heavy particle collisions depend on the plasma parameters [10][11][12]. The achievable extracted current of negative hydrogen ions is closely linked to the density of H − ions near the extraction aperture which is determined by the balance of the creation and destruction processes in the discharge. Therefore, improving this current requires an assessment of the actual conditions in the plasma in order to identify performance-limiting processes.
However, diagnostic access to plasma parameters especially in the region close to the extraction aperture is challenging, as the design of ion sources applied at accelerators is very compact. Furthermore, some diagnostic methods like Langmuir probes have problems withstanding the strong heat loads occurring on very short time scales [13]. Optical emission spectroscopy (OES) can ideally be applied at ion sources as only a small window to the plasma chamber is required but the evaluations only yield line-of-sight integrated results. Hence, it is barely possible to assign plasma parameters to the region close to the extraction aperture. They can be accessed more easily via spatially resolved simulations of the discharge dynamics. For the Linac4 H − source, the electromagnetic particle-in-cell Monte Carlo collision code NINJA [14] has been developed for this purpose. The code simulates the inductive RF heating of the discharge and tracks the particle motion in 3D3V (details are given in section 3). Similar codes exist for other ion sources, e. g. for the J-PARC source [15]. In order to check that the simulations reflect the actual discharge conditions and yield correct results, a detailed experimental benchmark is essential. However, such investigations have not been carried out for any of these codes yet.
The investigations presented in this paper address this task for the NINJA code. For a detailed benchmark, experimental access to the spatially resolved plasma parameters is required, which is not feasible at the Linac4 ion source. Therefore, dedicated measurement and simulation campaigns have been performed at the laboratory experiment CHARLIE. Although the RF power levels reachable at this ICP (1 kW max.) are much smaller compared to the Linac4 source, the experiment is ideally suited for the NINJA benchmark: it has a similar geometry than the Linac4 plasma chamber, it is capable of operating in CW within the desired pressure range of a few pascal, and most importantly, it is equipped with multiple diagnostics allowing the determination of the relevant plasma parameters and their spatial distributions. The detailed comparison of measurements at CHARLIE and the corresponding NINJA simulation results is contained in section 4. A further benchmark of NINJA results is performed directly at the Linac4 ion source with line-of-sight integrating OES measurements (see section 5). Based on the spatially resolved NINJA calculations, an assessment of the discharge conditions close to the extraction aperture with respect to the H − production and destruction mechanisms is carried out in section 6.
2. Experimental setup and applied diagnostics 2.1. The Linac4 H − ion source The plasma of the Linac4 ion source is generated via inductive RF coupling inside a cylindrical Al 2 O 3 plasma chamber (see figure 1). The RF amplifier operates at a frequency of 2MHz with a maximum RF power P set of 100kW. A typical RF pulse lasts about 900μs and the repetition rate is up to 2Hz. The RF amplifier is connected to the helical coil (five windings, embedded in epoxy to avoid RF breakdowns between the windings) via a matching network. For the RF power values in this paper, the power delivered to the ion source P del is given, i.e. the power reflected due to a mismatch between the ion source impedance and the 50Ω output of the generator is subtracted from the power set at the generator ( = -P P P del set refl ). The RF current flowing over the antenna, which serves as an input parameter for the NINJA simulations, can be measured with a current transformer. A magnetic cusp field in Halbach offset configuration (formed by alternating magnets with clockwise and counterclockwise magnetization) is generated by NdFeB permanent magnets in order to reduce the loss of electrons to the plasma chamber wall. In addition, permanent magnets installed close to the extraction aperture create a dipole filter field for reducing the electron temperature in that region and thereby the destruction rate of H − via electron stripping. When the source is operated in surface mode, caesium is evaporated in an oven (not shown in figure 1) and delivered to the plasma chamber via a transfer line. The conversion of atomic and ionic hydrogen to H − is intended to happen on the  45 -chamfered Molybdenum plasma electrode. Hydrogen gas is supplied to the ion source with a fast piezo valve opening only for a short time prior to the RF pulse. The gas diffuses through the plasma chamber and is pumped through the extraction aperture. Due to the transient nature of this pulsed gas inlet, it is very difficult to give an exact number for the pressure during plasma phase. In situ pressure measurements are also not possible due to the compact design of the ion source. Dedicated experiments concerning the gas particle propagation inside the vacuum setup of the ion source with respect to the gas valve settings have been carried out by replacing the discharge chamber with a T-shaped flange where a pressure gauge has been installed [16]. From these investigations, the gas pressure has been estimated to 3Pa for the measurement campaigns presented in this paper. As the gas dynamics happens on the scale of several ten milliseconds (whereas the RF pulse only lasts about 900 μs) the gas pressure can be treated as approximately constant during the RF pulse. Figure 2 shows the timing of a typical RF pulse. The timing reference t=0 is the start of the 'useful' part of the beam. Hydrogen is injected 2.5ms before t=0 as it takes some time until the gas density in the plasma chamber is in the desired range of 3 Pa. During an RF pulse there is always beam extracted as the high voltage supply of the extraction system cannot be switched on and off within a few microseconds. This leads to a unstable beam especially during the transient ignition phase of the discharge. Therefore, the RF pulse duration is chosen much longer than the 500μs actually required for the beam but the unstable head of the beam pulse and a short tail are chopped in the low energy beam transport section of the Linac4 [5]. For the pulse shown in figure 2, the RF pulse starts at −400μs and the discharge ignites between −350 and −300 μs. For ignition, a lower RF power is set in order to limit the co-extracted electron current during ignition. Afterwards, the RF power is ramped up to the desired value of = P 40 del kW. The 500μs of 'useful' beam are marked by the grey shaded area.
OES measurements can be carried out via three lens heads collecting light emitted by the discharge (acceptance angle • 3 ). They are separated from the discharge by a sapphire window and focus the plasma emission into optical fibres leading to a spectrometer. The different lines of sight (LOS) determined by the lens heads are tilted with respect to the central axis of the cylindrical discharge vessel: • 0 which means on-axis, • 19 and • 26 . For the measurements presented in this paper only the on-axis view port has been used. The utilized high-resolution spectrometer (1 m focal length, grating 2400 grooves per mm) is equipped with an ICCD camera and has a Lorentzian apparatus profile with a full width at half maximum of l D » 8 FWHM pm at a wavelength of 600nm. The acquisition of spectroscopic data starts at t=0 with typical integration times of 400μs, so the measurement is only taken during the 'useful' part of the beam. The intensity calibration of the system has been performed via an Ulbricht sphere serving as secondary radiation standard.

CHARLIE
The laboratory experiment CHARLIE (Concept studies for Helicon Assisted RF Low pressure Ion sourcEs) is actually intended to investigate the possibility of applying Helicon coupling for RF driven H − ion sources [17,18]. However, during the benchmark measurement campaign for the NINJA code, the experiment was operated as ICP using a helical coil made out of copper tube (6 mm diameter) with 5 windings. A sketch of the experimental setup is shown in figure 3. The discharge vessel consists of a cylindrical quartz vessel having a diameter of 10cm and a length of 40cm. It is attached to the vacuum system at the ends of the cylinder where gas is supplied at a constant flow rate of 5ccm from the left side and the pumps are attached to the right side. The pressure, which is adjusted via limiting the pumping speed, is measured gas-type independently by a temperature-stabilized capacitive pressure gauge before plasma ignition. For the investigations presented in this paper, the experiment has been operated in a pressure range between 0.5 and 5Pa what covers the range relevant for the Linac4 ion source. The coil is connected to the RF generator working at 1MHz with a maximum power of 1kW via a matching network. For all investigations presented in this paper, a perfect match of the load to the 50Ω output impedance of the generator has been achieved by adjusting the variable capacitors in the matching network, i.e. the reflected power equals zero and = P P del set . The impedance matching is monitored in situ with a voltage and current probe placed between the generator and the matching network. Furthermore, the V/I probe is used for a precise determination of P del . The RF rms current I ant flowing over the antenna is measured with an RF current transformer installed around the connection between the matching unit and the coil.  Concerning diagnostic methods for plasma parameters, OES can be performed at two different LOS, both shown in figure 3. One LOS is oriented parallel to and in 1cm distance to the cylinder axis with a diameter of 1cm and virtually no opening angle. It is used for the determination of plasma parameters as described in section 2.3. The second LOS is oriented laterally to the vessel and is movable along and perpendicular to the cylinder axis. This LOS is used for recording lateral intensity profiles of the atomic b H emission line at particular z positions. From these measurements, the radial emission profile of b H is determined via Abel inversion as described in [19]. Prerequisites of the numerical approach are a cylindrical symmetry of the emission and a negligible optical thickness of the emission line. Both requirements are fulfilled due to the geometry of the experiment and due to the low population density of the lower n=2 state of the b H line determined from the CR model Yacora H. The OES measurements at both LOS are performed with a high-resolution spectrometer (Gaussian apparatus profile, l D » 18 FWHM pm at 600nm) which is intensity calibrated with an Ulbricht sphere.
In addition to the line-of-sight integrated results yielded from OES, a movable probe is available. As the plasma has only little contact to metallic vacuum parts (serving as reference point for electric potentials) due to the dielectric discharge vessel, a common Langmuir probe cannot be applied. Instead, a floating double probe is used [20,21]. The probe is inserted through the right end plate at a radial distance of 1cm to the central axis and it can be moved parallel to the cylindrical axis in a range of 20cm. The probe tip consists of two identical tungsten wires (length 10 mm, diameter of 300 μm) which are separated from each other by 6 mm. From the probe characteristics the distribution of the ion density n i and, applying quasineutrality also the electron density n e =n i are derived [21] along the z direction. In general, also the electron temperature can be evaluated but the absolute value is very susceptible for RF distortions typically leading to an overestimation of T e . Therefore, only the relative trends of the obtained electron temperatures values along z are used for benchmarking the simulation.

Determining plasma parameters from OES
The OES measurements carried out both at the Linac4 test stand and the CHARLIE experiment cover the Balmer series of atomic hydrogen (H a to d H ) and the Fulcher-α transition of molecular hydrogen ( P  S + d a u g 3 3 , located between 590 and 650 nm). For evaluating the obtained data, the collisional radiative models Yacora H for atomic and Yacora H 2 for molecular hydrogen [22,23] are applied. These models balance all relevant population and depopulation processes for the particular states of H or H 2 yielding population densities. Processes including the reabsorption of photons have not been considered in order to facilitate the evaluation, which means the plasma is treated as optically thin. As input parameters for the models, the densities and temperatures of neutral (H and H 2 ) and charged particles (electrons and H + , + H 2 , + H 3 , H − ) are required. In order to determine plasma parameters from OES measurements, these input parameters are varied until both the absolute emissivities and the line intensity ratios of the measurements are matched (in the models a Maxwellian EEDF is assumed, the validity of this assumption is discussed later on). For the investigated discharges both at the Linac4 test stand and at CHARLIE, the dominating population processes arise solely from atomic or molecular hydrogen and processes involving the hydrogen ions play only a minor role. This means that the densities of H + , + H 2 , + H 3 and H − cannot be determined reliably from the OES measurements. However, this also has the large advantage that the evaluation is strongly facilitated as the free parameters of the models are reduced to the densities of H and H 2 besides the electron density and temperature.
The rotational and vibrational population of the hydrogen molecule and the gas temperature are determined from the Fulcher-α transition. From its spectrum, the first twelve emission lines (rotational quantum numbers Figure 4 shows an exemplary emission spectrum measured at the Linac4 test stand where the relevant emission lines are indicated by drop lines. For hydrogen, the single rovibronic emission lines can be well separated, hence the population of the rotational states can directly be obtained from the emissivities.
In low pressure hydrogen discharges a non-Boltzmann rotational distribution is typically present [24], which can be approximated by a two-temperature distribution [25]. The low-lying rotational levels are described by the cold part of this distribution according to the temperature T rot,1 reflecting the population via heavy particle collisions. The hot part of the rotational population describes the high-N levels according to T rot,2 . It is weighted relatively to the cold part by the factor β for fitting. There are several possible reasons behind the hot population: recombinative desorption of hydrogen atoms at the wall of the discharge vessel, dissociative recombination of + H 3 with electrons, or direct electron impact excitation [24]. For the high-power Linac4 discharge, the most likely process leading to the hot rotational distribution is collisional excitation by electrons. For the CHARLIE experiment, it is recombinative desorption of hydrogen atoms at the wall where a part of the binding energy is converted in rotational excitation. The gas temperature of the discharge can be determined by projecting T rot,1 from the excited P d u 3 state into the electronic ground state of the hydrogen molecule according to the rotational constants of the two states [26].
The vibrational population is determined from OES by summarizing over all measured rotational states within one vibrational level. The obtained vibrational distribution in the P d u 3 state is then fitted with a population calculated via a corona model. In this model, a Boltzmann vibrational distribution is assumed in the electronic ground state of H 2 which is transferred to the excited state by electron impact excitation. The dominant loss channel is spontaneous radiative decay. It has been proven that a corona model describes the equilibrium population of the P d u 3 correctly, as contributions from other population channels like cascades from higher lying levels can be neglected [27]. Finally, knowing both the rotational and vibrational distributions of the P d u 3 state, the integrated emissivity of the Fulcher-α transition can be calculated. It is required for the previously described evaluation of the CR model Yacora H 2 .

The NINJA code
The goal of the simulations is to investigate the plasma dynamics in the inductively coupled discharge of the Linac4 as a result of the RF heating process. This requires a self-consistent solution of the coupling between the electromagnetic field generated by the RF coil and the plasma response. NINJA is a fully implicit electromagnetic particle-in-cell Monte Carlo collision (EM-PIC-MCC) code designed for the simulation of RF-ICPs in cylindrical configuration. A detailed description can be found in [14], therefore only a brief summary is given here. The code is composed of three modules: an EM-PIC module for the solution of the EM field and the motion of charged particles, a MCC module taking into account the collision processes between the different species, and a neutral transport module aimed at tracking atomic and molecular species (the ground state of H 2 vibrationally resolved). The EM-PIC module is 2.5 D in cylindrical coordinates, i. e. the EM field calculation is performed in 2D assuming azimuthal symmetry, while the particle dynamics is fully 3D3V. The simulation domain used for modelling the Linac4 discharge is shown in figure 5. The plasma chamber itself is a cylinder of radius R=24mm and length L=138mm (see section 2.1). The calculation of the EM field is performed on a larger domain in order to take into account the coil and to avoid large reflections at the boundaries. The RF coil is simulated as a perfect conductor, with a rectangular cross-sectional area (as indicated in figure 5) covering the actual coil dimension. For the consideration of the external static magnetic cusp and filter fields, a B-field simulation with the 3D Tosca module of the OperaVectorfield software application is carried out beforehand and imported to the PIC-MCC via a field map. The external magnetic field is interpolated at the particle level to maintain its 3D features. For the simulations of the CHARLIE plasma, a similar domain with dimensions adjusted according to the actual geometry (see section 2.2) and without plasma electrode is employed.  The simulated plasma includes five charged species: electrons, H + , + H 2 , + H 3 and H − , all tracked kinetically in 3D3V. At the plasma chamber wall, absorbing boundary conditions are employed, removing from the computation all charged particles hitting it. The governing equations (Maxwell and Newton equations of motion) are solved with an implicit time integration scheme. This allows the cell size to exceed the Debye length and the time step to be larger than the plasma period while preserving the conservation of the total energy.
The MCC module includes a database of over 200 cross sections taking into account the most important processes in hydrogen plasmas (e. g. elastic collisions, ionization, electronic and vibrational excitation via collisions, see [14] for details and references). Electron-neutral, electron-ion and ion-neutral collisions are handled by a null-collision method [28] while Coulomb collisions are treated according to the binary method [29]. The neutral transport module kinetically tracks atomic and molecular populations (the latter vibrationally resolved). To cope with the large number of neutrals compared with the charged particles (the ionization degree is typically only few per cent) NINJA uses a variable weight scheme with a rezoning technique. In a given cell, neutral particles are merged and split in order to keep a constant number of particles per species per cell (typically 100). At the plasma chamber wall, pure reflection for molecules and a user-defined coefficient for recombinative desorption of hydrogen atoms, which is set to 0.5 for the Molybdenum plasma electrode and tó -2.5 10 3 for the ceramic walls [30], is employed.
As input, the code requires initial densities and temperatures for the considered particles. Furthermore, the RF current flowing over the antenna needs to be set. The value of the current, which is taken from measurements, is kept constant during the simulation. It should be noted that the specification of the RF current for the simulation includes an implicit assumption as actually, the RF power is specified in the experiment and the current adjusts to the RF circuit characteristics. However, for high-power discharges like the Linac4 plasma, a strong inductive coupling is present where the antenna current is directly proportional to the applied RF power: increasing the electron density by increasing the power level leads a higher antenna current as the equivalent plasma resistance stays constant [31,32]. In the NINJA code, the intermediate steps of calculating the equivalent plasma resistance from the plasma parameters and determining the resulting current afterwards are omitted for the sake of simplification. Including these steps would require the consideration of the full RF circuit in the simulation and a feedback loop between the plasma parameters and the characteristics of the RF circuit for each time step.
The simulation outputs cover the density and energy distributions for all simulated species, the plasma current, the EM field, and the vibrational distribution in the electronic ground state of H 2 . In addition, the code output contains information on its numerical performance and on the satisfaction of conservation laws.

Benchmark of the NINJA code at the CHARLIE experiment
The benchmark of the NINJA results by measurements at CHARLIE is carried out for a variation of the gas pressure between 0.5 and 5Pa at a fixed power of 520W. Additionally, the RF power is varied between 250 and 800W at a pressure of 3Pa. For all these experimental conditions, the cell size of the NINJA code is set tó 2.5 5 mm (D´D r z) and the time step to D =´t 5 10 11 s in order to resolve the gradients in the plasma skin depth, which is about 1.7cm, and to provide sufficient statistical sampling of the collision frequency, which is about 10 8 Hz, respectively. Independent on pressure or RF power, the plasma is initially loaded with a density of 10 17 m −3 uniformly in the whole plasma chamber and with an equal distribution between the positive ions H + , + H 2 , and + H 3 , using 8 million particles in total. The initial temperatures are 1eV for the electrons and 0.1eV for the ions. The value of the current flowing over the RF antenna is set to the one measured experimentally.
The plasma dynamics are simulated for 30μs and within the last 5μs, the electron and ion densities do not change any more, i. e. the discharge reached steady state as shown exemplarily in figure 6 for 0.5Pa pressure and 520W power (corresponding to an RF current of 34.0 A). Two weeks calculation time are required on a 16-core cluster for performing such a simulation. As already described, CHARLIE is operated at much lower power levels than the Linac4 plasma. Accordingly, the achieved electron densities are distinctly lower. Concerning the simulation, this has the direct implication that dissociation and recombination processes which determine the degree of dissociation and processes leading to a redistribution of vibrational states in the H 2 molecule are not yet in equilibrium after 30 μs. For the low electron density conditions present in the CHARLIE plasma, these processes would require a much longer simulation time in the range of milliseconds. As such long time scales cannot be simulated with the NINJA code with reasonable time effort, the dissociation degree and the vibrational temperature are not calculated self-consistently and the initial parameters are chosen according to the experimental results. In figure 6 it can be seen that the densities of H and H 2 stay virtually constant at the value initially chosen from the experimentally determined dissociation degree of about 17% during the whole simulation whereas the other particle densities evolve. In general, the simulation results do not depend on the initially chosen specific values for the densities and ion distributions with the exception of the H and H 2 densities and the vibrational distribution in the H 2 ground state.
For determining the steady-state values of the plasma parameters, an averaging over the last two RF cycles (which is equal to 2 μs) is performed for the time resolved values. Figure 7 shows the corresponding spatially resolved steady-state results for the electron density and mean electron energy obtained from the NINJA simulations for a pressure of 3Pa and an RF power of 520W corresponding to an antenna current of 28.5A. The electron density shows a pronounced maximum at the discharge centre below the RF coil where values of up to7 10 17 m −3 are reached and decreases towards the wall both in axial and radial direction. For the mean electron energy, the profile in z direction is also characterized by a maximum below the antenna. In radial direction however, a minimum on the cylinder axis is observed under the antenna and two peaks occur close to the wall of the discharge vessel (in positive and negative r direction). These peaks originate from the RF heating taking place in that area. Towards the end plates of the vessel, i. e. approaching z=0mm or z=400mm, the two peaks vanish as the heating only occurs under the antenna.
In order to allow for a direct comparison of the simulation results with those results obtained from OES, the line-of-sight averaging implicitly performed by OES has to be taken into account. For that purpose, those electron density or mean electron energy values are averaged, which are located inside the capture volume of the light collecting optics. The corresponding area is indicated by the top inlets in figure 7 showing the projection of the spatial distributions. Furthermore, a crucial point for evaluating OES measurements by CR models is the EEDF of the discharge. Within the whole plasma volume, the simulated EEDFs show virtually no deviation from a Maxwellian one, which justifies the application of a Maxwellian EEDF in the CR models.
For comparing the z profiles of the simulated and measured electron densities and energies, the floating double probe is applied. The probe measures at a distance of 10mm from the central axis and the distance between the two probe tips is around 6mm. For one specific z value, the simulation data is therefore averaged in radial direction between 5mm   r 15 mm. The probe measurements have been performed between 180mm   z 320 mm.  Measuring at lower z values than 180mm resulted in a strong distortion of the discharge by the long probe rod. Therefore, the values from z=80mm to z=170mm are obtained by mirroring the corresponding high-z range (it has been checked previously that the profile is symmetric). The left part of figure 8 shows the resulting z profiles of the electron density both for the NINJA simulation and the probe measurements where a pronounced central maximum and a decrease towards the end plates is present. The agreement between simulation and probe data is very good except for the central peak, which is higher in the simulation. This can be explained by the distortion of the RF heating when the probe is moved under the antenna what lowers n e .
The analogous comparison of the electron temperature is presented in the right graph of figure 8. As described in section 2.2, only the relative trends of the probe are reliable, therefore the maximum of the profile determined by the probe is normalized to the maximum obtained from the simulation. Also for T e , an excellent agreement of measurement and simulation (the electron temperature is derived from the mean electron energy after = á ñ T E 2 3 e e ) is observed. Concerning the r profiles of the electron density and temperature, a direct comparison can unfortunately not be carried out as the probe cannot be moved in r direction. Therefore, the validity of the simulated r profiles is examined indirectly by a comparison with measured radial emission profiles of the atomic b H emission line ( =  = n n 4 2 ) obtained from Abel inversion (see section 2.2). The simulated emission profile is calculated from the simulated profiles of ( ) is derived from cross section data [12]. In the left part of figure 9, the radial profiles of the electron density and temperature provided by NINJA are shown at the axial position of =  z 250 5 mm (directly right of the RF antenna) where the Abel inversion measurements are carried out. The electron density profile can be well described by the plotted Bessel profile, which is typically found in cylindrically symmetric discharges [33].
The right part of figure 9 shows the corresponding simulated emission profile together with the one evaluated via OES and Abel inversion. Both profiles show a steep increase of the emission at the cylinder walls  followed by a broad maximum in the vessel centre. There are small deviations present, as the profile determined from the NINJA simulation shows a small peak at r=25mm whereas the profile determined from Abel inversion shows a small peak in the discharge centre. However, it should be noted that Abel inversion is very sensitive on slight changes of the measured line-integrated emissivities resulting in the appearance or absence of small peaks in the radial emissivity profile. Therefore, rather the general shape of the profiles and not the small details should be compared. Bearing this in mind, a good general agreement is observed.
For all investigated variations of pressure and RF power, the basic spatial distributions of electron density, energy and emissivity of the b H line described above are similarly observed. Therefore, only the derived line-ofsight averaged values are considered in the following. Figure 10 shows the comparison of n e and T e at varying pressure between 0.5 and 5Pa for an RF power of 520W. For the electron density, the measured data obtained by probe and OES are in excellent agreement and depict a steady increase with increasing gas pressure. In general, the data provided by the simulation matches the absolute values closely but the trend with increasing pressure is opposite than the measurements. Concerning the electron temperature, both the measured and simulated results show a steady decrease with increasing pressure, in agreement to the expected behaviour in low-pressure discharges. However, the values obtained by OES are systemically lower by about a factor of two than the simulation. Figure 11 shows the benchmark results obtained for 250, 520 and 800W at a constant gas pressure of 3Pa. The electron densities measured via OES and probe agree well and exhibit the expected increase with increasing power. In contrast, the simulation shows no systematic trend with increasing pressure as the density remains virtually constant. There is no pronounced influence of the power on the electron temperature which remains virtually constant for both OES measurement and simulation (with the exception of the increased temperature at 800 W in the simulation). However, similar to the results of the pressure variation, the values determined from OES are systematically about a factor of two lower than those provided by the simulation.
The observed deviations between measurement and simulation is likely to arise from the assumption made in the simulation concerning the correlation between RF power and current flowing over the antenna. As described in section 3, NINJA assumes strong inductive coupling, which means that increasing the RF power results in a higher current flowing over the ICP antenna. This assumption is valid for high-power discharges like the Linac4 plasma. However, the discharges investigated at CHARLIE are operated at moderate RF power. This regime is better described by the characteristics of a loosely coupled ICP where the antenna current is not necessarily proportional to the input power [31,34]. Increasing the electron density by raising the RF power level results in a higher equivalent plasma resistance [32] leading to a reduced antenna current and hence to the opposite behaviour than assumed in the NINJA code. For CHARLIE, the measured currents show a decrease from 31.8 to 28.5A when the RF power is increased from 250 to 520W at 3Pa. When the power is raised further, the currents (which directly serves as input parameter for NINJA) stays virtually constant (see [35,36] for more details).
The RF current strongly influences the discharge properties as it specifies the magnitude of the RF field and has therefore a direct impact on both the electron temperature and density. The deviations between simulation Figure 10. Comparison of the simulated and measured electron density and temperature for a variation of the pressure at a fixed RF power of 520W at CHARLIE. For all parameters, axially averaged values are presented. and measurement at CHARLIE can be explained by the missing feedback between the plasma parameters and the RF circuit in the NINJA code. It is foreseen to include this feedback-loop in a next improvement step of NINJA as it requires the consideration of the full RF circuit and hence a major modification the code. Nevertheless, the conducted simulations predict the correct spatial distributions of the electron density and temperature for the CHARLIE discharges. It has not been shown in detail, but this is the case for the other operational parameters, too. Also the absolute values of n e and T e are matched within a factor of two what can still be considered a reasonable agreement.

Benchmark of the NINJA code at the Linac4 ion source
The comparison of NINJA simulations with the results from OES evaluations at the Linac4 ion source is carried out for RF powers of 30, 40 and 50kW what corresponds to 180, 200 and 220A of RF current flowing over the antenna. These power values have been chosen as the optimum source performance is typically obtained at 40kW both for surface and volume operation. For simulating these discharges, the cell size of the NINJA code is set to1 2 mm (D´D r z) and the time step to D =´t 2.5 10 11 s. The initial parameters of the simulation are chosen as follows: a plasma density of 10 18 m −3 consisting of 90% n e and 10% of H − , positive ion distribution 94% H + , 3% + H 2 , 3% + H 3 , electron temperature 1eV, ion temperature 0.1eV, gas temperature 300K. Figure 12 shows the temporal behaviour of the particle densities. The volume-averaged densities have reached a stable value after 20μs. In contrast to the simulation of the CHARLIE plasma, dissociation and recombination processes as well as the vibrational distribution of the H 2 ground state can be calculated selfconsistently as the reaction rates are much higher for the parameters present in the Linac4 discharge. This can be seen in figure 12 where the H density starts from zero and evolves to a steady-state value. The calculated results are independent of the exact value of the input parameters, but they have been chosen from the experience of previous results to have a quicker convergence of the simulation. For determining steady-state plasma parameters, the time resolved values are averaged over the last two RF cycles (which is equal to 1 μs). Figures 13 and 14 show the spatial distributions of the electron density and mean electron energy in the Linac4 plasma chamber for a coil current of 180A. It can be seen that both the electron density and energy peak in z direction at the position where the antenna is located (between z = 56 and 84 mm). Their values drop towards the end plate (visible in the right part of the figures 13 and 14) and the extraction aperture (visible in the left part of the figures) whereby the decrease towards the plasma electrode is much more pronounced due to the influence of the magnetic filter field. In radial direction, the maximum of n e is located on the axis of the cylindrical discharge vessel. The decrease towards the wall is much faster than the broad Bessel profile observed in CHARLIE due to the confining impact of the cusp field. Looking radially from the central axis towards the wall, the electron energy shows a peak similar to the CHARLIE plasma. However, the high-energy values are located much closer to the central axis what can also be attributed to the cusp field. Concerning the EEDF, virtually no deviation from a Maxwellian distribution is observed in the whole plasma volume, which justifies the assumption of a Maxwellian EEDF for OES evaluations of the Linac4 plasma.    values are considered in the following. A comparison of n e and T e obtained from OES and the NINJA code is shown in figure 15. For both simulation and measurement, a good agreement both in the absolute values as well as in the observed trends is evident.
Concerning the rotational population obtained from OES, also the Linac4 discharge exhibits the twotemperature distribution. Furthermore for this plasma, the share of the hot rotational population can be quite high [37] making it difficult to determine T rot,1 and therefore T gas accurately as only very few levels are dominantly populated according to the cold population. A dedicated OES measurement campaign using deuterium instead of hydrogen revealed that the gas temperature is virtually equivalent to the ambient temperature of 300K for all chosen parameter settings [37]. Due to the higher mass of D 2 the energy difference between the particular rotational levels is smaller what makes the cold distribution much more evident compared to H 2 . Also for the simulation, the gas temperature remains at the input value of 300K. This can be explained by the short plasma pulse duration of 900μs preventing the heavy particles from heating up considerably.
For the vibrational distribution in the H 2 ground state, the simulation demonstrates that a two-temperature distribution evolves [38], which is also observed in other low pressure discharges [39,40]. The cold temperature is in the range between 3500 and 4500K, which is in agreement with the measurement. The hot temperature describing the high-v states reaches values of more than 10 000K. This high population of the levels with high vibrational quantum number is especially beneficial for the volume production of negative hydrogen ions. However, these states cannot be probed with OES as only the levels ¢ =  = v v 0, ... 3 have a sufficient signal-tonoise ratio.
Finally, figure 16 shows the densities of atomic and molecular hydrogen determined via OES or simulation. Again, the results from OES and NINJA agree well although the simulation underestimates the H density at 30kW and the observed trends are not reproduced: the simulation shows an increase of the H density froḿ 1.2 10 20 to2.5 10 20 m −3 with increasing RF power whereas the values obtained from OES stays virtually constant around2.4 10 20 m −3 . Concerning the H 2 density, the picture is equivalent with a decreasing trend for the simulation. These slight deviations may arise from the assumption of a constant rate of recombinative desorption of H at the wall in the simulation. The coefficient γ describing this rate determines the H density to a large extent as volume recombination plays only a minor role due to the small vessel dimensions. The value of γ is influenced by many parameters, for example by the coverage of the wall by adsorbed H atoms. Therefore, the exact value may vary for different operational parameters, but the surface parameters cannot be accessed in full detail making a self-consistent consideration of γ in the simulation impossible.
In general, the degree of dissociation is in the range of 15%-30% what is rather low for such high power levels. This can be explained by the large coefficient of recombinative desorption of hydrogen atoms at the Molybdenum plasma electrode. As it is about three orders of magnitude higher compared to the coefficient of the ceramic plasma chamber wall (see section 3), it dominates the surface assisted recombination rate of hydrogen atoms to H 2 molecules despite the smaller area leading to a low degree of dissociation. Figure 15. Electron density and electron energy evaluated from OES and simulated by the NINJA code for a variation of the RF power at the Linac4 ion source.

Assessment of simulation results for the Linac4 ion source
The performance of the ion source is largely determined by the rates of the particular H − production and destruction processes in the plasma region close to the extraction aperture as the negative hydrogen ions can only move about 1 cm in average before they are destroyed. The spatially resolved particle density and energy results obtained from the NINJA simulation can now be used for assessing the processes which limit the achievable H − density in that region with respect to the desired operational mode of the source, volume or surface. For evaluating the processes, the NINJA simulation results are averaged within the chamfered area of the plasma electrode between z=122mm and z=126mm in order to smooth the data. All obtained plasma parameters like the electron density or the dissociation degree show an increase with increasing RF power, only the average electron temperature stays virtually constant at 1.5eV. This value is much lower than the bulk T e (see figure 15) what underlines that the magnetic filter field actually cools the electrons to the desired value.
Concerning the volume operational mode of the source, the averaged plasma parameters are used as input for a zero-dimensional rate balance model, which allows for a quick and easy assessment of the processes. The production of negative ions is considered via dissociative attachment and their destruction by collisions with electrons, positive ions, hydrogen atoms or hydrogen molecules (a detailed description of the model can be found elsewhere [11]). The assessment of these processes for the investigated three RF powers predicts the highest equilibrium H − density and therefore the optimum source performance for 40kW. This reflects the experimental observation perfectly. The model demonstrates that the increasing performance at low power levels results from an enhanced H − production rate due to the increased vibrational temperature. Above 40kW however, the balance shows that the increased dissociation degree reduces the H − density again twofold: first, it reduces the number of molecules what lowers the production rate despite a higher vibrational temperature. Second, the destruction rates by associative and non-associative detachment are increased due to the increased density of atomic hydrogen. Concerning the destruction processes of negative ions in general, the dominant processes are collisions with atomic hydrogen. They contribute to the total destruction rate with a share of about 65%. In addition to that, mutual neutralisation of H − ions due to collisions with positive ions has a share of roughly 25%. Collisional detachment of H − with H 2 is negligible with a share of less than 1%. Finally, the detachment of negative ions by plasma electrons contributes by 10% to the total rate. This low contribution is a direct impact of the filter field.
When the source is operated in surface mode, caesium is evaporated into the plasma chamber. Comparing the OES measurement results with and without caesium present in the source showed no influence of caesium neutrals or ions on the discharge characteristics. Therefore, the results of NINJA calculations, where caesium particles are not considered, can also be used for caesiated plasmas. However, the assessment of the surface operational mode requires including the conversion of atomic or ionic species into negative hydrogen ions on the Molybdenum plasma electrode. For a correct consideration of these processes, a very high spatial resolution especially in the sheath region between the plasma and the plasma electrode is necessary. Such a high resolution cannot be provided by the NINJA code, as it has to cover a very large simulation domain. For addressing this task, specialized beam formation codes simulating only a small area near the extraction aperture are required (e.g. [41,42]). In general, for surface operation, an increased RF power leads to an enhanced destruction rate of H − by collisions with H atoms or positive ions due to the increased degree of dissociation and ionization. However, the production rate of negative hydrogen ions is increasing simultaneously. As the production rate is furthermore strongly dependent on the work function of the conversion surface, which is determined by the caesium coverage, detailed balancing with a simple model cannot be carried out.

Conclusion
The electromagnetic particle-in-cell Monte Carlo collision code NINJA has been developed for the Linac4 H − ion source at CERN for gaining insight in the spatial plasma characteristics. A benchmark of the code has been carried out at the experiment CHARLIE as it is equipped with multiple diagnostics allowing for a detailed comparison of the NINJA simulation results with the measurement.
In general, a very good agreement between simulation and measurement is obtained at CHARLIE for the relative spatial distributions of the electron temperature and density. Concerning the absolute values, the maximum deviation of n e is about a factor of two but the relative trends for pressure and RF power variations are not reproduced by the code. In contrast, the trends are well reproduced for T e but the absolute values are systematically about a factor of two higher than the measured ones. The observed deviations can be explained by the assumption of strong inductive coupling in the NINJA code. For CHARLIE, this assumption is not fully valid as it operates in the loosely inductive coupling regime due to the moderate RF power levels. A correct consideration of this regime would require an implementation of the back coupling between plasma and RF circuit for each simulation time step, which is foreseen for a next improvement step of the NINJA code. Nevertheless, the general agreement between the code and the measurement can be considered good. Furthermore, the Linac4 ion source operates in the strong inductive coupling regime as much higher RF power levels compared to CHARLIE are applied.
Accordingly, a very good agreement between the NINJA simulation results and OES measurements with respect both to the absolute values and the observed trends of T e , n e , T gas , and the densities of atomic and molecular hydrogen for varying RF power is obtained for the Linac4 plasma. Based on the spatially resolved results obtained from the NINJA code, an assessment of the production and destruction processes of H − has been carried out for the plasma region near the extraction aperture. When the source is operated in volume mode, the assessment predicts the highest H − density in the region of the plasma electrode for an RF power of 40kW. This is in excellent agreement with the experimental observation where the best source performance is achieved at this power level.
When the source is operated in surface mode, the H − production is based on the conversion of atomic and ionic species at the caesiated surface of the Molybdenum plasma electrode. This conversion process is not considered in the NINJA code, as it cannot resolve the small-scale physics happening on the surface. Such a task can only be fulfilled by specialized beam formation codes, which in turn require the densities and energies of the various species in this plasma region as input, which are not known a priori. As it has been demonstrated in this paper, NINJA can provide exactly this data reliably also for caesiated discharges as the caesium particles do not alter the plasma parameters. Therefore, the NINJA code establishes the link between the experimentally chosen operational parameters of the ion source and the plasma properties required for the beam formation calculations.
In summary, it has been shown that the NINJA code allows an assessment of the plasma physics taking place in the Linac4 ion source. The code can now be used for systematic and predictive investigations concerning the improvement of the source performance by studying the impact of changing operational parameters or of the ion source design. It should also be noted, that the code is not limited to simulation of the Linac4 ion source but it can also be applied to other high-power inductively coupled ion sources used in the accelerator community.