The particle tracking code BBCNI for large negative ion beams and their diagnostics

Heating and current drive in the next generation tokamak ITER requires the use of large and powerful neutral beams, generated by a precursor ion beam from an ion source around 1 m × 2 m in cross-section. To avoid energy losses and component damage, strict requirements are placed on the divergence and uniformity of this ion beam, which is comprised of many individual beamlets. Understanding, controlling, and predicting the behaviour of these large ion beams requires knowledge of these individual beamlets and their interactions with one another. This is hindered by available experimental diagnostics on these large beams typically only having access to volume averaged information. A forward simulation of beam diagnostics would allow the connection of experimental results with otherwise unobtainable individual beamlet properties. The particle tracking and ray tracing code Bavarian Beam Code for Negative Ions was developed for this reason, and takes into account the interaction of individual component beamlets with whole-beam diagnostics to produce synthetic data that can be compared with experimental results. In this work a significantly reworked and upgraded version of the code is presented and example results are given and analysed for the ITER relevant test facility BATMAN Upgrade. It is shown how the simulation can recreate experimental results, and that one must consider the whole beam in order to do so. The impact of beamlet mixing on beam emission spectroscopy results is shown, as is the importance of long range magnetic fields on the beam transport. The capabilities and limitations of the code are discussed with a view toward application to ITER size ion sources.


Introduction
The operation of large tokamaks, such as the ITER facility currently under construction, calls for correspondingly large heating and current drive capabilities. Neutral beam injection (NBI) plays a key role for both heating and current drive in all large tokamaks, as the injection of high energy H or D atoms allows the deposition of both energy and momentum into the tokamak plasma. 33 MW of NBI power is foreseen to be installed across two beam lines at ITER [1]. To achieve this high power, the beams are required to have a high particle energy of up to 1 MeV for D and 870 keV for H. Due to the achievable neutralisation efficiencies, a beam with such a high energy will need to be created from a precursor beam of negative ions [2]. In order to ensure a high efficiency of energy transfer into the tokamak plasma, and to avoid damage to components along the beam line, strict requirements are placed on the beam divergence (<7 mrad) and beam homogeneity (>90%) [2]. The power and energy requirements correspond to extracted current densities of 28.5 mA cm −2 of 14 mm diameter in the first grid, corresponding to a total extraction area of roughly 0.2 m 2 .
These extreme parameters are required to be sustained for 1 h in D and 1000 s in H, at a low source pressure of 0.3 Pa. The low pressure is necessary to reduce the amount of socalled 'stripping', the neutralisation of partially accelerated ions by background H 2 particles within the acceleration system, which can significantly reduce the efficiency of an NBI system, and may cause excessive heat loads and damage to beamline components from high energy electron impact. The generation of a negative ion beam with the required properties is not a simple task, and no facility yet exists that can concurrently meet all of the requirements. Although recent progress has been made towards the necessary high ion energies [4] and current densities [5], the control of beam divergence and homogeneity for large sources at high ion source powers and low pressures is still a significant challenge [3].
The beam (or beamlet) divergence is a measure of the spread in particle velocity perpendicular to the beam (or beamlet) direction. This is usually assumed Gaussian in distribution, and a divergence in degrees can be calculated from the e-folding half width ( ) s º 2 of this Gaussian [6]. This can be applied to both individual beamlets, or to the beam as a whole, for which different effects need to be considered. Within a single beamlet, the divergence is reasonably well defined, and the main difficulty with this approach comes from any deviations that the velocity distribution may have from a Gaussian. However, when measuring the divergence of a whole beam, one needs to consider the impact of measuring multiple beamlets, each of which may have unique particle spatial and velocity distributions due to inhomogeneities, magnetic fields, or grid tilting [7] and beamlet steering [8], as described in detail in section 2. The divergence measured for the whole beam will therefore always be larger than the minimum value for a single beamlet, but the precise relationship between the two will depend on the properties of the ion source and the operational parameters, and no diagnostic is available that gives this information.
The whole-beam divergence depends on a number of key source parameters, including the ratio of the voltages between each of the grids in the acceleration system. This ratio determines the shape of the electrostatic lenses formed between the grids, and thus the trajectories ions take as they are accelerated. The normalised perveance also has a large impact on the divergence, which is defined as the perveance P = I U ex ext 3 2 normalised to the space-charge limited maximum possible perveance for planar grid surfaces [9], where I ex is the total extracted ion current and U ext is the extraction voltage across the first gap. The perveance is important for the divergence because the space charge of the ion beam affects the shape of the electrostatic lenses, in turn affecting the formation of the beam. In both the normalised perveance and ratio of extraction to acceleration voltage, there is an optimum value (depending on the system) for which divergence is at a minimum. This latter relationship of divergence with the perveance is why the homogeneity of the beam also influences the measured divergence. Unlike sources for positive ions, negative ion sources require the use of transversal magnetic fields to operate, as described in more detail in section 2. These fields cause the plasma within the source to drift in a perpendicular direction, resulting in inhomogeneous plasma densities in front of the extraction system [10]. This has the effect of increasing the extracted current for those apertures where the plasma density is higher, causing corresponding inhomogeneities in the overall beam current density. Due to the relationship between divergence and perveance, the divergence of these beamlets is also affected [11].
To reach the high extracted current densities required, caesium, as a low work function material, is used to increase the surface production of negative ions and thereby their density near the extraction system [12]. It is evaporated into the source, where it is redistributed onto surfaces by the plasma, thereby reducing the surface work function [13,14]. The plasma inhomogeneities mentioned above therefore also have a second order effect of causing non-uniform distribution of caesium, which affects the surface production of -H or -D from the atomic species, and subsequently the extracted current and divergence of each aperture. Despite these issues, low values of divergence have been measured at a number of facilities. A 900 keV negative ion beam was extracted from 15 apertures with an estimated beamlet divergence of less than 7 mrad (0.4°) at the MAMuG facility [15]. The NBI system on the Large Helical Device has demonstrated large beams with an estimated 4-6 mrad (0.23°-0.34°) beamlet divergence at 120 keV, with 5 grid segments each containing 168 apertures [16]. These measurements were both performed using structured carbon based calorimeters. Due to the high power densities involved (often above 1 MW m −2 ), this type of diagnostic cannot be used on the beam at full power, energy, and duration. Instead, diagnostics are limited to non-invasive diagnostics, or diagnostic calorimeters large distances (on the order of metres) away from the acceleration system [17]. Other smaller or more accurate diagnostics would be destroyed by the full beam power and duration. An overview of possible diagnostic systems is given in [17].
This limitation poses issues for the understanding of diagnostic results of large and high power ion beams, as only spatially averaged information is available [18]. Calorimeters can measure absolute beam powers, with limited spatial resolution, and only after beamlets are well merged [19]. Beam Emission Spectroscopy (BES) gives access to the the divergence of beam particles, but only within a defined line of sight. For spatial information, multiple lines of sight can be used to measure vertical or horizontal profiles of whole-beam divergence and relative beam intensity.
In order to accurately interpret the experimental data available for large, high power negative ion beams, and to link it with the formation of beamlets, a simulation is required that can forward calculate whole-beam diagnostic results from the generation of individual beamlets. The simulation tool EAMCC for synthetic calorimetry has been applied to smaller sources [20,21], or small sections of large sources [22], and includes the beamlet formation process. However other simulation work on negative ion beams is limited to single (or small numbers of) apertures, and only for short distances after the grid system, due to the computational intensity of solving Poisson's equation. Presented in this work is the Bavarian Beam Code for Negative Ions (BBCNI), a particle tracking and ray tracing code originally developed at IPP a number of years ago [23,24], but recently heavily upgraded. The code is able to simulate large ion beams and synthetic beam diagnostics, including electrical measurements, calorimetry, and BES. It is the first such tool that can provide beam emission spectra from particle level information, to give a full suite of synthetic diagnostics on a negative ion beam. The simultaneous modelling of the beam and its diagnostics can be used to bridge the gap that exists between the experimental parameters of ion sources and the information that can be obtained from their beam diagnostics. The improved understanding of the beam divergence, uniformity, and transport will allow for the suggestion of specific improvements to the ion sources, and help to bring large negative ion sources towards their targets for ITER and future projects.

Generating and diagnosing large negative ion beams
The generation of a large -H or -D beam requires the creation of a H 2 or D 2 plasma. There are a number of suitable methods for plasma generation. The issues faced are similar between systems, but this work will focus on those using RF driven plasmas, as is foreseen to be used for ITER [1]. A top-down schematic of an idealised RF driven negative ion source as used in an NBI source for fusion applications is shown in figure 1. In these sources, plasma is generated using one or more RF 'drivers'; cylindrical chambers with an RF coil around the outside for inductive RF coupling into the plasma. Following generation of the plasma, it flows into an expansion region, then on to a series of ion extraction and acceleration grids with multiple individual apertures.
A magnetic filter field is used in the expansion region to reduce the electron temperature from around 10 eV near the drivers, to roughly 1 eV near the extraction system, which reduces the electron impact destruction of negative ions. The electron density is also decreased from 10 18 to -10 m 17 3 , which reduces the co-extracted electron current [25]. This filter field is generated with either permanent magnets mounted on the sides of the expansion region, or a DC current driven vertically through the first grid ('Plasma Grid', PG), or a combination of the two [26]. In this case, this field has the main component horizontally across the expansion region, perpendicular to the direction of plasma expansion. The field strength near the PG is on the order of a few mT, which is enough to magnetise the electrons fully and reduce their transport perpendicular to the magnetic field. A side effect of this field is that the plasma within the expansion region is no longer homogeneous, but typically drifts in the vertical direction [27,28].
As electrons are extracted alongside the negative ions, permanent magnets embedded in the second grid, the Extraction Grid (EG), are arranged to create a magnetic 'deflection' field whose main component is in the vertical direction, and alternates row-by-row. This deflection field causes both ions and co-extracted electrons to be deflected sideways. Co-extracted electrons impact on the surface of the EG, where their current can be measured. Due to their higher mass, negative ions undergo a significantly smaller, but not negligible, deflection than the electrons, and pass through the EG to be further accelerated. Electrical current measurements from the components downstream of the EG can be used to estimate the total extracted negative ion current [12].
With the combination of filter and deflection fields, the net magnetic field within the extraction systems has components in both horizontal (y) and vertical (z) directions, both of which are perpendicular to the direction of travel of the beamlet (x). These fields therefore cause deflection of the beamlets as they are accelerated, in both the horizontal and vertical directions. Figure 2 depicts this situation for the extraction system of the BATMAN Upgrade test facility (BUG) [29], which uses an RF driven IPP Prototype Source [30]. This figure shows a cross section of a single aperture, and horizontal and vertical slices through the total magnetic field within two of the apertures. Also shown is the layout of the 5×14 apertures on the PG, roughly equivalent to one beamlet group of the planned ITER source. Neither the filter nor deflection fields are uniform across the whole source, and there are differences between the fields seen by each beamlet. These differences are most noticeable after the Grounded Grid (GG) in figures 2(c) and (d), although there are changes throughout the whole volume. Beyond the GG the magnetic components in the vertical plane are completely different between the two apertures. The components in the horizontal plane show a shift in figure 2(c), which brings the net field away from being truly perpendicular with the beam direction. Although the individual differences are hard to quantify, this disparity between apertures results in each beamlet having a slightly different direction, particularly toward the edges of the grid system where non-uniformities in the field are highest. As a measure to increase the extracted ion current, caesium is evaporated into the plasma expansion region. When deposited onto the plasma facing surface of the PG, the work function of this surface is reduced, which enhances the surface production of -H (or -D ) ions [12]. As the caesium is redistributed by the presence of the plasma [13], inhomogeneities in the plasma can cause non-uniformity in the surface work function and subsequently in the efficiency of surface ion production [14].
The possible combination of plasma inhomogeneities and nonuniform surface production efficiency leads to each beamlet potentially having a different ion current and, consequently, divergence [10]. As well as posing issues for meeting homogeneity targets, this creates difficulties with diagnosing a beam comprised of potentially hundreds of beamlets, as it is typically only possible to perform measurements on collections of beamlets. Whilst diagnostics do exist that are able to resolve individual beamlets through calorimetric measurement [15,31,32], these are necessarily placed within a few cm of the extraction system, where beamlets have not yet merged, and power densities can be in excess of 10 MW m −2 [33]. Probe type diagnostics are not able to survive such high heat loads, and so they require either a modification of the ion source or a reduction of the ion current below the optimum to operate.
Options for diagnostics on large powerful ion beams are therefore somewhat limited. A diagnostic that can be implemented relatively easily on any source is BES [34]. This non-invasive diagnostic measures the Doppler shifted H α or D α emission from fast excited that are formed from beam particles interacting with the background gas. The amount of Doppler shift for photons from beam particles depends not only on their energy but also the angle of their velocity vector with the BES collection optics. By neglecting the small spread of energy of the beam particles, the width of the Doppler peak can be related to their angular distribution, and so a divergence estimated [6]. The total amount of light under the Doppler peak can also be used to measure the relative intensity of the beam within each line of sight. Light collected at wavelengths between the Doppler shifted and unshifted peaks can further give estimates of how much stripping occurs in the grid system. BES is typically measured for a narrow line of sight 1-2 cm in diameter. In BUG there are two separate arrays of BES lines of sight, crossing the beam centre either 26 cm or 1.29 m from the GG, with either 5 or 11 vertically arranged lines of sight respectively, each at an angle of 57°.
As mentioned above, each beamlet may have different properties, and may be travelling in slightly different directions. As light can only be collected from many beamlets at once, the Doppler peaks for each beamlet will be mixed into a single spectrum. The measured Doppler peak will therefore be a summation of many smaller peaks, each having a different width and central wavelength. With only the information in the measured spectrum itself, it is not possible to know if the divergence obtained through fitting of this peak is affected in some way by the beamlet mixing that occurs in the line of sight. This issue becomes particularly relevant for the very large beams planned for ITER [1], where dozens of beamlets may intersect each BES line of sight.
The other main option for these ion beams are calorimetric type diagnostics. Calorimetry can be performed on beam dump systems, such as the one 2.1 m from the GG in BUG [35]. Localised measurements of the temperature of the beam dump calorimeter can be performed using either thermocouples embedded in the surface, or by 2D infrared imaging. Combined with water calorimetry of the beam dump for the total deposited energy, these temperatures can be used to estimate the power density of the beam with some spatial resolution, typically on the order of centimetres. These measurements can give information on the size and location of the whole beam several metres downstream of the grid system, which can then be used to estimate the average beam divergence and net deflection.

Large beam simulation with BBCNI
The requirements for accurately modelling large ion beams and their diagnostics limit the available methods of simulation. The system is required to be fully 3D, in order to encompass the complex 3D electric and magnetic fields [36]. Each beamlet may experience different magnetic fields [37], as seen in figure 2, such that there are differences in beamlet formation between apertures. There may be gradients in the extracted current density across the extraction area [10], leading to further differences between beamlets. As the magnetic field and extracted current density both affect the optics of the beamlet during formation, full ion optics simulations need to be carried out for each aperture; the result of a single simulation cannot be said to apply to all the beamlets. Synthetic diagnostics of the resulting beam need to be included metres downstream of the grid system to realistically recreate the merger of many individual beamlets, as well as any additional beam transport effects due to long range fringe fields. Thus the physical size of the simulation domain, shown in figure 1, is on the order of metres, while capturing the details of the fields and beamlet formation within the apertures requires sub-millimetre resolution. Particle energies may range from 1 eV before extraction [38] to 1 MeV after acceleration [39], depending on the system.
As no ensemble assumptions about velocity or spatial distributions can be made, a self-consistent PIC simulation would be the most accurate choice, one that calculated the space-charge of the beam and interacted particles with it. However simulating such a large volume with sufficient accuracy to steady state would result in unmanageable requirements for data volumes and computation time. A particle tracking code is considered a reasonable compromise for the beamlet formation and transport, and can be supported by single aperture ion extraction calculations for solving the electric fields including the space charge contribution of each beamlet [40]. This is the form that the beam simulation of BBCNI has originally been designed with [24].
The version of BBCNI presented in this work is the result of extensive upgrades and improvements. Approximately 90% of the code base has been rewritten in a computationally efficient and modular fashion, with the key ideas kept or improved upon where appropriate. Particular attention was placed on the replacement of numerical solvers with analytical solutions where possible, which has helped the new version to achieve a two order of magnitude increase in speed. Distinct and easily replaceable modules now exist for all aspects of the simulation, which allows for easier development and more flexibility for future upgrades.
This upgraded version of BBCNI can be split into two main parts: a beam simulation including particle trajectory integration, Monte Carlo reactions, collisions with physical objects like grids and walls, and weighted photon generation; and a suite of synthetic diagnostics which analyse the resulting data. Records of generated photons are used to create simulated BES spectra for each line of sight. Particle records provide information for calorimetric and electrical diagnostics, as well as information unavailable in the experiment, such as beam neutralisation fraction and particle energy spectra.
In a full simulation H − (or D − ) ions are started in the plasma close to the PG, from where they are accelerated through pre-supplied 3D electric and magnetic field maps. This ion type also defines the species type used as the background gas and other particles; H and D may be interchanged in this document, unless explicitly mentioned. It is important to state that the electric fields are calculated for each aperture in isolation using ion optics calculations that consider the voltages applied to the grids and the current density extracted from the aperture. These calculations are performed only up to the end of the grid system. This means that although beamlet formation is modelled accurately, beamlet-beamlet interaction is not considered, whereby particle trajectories are affected slightly by the space charge of neighbouring beamlets. After the grid system, the beam is assumed to be fully space-charge compensated, whereby the creation of positive charges through interaction with the background gas counteracts the space-charge of the negative ion beam. This is considered a reasonable assumption if the background gas pressure is sufficient (around 10 −5 mbar) [41].
After defining the initial positions of around 10 5 particles for each aperture, they are tracked until hitting a surface or interacting with the background gas, which can range in pressure from 10 −3 to 1 Pa. Background gas and magnetic field information is supplied for the whole domain, so that their effect on the entire beam line is considered, including the long range residual filter field. Ions and fast atoms may interact with the background gas at any point to produce fast and slow H atoms, fast H + ions, or slow + H 2 . The reactions currently included in BBCNI and the sources of their cross sections are given in table 1. Reactions 7 and 8 are used only to calculate weights for the photons generated by fast -H or H. Unfortunately the cross sections available for these reactions are not well defined, with uncertainties of over 10% in many cases [42,43]. Particularly problematic is the need to extrapolate some cross sections to lower energy collisions below a few hundred eV due to lack of experimental data. Even less data is available for explicit D reactions, and so the same cross sections are used for both isotopes.
Reactions that occur within the grid system have a special significance for two reasons. Firstly, if a H − reacts to create a fast H atom and an e − , it is termed 'stripped' and the resulting H atom has less than the full energy of the beam. The resulting electron is further accelerated, and can go on to damage beamline components. There is therefore a maximum allowable amount of stripping, which depends on the system, but for the ITER systems is estimated at 30% [1]. This value depends not only on the beam energy, but also the background gas density within the grid system, and therefore the source pressure, the geometry of the grids and any temperature gradients that are present. The second significant process is the creation of positive ions within the grid system, which can be accelerated in the opposite direction to the rest of the beam, and go on to hit the backplate of the source wall. These 'backstreaming' ions are significant for their potential to cause damage to the source walls [8,45] as well as their contribution to caesium redistribution from the backplate [14,46]. As well these reactions, the generation of photons is considered for the synthetic BES diagnostic. In reality, both H − ions and fast H atoms interact with the background gas to form excited atoms, which can then decay through spontaneous emission to emit H α photons at 656.279 nm [47] (D α at 656.101 nm [47]). However, in order for sufficient photon generation within the collection volume of the BES optics via Monte Carlo reactions of the correct type, particle numbers would be two or more orders of magnitude higher than current simulations, rendering computation unfeasible for now. Therefore both H − ions and H atoms produce weighted photons corresponding to the fast and slow particles without the need for a specific reaction to occur. The weight of each photon is interpreted as an intensity value, and is calculated from the gas density at the particle location and the cross section of the relevant reaction from table 1.
The important modules used in BBCNI and their information connections are laid out in figure 3, with the main inputs on the left, and outputs on the right. As stated above, the electric and magnetic fields through which particles are tracked need to be supplied from external calculations. Electric fields are currently supplied by calculations performed using the IBSimu code [40], which simulates the area around a single aperture. Magnetic fields are supplied using either the PerMag code [48] for permanent magnets, or ANSYS ® Workbench TM for fields generated by flowing DC currents through the plasma grid. Both of these codes are able to provide magnetic field values on a 3D mesh of practically any size. Design drawings and specifications are used to provide information about the source geometry and diagnostics of the relevant experiment. Gas densities are estimated from pressure and temperature profiles generated by an inhouse gas flow simulation. Using equations for gas conductance through a cylinder and estimates of gas temperatures at each grid, 1D profiles of neutral gas pressure and temperature are generated along the centre of an aperture. This is then applied to all apertures of the system.
The input modules allow for a great amount of flexibility in the provided data, and in principle any ion source can be modelled. For each type of field, the user can specify a number of different files, each with a 'priority', that are all considered part of the same field. The data in these files does not have to match in resolution or the volume of space they describe. Any spatial overlap between the information in these files is handled by the priority system. This is necessary as, for example, the filter field within the extraction system can have features on the sub-mm scale. Including such detail for the whole domain would be unnecessary and highly impractical, as variations in the long range field are typically important only over of tens of cm. As the filter and deflection fields for ion sources are often generated independently, data for these fields can be provided separately. Due to the possible differences between apertures, as discussed above, individual 3D electric field maps can be provided for each aperture. The fields are stored and interrogated separately, with linear interpolation used on each vector component individually. Data files are assigned a priority to resolve any overlapping that may occur, to provide a single result for the vectors of E and B at any point in space.
Particles are given initial conditions sampled from userspecified spatial and velocity distributions. Currently ions start on a disk at a user-defined distance from the PG, as an approximation of volume produced negative ions, but there are future plans to incorporate an estimation of surface production using a defined flux from the surface of the PG. Particle initial velocities are defined by a uniform initial energy, to which parallel and perpendicular temperatures can be added. Particles are moved through the E and B fields by a Boris-style [49] particle integrator with both adaptive step size and full relativistic treatment [50]. After each integration step, the trajectory is checked for collisions with physical objects in the beamline-grids, walls, or the calorimeter-and with non-destructive diagnostic planes, for the recording of all particles passing through a particular 2D plane. If one occurs, the particle trajectory is interpolated to the collision point and these properties recorded.
Following this, a Monte Carlo routine checks if a reaction with the background gas occurs. If a reaction occurs, it is recorded, and the trajectories of the resulting particles integrated as before. Slow particles created outside the grid system can take an unreasonable number of integration steps to collide with an object, so the option exists to kill those below a certain energy immediately after creation to save computation time. As part of the reactions module, both H − and H atoms produce weighted photons that are recorded if they are inside the acceptance cone of one of the simulated BES lines Necessary input data is shown on the left, processing modules in the centre, and outputs on the right. A distinction is made between the raw records of particle locations and properties, and the results of post-processing these records. Also shown is the option to record full particle trajectories. of sight. The simulation is also able to record entire particle trajectories, although this comes with a significant increase in computation time and disk space, generating up to 1 GB per aperture for just 5k initial particles. The particle tracking, collisions, and reactions are parallelised aperture-by-aperture using OpenMP, for a simple yet robust scalability.
Once the integration of all particles and reaction products is completed, the records of particle destinations, generated photons, and reactions are passed to the post-processor for simulating the diagnostics. A specific module exists for the processing of photon records into beam emission spectra as would be seen by a BES diagnostic. Each photon record is treated as a point source radiating uniformly in 4π sr, with an intensity given by the weight calculated during the particle integration. For each source that lies within the viewing cone of a BES line of sight, the entrance of the collection optics is sampled randomly to generate a collection of incoming photons to be considered; currently 1000 sample points are used per photon source, but this is user-defined. For each of these photons, the relativistic Doppler shift is calculated using the velocity of the photon source and the vector between the source and the collection optics. The intensity of the photon is scaled (on top of the weighting described earlier) according to the distance between the source and the collection optics, as well as due to relativistic beaming [51]. This scaled intensity is then added to a spectral histogram according to the Doppler shifted wavelength. This process is also parallelised by aperture, such that individual spectra are generated for each. The full simulated spectrum is the sum of these individual contributions, however being able to trace particular spectral observations to specific apertures is of great advantage to understanding diagnostic results.
Further, simpler, statistics are produced for the total energy, charge, and particle numbers deposited onto each object in the simulation domain. These values are resolved by particle type, and allow direct insights into global measures such as ion neutralisation fraction and grid currents. The particle destination records are used to create 2D histograms of the surface density of energy deposition onto the calorimeter and grid system. These can be used for the detection of hotspots and the validation of the code with experimental results. Energy spectra for particle populations are generated, for the analysis of backstreaming positive ions and stripped particles.

Output capabilities
In order to demonstrate some of the capabilities of BBCNI, the RF driven prototype source used at BUG was simulated. The parameters were chosen to match measured values for a particular experimental shot undertaken during the conditioning phase of operation. The extraction potential was 5.7 kV and the acceleration potential 25 kV. The magnetic fields were as shown in figure 2, with = I 2.2 kA PG , which corresponds to a field strength of around 4.6 mT just in front of the PG. The source pressure was 0.4 Pa, reducing to 2.5×10 −4 Pa in the tank. A current density of 12.8 mA cm −2 was assumed to be extracted homogeneously from the 5×14 apertures of the grid system. 10 5 initial particles were used per aperture, for a total of 70 million initial particles. Using 16 cores of a cluster based on Intel ® Xeon ® E5-2670 CPUs, the simulation ran for a total of 5 h40 min, of which around 30 min was spent in post processing. The particle records alone sum to roughly 14 GB of data files, which is why pre-processed statistical analyses are used.
One of the simplest analyses of the results is shown in figure 4, which shows the beam power density in various locations. The 2D map at the diagnostic calorimeter given in figure 4(a) shows clearly the downward deflection of the beam as a result of the fringe filter field present in the grid system, as well as the elongated shape due to the 5×14 aperture pattern. As this calorimeter is just over 2 m from the grid system, the beam needs a downward angle of only 2°to give the amount of shift seen. The merging from the initial 14 beamlet rows into a single beam can also be seen in both parts of figure 4, and is shown to occur within 50 cm of the last grid. For these system parameters the divergence of the beamlets is roughly 2°; BUG has demonstrated divergences of less than 1°, which would lead to a slower merging of the beamlets, as well as a smaller overall width for the 2D profile observed in figure 4(a). Figure 5 shows examples of both experimental and simulated BES spectra from a vertically central line of sight, which crosses the centre of the beam 1.29 m from the GG at an angle of 57°. The locations of the peaks in the simulated data are well in agreement with the experiment, although there is a slight mismatch in the peak widths. After generation of the simulated spectrum the data was convoluted with a Gaussian of σ=12 pm to simulate the apparatus profile of the experimental system. This indicates that other broadening mechanisms are the cause of the large peak widths observed in the experiment. The whole-beam divergence as evaluated from the simulated data is 2.1°, compared with 2.7°for the experimental results, due to the difference in peak widths. The wider unshifted peak in the experimental data suggests that the temperature of the emitting particles is higher than the 300 K assumed for the background gas. As the background gas pressure in this example is a relatively low 0.4 Pa, the stripping peak is significantly less intense than the other two; although it can be seen in the simulated data with a change of scale, it falls below the noise floor in the experimental spectra.
As discussed in sections 1 and 2, there are a number of possible reasons for the discrepancy between the experimentally measured divergence and the simulated one, due to either properties of the diagnostic system or the beam itself that are not considered in the code. Effects such as the influence of surface produced ions, or the impact of inhomogeneous extraction, may lead to differences not just in the divergence of each beamlet, but also non-Gaussian velocity distributions, and can change the direction in which the beamlets propagate, due to interaction with the magnetic fields. These types of effects cannot be measured directly in the experiment, but can be in the simulation. Not only can BBCNI determine the potential impacts on diagnostic results, but can also give information on each beamlet individually.
An example of beamlet level information is given in figure 6, which shows the fraction of total intensity contributed by each aperture to the Doppler peak observed by the synthetic BES observing 1.29 m from the GG. The approximate acceptance cone for the vertically centred collection optics is sketched over the figure for reference.
Two effects are seen in this data, beyond what one would trivially expect. Firstly there is a horizontal difference in the contributions across each row. This has a relatively trivial explanation in that the apertures to the left of the figure are closer to the collection optics. The line of sight then subtends a larger solid angle of the simulated light emission from these beamlets. The second observable effect is the downward shift due to the filter field, as also seen in figure 4. This is apparent here in that the highest contributions come from apertures above the collection optics. As the line of sight is observing a point 1.29 m from the grid system, the net downward shift is smaller than that observed at the calorimeter.
If the net downward shift is analysed for simulated data from both the BES (1.29 m from the GG) and calorimeter (2.1 m from the GG), one finds that the shift does not increase linearly with the distance-the net beam angle must change to some degree after the grid system. This is shown to be the case in figure 7, which gives more detailed results from the BBCNI calculation. The bottom half of figure 7 shows a sampling of particle trajectories projected onto a vertical plane through the centre of the beam. Individual particle trajectories show how much overlap there is between beamlets, while average beamlet centrelines show that each beamlet is affected approximately equally by the long range magnetic field. The mean vertical beam angle is shown in the top section of figure 7, the differential of which follows exactly the horizontal field strength, also plotted. It should be noted again that beamlet-beamlet interaction is not considered in these results. It would be expected that, with its . Shown is (a) a 2D simulated calorimeter power density profile 2.1 m after the GG with a total power of 40.4 kW and (b) the evolution of the vertical power density profile of the beam as a function of the distance from the grid system. Due to a high divergence of around 2°for these parameters, all 14 beamlet rows have merged into a single beam after around 50 cm. The vertical deflection of the beam due to the filter field is also apparent. Figure 5. Example of BES spectra from the BBCNI simulation (line) and the BUG Experiment (boxes). The important peaks from various particle populations are labelled, giving the Doppler peak from full energy excited particles, the stripping peak from partially accelerated particles, and the unshifted peak from excited background gas. inclusion, the average beam angle would remain the same, but the top-and bottom-most beamlets would be deflected slightly outward, broadening the beam.
The results presented so far show how performing diagnostics on large ion beams presents challenges in terms of analysing the resulting data. In particular, the differences between the properties of individual beamlets, and emergent properties of the beam, is not always clear, particularly if the beamlets themselves are not Gaussian in nature. The divergences calculated from the spectra shown in figure 5 assume that the beam can be treated as a single entity, when in reality each beamlet can have unique properties. This is not likely to be an issue in smaller ion sources that use just a handful of apertures, or even just one. However in larger sources, particularly those on the scale of metres in width as planned for ITER, differences in beamlet properties occur more gradually over the grid system, and it is not possible to disentangle results from individual beamlets. This information is readily available in the simulation, and so these effects and their impact on the experimental results and trends can be investigated. The accurate recreation of diagnostic data relies not just on what happens during beamlet formation, but also on the transport of the beamlets to the location of the diagnostics, and the consequent beamlet merging that occurs. The consideration of these processes is shown by BBCNI to give good recreations of experimental diagnostic data, which can be used to further understand the generation and transport of ITER-relevant large negative ion beams.

Conclusions
The results in this work show the importance of a full ion beam simulation that includes all the beamlets individually; it is not sufficient to obtain results from models of single beamlets and project to the position of interaction with the diagnostics, nor is it sufficient to rely on whole-beam results as representative of the individual beamlets.
Experimental diagnostics that are available on such systems do not provide enough information to disentangle data acquired from multiple beamlets, which even under homogeneous extraction conditions, will be travelling in different directions due to the magnetic fields present. Thus the required distinction between beamlet and whole-beam level information is not available using experiments alone.
In order to give the required traceability and allow beamlet level information to be obtained, the code BBCNI has been developed to provide a full-beam simulation including interaction with synthetic beamline diagnostics. These synthetic diagnostics have the important addition of traceability back to individual apertures or particles. This allows the investigation of the impacts on diagnostic results of beam inhomogeneity or of changes occurring during beamlet formation.
The successful recreation of beam emission spectra and their key features, and the ability to trace the contributions that each aperture makes to such a spectrum shows the value of BBCNI. The code is planned to be used both in conjunction with experiments for enhanced data analysis and as a stand-alone tool for fundamental investigations. There is also scope for use as a predictive tool, particularly for assisting the design of new beam diagnostics, or investigating the impact of different filter field configurations on the beam.
Despite simulation results comparing well with experiment across a range of parameters [52], there are still improvements to make and considerations to include in the code that will greatly increase the applicability. Most importantly, inhomogeneous extracted current caused by magnetic drifts, possibly on both whole-beam and beamlet scales, and thereby gradients in the beamlet divergence, has been shown to have a significant impact on the beam properties [18,25,53], and is therefore the next planned upgrade to BBCNI.