Present status of numerical modeling of hydrogen negative ion source plasmas and its comparison with experiments: Japanese activities and their collaboration with experimental groups

The present status of kinetic modeling of particle dynamics in hydrogen negative ion (H−) source plasmas and their comparisons with experiments are reviewed and discussed with some new results. The main focus is placed on the following topics, which are important for the research and development of H− sources for intense and high-quality H− ion beams: (i) effects of non-equilibrium features of electron energy distribution function on volume and surface H− production, (ii) the origin of the spatial non-uniformity in giant multi-cusp arc-discharge H− sources, (iii) capacitive to inductive (E to H) mode transition in radio frequency-inductively coupled plasma H− sources and (iv) extraction physics of H− ions and beam optics, especially the present understanding of the meniscus formation in strongly electronegative plasmas (so-called ion–ion plasmas) and its effect on beam optics. For these topics, mainly Japanese modeling activities, and their domestic and international collaborations with experimental studies, are introduced with some examples showing how models have been improved and to what extent the modeling studies can presently contribute to improving the source performance. Close collaboration between experimental and modeling activities is indispensable for the validation/improvement of the modeling and its contribution to the source design/development.

The H − production mechanisms in such H − plasma sources are classified into mainly two mechanisms, i.e., (i) volume production (VP) [11], and (ii) surface production (SP) [12]. In the former, the H − VP is mainly conducted by the following two-step processes, (1) the production of vibrationally excited H 2 molecules by energetic electrons (E>20 eV) (EV process), and (2) H − production by slow electrons (E1 eV) (DA process). Therefore, the control of the electron energy distribution function (EEDF) is very important for efficient H − VP. In the latter mechanism (H − SP), H − ions are produced by hydrogen neutrals and positive ions which impinge on a low work function surface, e.g., Cs covered metal surface.
Various types of H − ion sources have been developed, depending on the purpose and application [1]. Among them, we focus on the multi-cusp arc-discharge source and radio frequency inductively coupled plasma (RF-ICP [13]) source. In the former type of H − sources, DC arc-discharge plasmas are produced between the Not only understanding the H − source plasma characteristics above, but also understanding the H − extraction process from the H − sources is one of the most common and critical issues of H − source development. The understanding of H − extraction physics is very important in helping improve the extracted H − beam quality. As has been reviewed in [16], modeling studies with electrostatic particle in cell (ES-PIC) are useful in helping understand the electric potential structure near the extraction aperture and to clarify the process leading to the extraction of H − ion from the source.
A series of 2D3V ES-PIC (two-dimensional in real space and three-dimensional velocity space, ES-PIC) modeling has been done in [36][37][38][39][40][41][42]. A reduced-size-scaling has been adopted to save the computational cost by keeping important normalized parameters the same as those in the real system [36,[38][39][40][41][42], e.g., normalized Larmor radius r L, j Lj * r = / collisionality parameter L, jk jk * l l = / where L is the characteristic scale length of the real system, while r Lj and jk l are the Larmor radius and collision mean free path, respectively. The symbols j and k denote particle species. Recently, almost the same 2D3V ES-PIC modelings [32,43], but with the different scaling (i.e., reduced-density scaling) have been done with a more detailed model of the collision processes with Monte-Carlo method for collisions (MCC).
The invariant property between the real system and the numerical scaled system has been discussed in [32]. Above, different scaling schemes (reduced-size scaling and reduced-density scaling) have been compared with each other and they are shown to be almost equivalent, if a certain relationship is satisfied for the size and density scaling factor. As mentioned above, j * r and jk * l are kept the same as in ) at least for a simple 1D problem. This property has been pointed out by the kinetic study in [44] by numerically solving the plasma-sheath equation derived by the Vlasov equation and Poisson equation.
The simple ES-PIC model in [36][37][38][39][40][41][42] mentioned above has been improved step-by-step through the model validation by comparisons with various experiments and other theoretical/numerical results [44][45][46]. The models give useful insights into the H − extraction process, and also help identify basic and intrinsic mechanisms of various interesting phenomena in the qualitative manner. For example, these 2D ES-PIC modeling studies have given a possible answer to the following questions: (electrons are magnetized, while H + and H − ions are not magnetized) in the extraction region near the single extraction aperture. Why? [37] This trend has been observed experimentally in the Camembert III H − source under the Cs-free H − VP operation [47].
(2) The extracted H − current I H − has a peak with respect to the applied plasma grid (PG) bias voltage. Why there exists an optimum PG bias voltage for the extracted H − current I H −? [38] This trend has been observed in the Camembert III H − source under the Cs-free H − VP operation [48].
(3) The beam halo component has been observed in the giant H − source experiments under the Cs seeded H − SP operation [41,42,49,50]. Why?
(4) Why are so-called 'ion-ion plasmas' formed in front of the extraction aperture and the PG under Cs-seeded H − SP? [39-42, 51, 52] (5) From where are most H − ions extracted under Cs-seeded H − SP operation? Directly from the surface of the PG or through the bulk plasma?
Recent experimental observations in the NIFS (National Institute for Fusion Science) R&D (Research and Development) H − source for fusion plasma heating [51,52] raised questions (4) and (5). Interesting and important points in these experiments are summarized as follows: (i) the ion-ion plasma layer which consists of only H + and H − ions almost without electrons is extended over the extraction region with a relatively large thickness (∼10-20 mm), (ii) extraction grid (EG) voltage affects H − ions existing far away (∼10-20 mm) from the aperture.
Modeling studies in [36][37][38][39][40][41][42] give useful insights in helping to solve these questions, but problems (4) and (5) above have not yet been completely understood. Understanding of the ion-ion plasma and also plasma meniscus (H − emitting surface) shape in such ion-ion plasma are very important in helping to obtain good beam optics and to reduce the beam halo component.
As mentioned above, these 2D ES-PIC modeling studies are quite useful for the basic understanding of H − extraction physics, but are still far from the perfect. For example, the extraction aperture has been modeled as an infinite slit instead of the hole in these 2D3V ES-PIC modeling. Recently, full 3D ES-PIC models are being developed with the realistic dimension without the size scaling above in [53][54][55]. These full 3D ES-PIC models are still on the way towards further model improvements.
Recently, the 2D ES-PIC model in [36][37][38][39][40][41][42] above has been extended to the 3D ES-PIC model and simulations of the ion-ion plasma have been carried out [56]. The 2D ES-PIC results of the formation mechanism of the ion-ion plasma have been validated by the 3D ES-PIC model. The results will be reviewed in section 4.1. In addition, a systematic comparison of the results has been made between 2D ES-PIC and 3D ES-PIC modeling in section 4.2.
Furthermore, based on the 3D ES-PIC model above, an integrated 3D model from the plasma meniscus formation to the beam acceleration [57] is in development for the purpose of more quantitative comparisons with the experimental results of H − beam quality, such as the fraction of the beam halo component to the core beam component, beam emittance and so on. These results will be discussed in section 4.3.
Finally, in section 5, the present status of the model development and the understanding of the H − source plasmas and beam extraction have been summarized together with the open problems.

Modeling of multi-cusp arc-discharge sources
In this section, modeling studies of multi-cusp arc-discharge H − ion sources and their comparisons with the experiments are reviewed with some new results. In order to understand the basic physics and to improve/ optimize the source performance, the KEIO-MARC code (Kinetic modeling of Electrons in the IOn source plasmas by the Multi-cusp ARC-discharge) [26][27][28][29][30] has been developed. First, we briefly summarize the model implemented into the KEIO-MARC code in section 2.1. Then, in section 2.2, the code application to the compact H − source for medical application will be presented. The focus is placed on H − VP in the magnetic filter configuration. The effects of the non-equilibrium feature of the EEDF will be discussed on the vibrational excitation of H 2 molecules and the resultant H − VP. Next, in section 2.3, the so-called 'giant' H − sources for the plasma heating in magnetic fusion application will be discussed. Much attention will be placed on the Cs seeded H − SP. The plasma spatial non-uniformity/asymmetry experimentally observed in such giant sources will be discussed with the modeling results by the KEIO-MARC code. [26][27][28][29][30] The KEIO-MARC code is based on a test particle Monte-Carlo model in 3D real and velocity space (3D3V). Realistic 3D geometry and 3D magnetic field configuration of multi-cusp arc-discharge sources can be taken into account.

Brief summary of the physics model implemented into the KEIO-MARC code
Under the assumption of quasi-neutrality, we focus on the electron dynamics in the arc-discharge source. As will be discussed below, the EEDF is crucially important for H − production in H − sources. In order to obtain the EEDF, the KEIO with the Boris-Benemann version of the leap-frog method [15], where m , e v , e q, E and B are the electron mass, velocity, charge, electric field and magnetic field, respectively. The last term on the right-hand side of equation (1) expresses the momentum change by collision processes. Debye length becomes 10 10 m l~---for the typical plasma parameter range (electron density n e ∼10 16 -10 18 m −3 , electron temperature T e ∼1-5 eV) of H − multi-cusp arc-discharge sources. The quasi-neutrality condition may be valid in most of the source volume except for the thin sheath layer (a few Debye length) near the anode wall and the cathode filaments. The local electric field may be neglected in equation (1) under the assumption of the charge neutrality. This assumption is rather crude, therefore careful model validations are needed and we have made many comparisons with experiments, not in one discharge/source condition, but in various discharge/source conditions up to now. These results will be shown later in sections 2.2-2.3.
The sheath potential drop is also not exactly calculated, but roughly estimated from the following simple analytic formula [44] (equation (2) below) by assuming the current density at the arc-discharge wall is small due to the large wall surface area, where k , B e(= q | |), and m i are the Boltzmann constant, unit charge, and ion mass, respectively. In equation (2), T e is estimated from the average kinetic energy of the test electrons in the cell in front of the wall. When test electrons reach the chamber wall, the electrons with kinetic energy K e larger than the sheath potential drop K qV e sh > ( | |) are absorbed at the wall, while those with low energy K qV e sh < ( | |) are reflected. In the KEIO-MARC code, collision processes are modeled by Monte-Carlo techniques. They are mainly classified into two categories. The first one is the collision processes of electrons with hydrogen particles (H atoms, H 2 molecules, H + , H , 2 + and H 3 + ions). Totally, about 540 collision processes can be included with tabulated data for the collision cross-section. Among these processes, the main collision species are summarized in table 1. Also, cross-sections of the reactions No. 1-5 are shown in figure 1.
The second category is the electron-electron collision, and electron-ion collision, i.e., Coulomb collision. The collisions in the former category with hydrogen particles are modeled by the 'null-collision method' (NCM) [62], while the 'binary collision model' [63] is applied to the latter, i.e., Coulomb collision.
On the basis of the model explained above and using the main input parameters listed below, the KEIO-MARC code calculates the EEDF in each numerical cell in the discharge chamber: (1) 3D geometry of the H − ion source shape and dimensions of the source chamber, locations of the cathode filaments, (2) 3D configuration of the magnetic field, (3) discharge conditions (arc power P arc or arc current I arc , arc voltage V arc and gas pressure P H2 ). (4) Background particle (H, H 2 , H + , H , 2 + and H 3 + ) density and temperature, and (5) numerical parameters (such as spatial-cell volume V , D time step t, D and the weight of Monte-Carlo test particle n w ).
The initial energy of the test electrons emitted from the filaments is set to be V arc by assuming the cathodesheath acceleration of electrons after passing through the thin sheath layer. The number of the test electrons emitted from the filaments per each time step is determined by the arc current I arc (=P arc /V arc ). In the simulation, not only the primary electrons from the filaments, but also all the trajectories of the secondary electron produced by ionization inside the discharge volume are followed after their production. The time step is updated until the total number of test electrons N T in the system (i.e., inside the arc chamber) reaches a steady state. After reaching the steady state, the EEDF in each cell, and electron density and electron temperature, and other macroscopic physical quantities can be calculated with the EEDF.  The model implemented in the KEIO-MARC code has been validated by comparisons of the modeling results with those in the experiments. It should be noted that the validations have been done not in one discharge/source condition, but multi-source comparisons have been done up to now. In other words, the KEIO-MARC code has been already applied to many arc-discharge sources in different applications, e.g., QST (National Institute for Quantum and Radiological Science and Technology) 10A Ion source [21,25], JT-60SA Kamaboko source [22,23], and NIFS (National Institute for Fusion Science) 1/2 R&D (Research and Development) source [51] for the plasma heating in magnetic fusion devices, and SHI (Sumitomo Heavy Industry) compact source [17,18] for the medical application as mentioned in section 1.
As discussed below, the non-equilibrium features of the EEDF are very important for the both volume and surface H − production. [17,18,64,65] Here, we focus on the modeling of Cs-free H − VP. As mentioned in section 1, the H − VP is mainly through the following two step processes, (1) the production of vibrationally excited H 2 (v) molecules (EV process) by fast energetic electrons (E>20 eV), and (2) H − production (DA process) by slow electrons (E<1 eV). The control of the EEDF inside the chamber by the magnetic filter concept is essential for efficient H − VP. In order to enhance/optimize the extracted H − beam current, it is necessary to understand the relationship between the H − production and the design/operating parameters of the H − source, such as arc-discharge power, magnetic field configuration, and so on. For this purpose, a systematic study by the KEIO-MARC code and comparisons with the experimental results have been done for the SHI H − source.

Compact H − sources for medical applications
In the following, we will discuss mainly the filter effect on the EEDF calculated by the KEIO-MARC code. Also, using the EEDF obtained by KEIO-MARC code, we have made simple evaluations of n H2(v) (the density of the vibrationally excited molecules) and n H − (the H − density) produced in the source. The main focus is placed on the dependence of the filter effect and these quantities (n H2 (v) and n H −) on the arc voltage V arc in this paper. The dependence on other parameters (e.g. arc power P arc , arc current I arc , magnetic filter field strength, B filter ) has been also studied and discussed in detail in [18,64,65]. Figure 1 shows a schematic drawing of the SHI H − source. The source consists of a cylindrical plasma chamber made of copper, multi-cusp permanent magnets, and three electrodes system. The SHI source has only a single extraction aperture. Hydrogen plasma is generated by the arc-discharge between two tungsten filaments and the chamber wall. The internal diameter and length of the chamber are 100 mm and 160 mm, respectively. The source is compact in comparison with so-called giant H − sources for the magnetic fusion application discussed in section 2.3.

SHI H − source
The filter magnets are in the extraction region of the chamber to create the magnetic filter field whose direction is perpendicular to the central axis of the chamber. The peak magnetic field in the extraction region has been ∼110 Gauss. As mentioned above, the magnetic filter divides the chamber volume into 'a high energy electron region (driver region)' for the EV reaction (see table 2) and a 'low energy electron region (extraction region)' for the H − VP (DA reaction: see table 2).  [58] Note. eV: vibrational excitation of molecules by electrons, EV: vibrational excitation of molecules by energetic electrons through electronic excitation and subsequent radiative decay, Diss: dissociation, DA: dissociative attachment, IC: ion conversion, ED: electron detachment, MN: mutual neutralization, AD: associative electron detachment.

EEDF analysis by the KEIO-MARC code
Taking into account the real 3D geometry and 3D magnetic configuration of the SHI H − ion source in figure 2, a systematic parameter survey for the dependence of the EEDF on the discharge parameters (P arc , I arc , and V arc ) has been done. The hydrogen gas pressure has been set to be 1.6 Pa, which is estimated from the actual value of gas flow rate measured experimentally. The internal space of the chamber is divided into small cubic cells (2 mm×2 mm×2 mm). The wall boundary is specified by the cell located at the chamber wall. The background gas parameters are roughly estimated from [25,[66][67][68] and fixed (n 1.5 10 m , ). In the early phase of the code development, the background ion parameters are assumed to be fixed (n 4 10 m , ), but the background ion density can be calculated iteratively consistent with the electron density at present. Only the ion densities of H 2 + and H 3 + (or their fraction) are fixed. The results of the EEDF presented here are not so much dependent on these assumptions at least in the qualitative sense. Figure 3 shows typical 2D profiles of the electron density on the cross-sectional area perpendicular to the Zaxis. Figure 4 compares the EEDF dependence on the arc voltage V arc (a) in the driver region at Z=80 mm and (b) in the extraction region at Z=16 mm. Both results have been calculated with fixed arc current I arc (=20 A) for simplicity in order to understand the pure effect of the V arc . It should be noted that the maximum value of the vertical axis in figure 4(b) is one order smaller than that in figure 4(a).
The effect of the magnetic filter can be clearly seen in figure 4. High energy electrons are significantly reduced by the magnetic filter field in the extraction region. The electron temperatures (T e ∼2 eV) of the low energy Maxwellian part in the extraction region in figure 4(b) is reduced compared with those (T e ∼4.5 eV) in the driver region in figure 4(a). Although T e in the extraction region seems to be slightly higher than the  experimental observation by the probe measurement, the filter effect is reasonably reproduced by the KEIO-MARC code in 3D realistic magnetic configuration.
As seen from figure 4(a), the high energy tail part (E>20 eV) exists in the EEDF in the driver region. The EEDF consists of mainly two components, i.e., the low energy thermal component and the high energy tail component. This characteristic feature in the driver region has been usually observed in the experiments [25,28,68].
Characteristic feature of the EEDF shown in figure 4 is mainly explained by the energy dependence of the collision processes. As mentioned above, electrons emitted from the filament with the relatively high initial kinetic energy E e V The high energy electrons lose their energy by these inelastic collision processes step by step after their emission from the filaments. If their energies become less than E 15 20 eV , th » -( )then they hardly undergo the inelastic collision process because of the energy threshold E th of the inelastic collision.
Instead, electron-electron Coulomb collisions become the dominant collision species, because the Coulomb collision does not have the threshold and the cross-section of the Coulomb collision (i.e. Rutherford cross-section) is proportional to E . 2 -Therefore, Coulomb collisions are more effective for the low energy electrons concerning their energy relaxation process. As a result, the low energy part of the EEDF relaxes towards an equilibrium energy distribution (Maxwellian distribution) due to energy exchange by electron-electron Coulomb collisions mainly among the low energy electrons.
As will be discussed later, the high tail component (E>20 eV) of the EEDF in figure 4(a) plays a key role in producing vibrationally excited molecules H 2 (v) with a high v-number, which are in turn important for H − VP. The tail becomes longer, as V arc increases. Therefore, it is favorable to increase V arc from the viewpoint of enhancing EV reaction (table 2) and the production of H 2 (v) with high v in the driver region.
However, the electron temperature T e of the low energy Maxwellian part in the extraction region in figure 4(b) becomes larger, as V arc increases. In addition, even a high energy tail still exists in the extraction region. Especially, it is visible for the cases of higher V arc with V arc >80 V. This relatively high energy component could contribute to the destruction of H − ions. This may lead to the decrease in the net H − density.
These results suggest that the increase in V arc and the resultant increase in the vibrationally excited molecules H 2 (v) by energetic electrons (EV reaction: see table 2) in the driver region does not always lead to monotonic increase in H − production in the extraction region. In fact, as will be discussed later, the extraction H − beam current tends to be saturated with V arc increasing in the experiment.  Therefore, the H 2 (v) population with respect to the vibrational excited number density is very important. It is estimated by the following system of the rate equations for each vibrationally excited state (0v14), where n e and n H+ are the electron and ion density, respectively. The densities of hydrogen molecules at different vibrationally excitation levels v and v′ are denoted by n v and n v′ , respectively. The symbols τ and γ represent the confinement time of H 2 (v) and the reflection rate at the wall boundary. The parameter γ is taken to be 0.5 for simplicity. And the parameter τ is also simply determined by L/v th,H2 (v), where L is the radius of the plasma chamber and v th,H2 (v) is the thermal speed of the H 2 (v). In addition, R denotes the reaction rate coefficient for each reaction process. The symbol (e.g. EV) in the parenthesis for R stands for reaction species listed in where, σ reac (E) and v e (E) represent the cross-section for each reaction and the velocity of the electron, respectively. As noted above, the non-equilibrium feature of the EEDF shown in figure 4(a) in the driver region significantly contributes to producing highly excited molecules with large v (v 7 8  -). In other words, the high energy component of the EEDF is indispensable for efficient production of H 2 (v) with high v, because the crosssections of the EV reaction (see table 2) become larger for the energy range (E>20 eV) and has a peak around E∼50-60 eV.
In addition to the basic study in [64] above, a more wide range of parameter survey for V arc has recently been done with a stronger filter field (the peak strength B=150 G) to prevent higher energy electrons from crossing filter field. By taking V arc as a parameter, the calculated results of H 2 (v) density from equation with v being higher than v=4, respectively. The H 2 (v) densities are clearly saturated for larger V arc with V 80 V. arc  This can be explained mainly by the energy dependence of the EV cross-section. The cross-section of the EV reaction is relatively small in the higher energy range (E is large than 100 eV). As a result, those electrons with such a high energy cannot contribute to the production of H 2 (v) effectively. This is one of the possible reasons for the saturation of the H − production for higher V arc discussed below. Figure 6 shows the dependence of the H − density in the extraction region on V arc . The H − densities are roughly estimated by the following simple 0D particle balance equation for the H − ions:

H − production
where the numerator is the H − production rate by the DA reaction between the vibrational excited molecules and electrons, while the first three terms in the denominator are related to the reaction loss due to the mutual neutralization (MN), associative detachment (AD), and electron detachment (ED), respectively, in table 2. Their reaction rate coefficients are denoted by R. The last term is related to the H − transport loss and H tis the confinement time of H − ions in the chamber.
In equation (5), the H 2 (v) density n v is estimated by equation (3). We assume that n v is uniform in the chamber, because H 2 (v) molecules are free from the magnetic filter and the chamber is relatively compact. In addition, the electron density n e and the positive ion density n H + in the extraction region in equation (5) are estimated by the KEIO-MARC code. Furthermore, we have used the same value of neutral atom density n H as in the EEDF analysis above by the KEIO-MARC code (see subsection 'EEDF analysis by KEIO-MARC code' above).
The characteristic feature in figure 6 is possibly explained by the balance between the H − production and the H − destruction. As for the production, the production rate increases with V arc for the smaller V arc range, while it is saturated for larger V arc range because of the saturation of the H 2 (v) density n v for larger V arc as has been already discussed above (see figure 5).
As for the H − destruction process in equation (5), the H − destruction by the ED reaction in table 2 is relatively important under the condition of the present study, especially for the high V arc regime. The electron temperature in the extraction region becomes higher with V arc increasing, and the amount of the high energy part of the EEDF increases in the high V arc regime (see figure 4(b)) even in the extraction region. As a result, the H − destruction by ED reaction becomes important and increases with V arc . This increase in the H − destruction and the saturation of the EV reaction are the possible reasons for the dependence of H − density on V arc and the roll-over in the high V arc regime. In figure 6, experimental results of the extracted H − current are also shown. They tend to be saturated for higher V arc in the experiments as has been predicted by the KEIO-MARC results. From the KAIO-MARC modeling, it has been suggested that further increase in the extracted H − current cannot be expected by increasing V arc .
In this section, we have discussed vibrational excitation of H 2 (v) and the H − VP with the relatively simple 0D model for H 2 and H − to illustrate that the non-equilibrium feature of the EEDF is very important to understand the H − VP process. For more quantitative comparison with the experiments and also for a more design oriented study, transport analysis of the H 2 and H neutrals is also very important. Especially, a detailed transport analysis of H 2 (v) from the driver region to the extraction region becomes important for H − VP [71,72]. As pointed out in [73], it is necessary to deal with molecules in the different vibrational state with v as a different species in such transport analysis, because the time scale of the vibrational excitation possibly becomes comparable to the transport time scale of a H 2 molecule.

Giant H − sources for fusion plasma heating
In this section, we focus on the so-called giant H − sources for plasma heating in magnetic fusion devices. As mentioned in section 1, one of the critical issues for such giant sources is the spatial non-uniformity of the H − SP on the large extraction area. The non-uniformity leads to excessive heat power loadings onto the acceleration grids and breakdown, which limits the injection power of the beam, e.g., in JT-60U N-NBI system [20].

Modeling of QST 10A source [25-30]
The origin of spatial non-uniformity in multi-cusp arc-discharge sources has been extensively studied by using the QST (National Institute for Quantum and Radiological Science and Technology) 10A source (see figure 7). Not only the experimental studies [21,74], but also the modeling studies [25][26][27][28][29][30] and their comparison of the experimental results have been done so far.
2.3.1.1. Summary of the experimental set-up of QST 10A source Before discussing the modeling results, we briefly summarize the 10A source and experimental set-up. The source is a box type source as shown in figure 7(a). The dimension of the source chamber is 240 mm in the transverse direction (X), 480 mm in the longitudinal direction (Y), and 200 mm in depth (Z). The cusp magnets produce a magnetic field near the anode walls to confine the plasma. The filter magnets are installed near the PG to prevent incoming of fast electrons which have a large probability to destroy negative ions produced on the PG surface. The discharge conditions have been fixed as following values; arc power: 10 kW, arc voltage: 60 V, arc current: 166.7 A, and H 2 gas pressure: 0.30 Pa. Langmuir probe measurements have been done along the longitudinal (Y) direction both in the driver region and extraction region.

Modeling studies of the 10A source and the origin of the non-uniformity
In [25], it has been first shown by neutral transport modeling that the local enhancement of the H atom production by H 2 dissociation due to high electron temperature T e is important to understand the nonuniformity. The local enhancement of H 2 dissociation in the high T e region, and the resultant non-uniformity of the H flux onto the PG is a possible origin of the H − beam non-uniformity observed in the 10A source.
It has been also pointed out in [25] that the spatial profile of surface H − production is mainly determined by the profile of the first flight atomic H flux (FFAF) onto the PG, namely the atomic flux which has reached the PG without undergoing energy relaxation due to collision processes or reflection on the walls. This feature is possibly explained by the following reasons: (1) H atoms produced by H 2 dissociation, i.e., the Franck-Condon neutral atoms have a relatively large energy (>2 eV), (2) the probability of the surface H − production becomes larger for H atoms with higher velocity from the Rasser and Seidel formula [75,76]. The profile of the FFAF on the PG surface is strongly influenced by the spatial distribution of the H atom (H 0 ) production rate in the source chamber.
In [25] above, the H atom birth profile due to H 2 dissociation has been calculated with a given T e profile obtained by Langmuir probe measurements. The high T e component has been taken into account by two temperature fit of the probe I-V characteristics curve. High T e component (>20 eV) has been shown to play a key role for promoting the H atom production and the resulting H − SP, because the cross-section of H 2 dissociation becomes large for the high energy range (E>20 eV) of electrons (see figure 1).
Subsequently, more systematic studies of the EEDF in the 10A source has been done by the KEIO-MARC code in [27][28][29][30]. The calculated electron temperature and density distribution from the EEDF well agree with those obtained from Langmuir probe measurements [28].
It has also been shown that the calculated EEDF in the upper part of the source chamber in figure 7(b) (positions at A and B) consists of a two energy component, as shown in figure 8. Namely, the EEDF consists of (i) The localization of the high energy electrons in the upper part is mainly due to the effects of the magnetic drift in the magnetic field geometry produced by the multi-cusp and filter filed magnet in the QST 10A source. The magnetic drift velocity becomes larger with increasing electron energy. [ / ]. Therefore, fast primary electrons, after being accelerated by the cathodesheath layer, tend to have a large magnetic drift velocity and tend to be localized in the upper region.

Direct comparison of the modeling results with the spectroscopic measurement
Recently, more direct comparisons of the modeling results with the spectroscopic measurements have been done. Based on the EEDF calculated by the KEIO-MARC code, we can calculate the H atom production rate, namely, the number of the H atoms produced by the H 2 -dissociation per unit volume and per unit time. The H 2 -dissociation reactions, taken into account in the analysis, have been already shown in table 1 (Reaction No. 1-No.5) [29]. The rate coefficients are calculated by the EEDF, which includes the high energy tail, as in the same manner as those by equation (4). The H atom production rates r S diss ( ) in each spatial cell in 3D space have been calculated with the rate coefficients. The 2D spatial profile of S diss in the X-Y plane is shown in figure 9(a). As seen from figure 9(a), S diss is localized in the upper part of the source chamber where the high energy electrons are localized.
Starting from these spatial profiles of r S diss ( ) in the source chamber, we can calculate the H atom density profiles from the neutral transport calculation with the following Boltzmann equation, is the velocity distribution function of the H atom in the phase space (r, v H , t). In addition, and m H are the velocity, energy and mass of the hydrogen atom, respectively. The produced atoms are assumed to be Franck-Condon atoms which have the initial energy E=2.15 eV and have an isotropic velocity distribution. In equation (6), S ion and S cx are the ionization loss of the high energy Franck-Condon neutrals and the net source rate due to the charge exchange reactions. The detailed explanation of S ion and S cx are left in [30] because of the page limitation, but the main assumption to evaluate these terms is that the ion distribution function is assumed to be spatially uniform Maxwellian distribution function with ion temperature 1000 K and density 10 18 m −3 . As mentioned in [67], the ion temperature is considered to be in the range between gas temperature (1000 K) and the Franck-Condon energy of atoms.
Using the above assumptions, we solve the Boltzmann equation (6) by directly following H-atom test particle trajectories. For the collision term on the right-hand side of equation (6), we applied the Monte-Carlo NCM [62]. In addition, the particle and energy reflection coefficients for test particles which reach the wall are included in the analysis in the same manner as those in [25,77,78]. Moreover, test particles which reach the extraction aperture are removed from the calculation.
As mentioned above, the production rate of the test particles and the initial position are determined from the production rate S PROD in the real system. Each test particle has a weight ω Η =1.0×10 10 . The source chamber of the QST 10A source is divided into a total 2880 numerical cells with volume ΔV=20 mm×20 mm× 20 mm. After reaching a quasi-steady state, the total number of test particles in each cell, N T , is counted to calculate the H atom density as n H =N T ×ω Η /ΔV. The size of the cell is smaller compared with the mean free paths of ionization and charge exchange of H atoms (both are up to 0.1-1 m) which have the largest reaction rates. The calculation result of H atom density distribution from equation (6) is shown in figure 9(b).
From the H atom density distribution obtained from the transport analysis, the H α line intensity is calculated by the collisional radiative (CR) model [79,80]. As shown in [79,80], the H α line emission rate ε Hα (m −3 s −1 ) in the driver region of negative ion sources is mainly described as, n n X n n X , 7 respectively. The coefficient A H (3, 2) represents the spontaneous transition probability of atomic hydrogen from the excited level 3 to 2. The population coefficients R s (p) (s denotes the particle species; s=H, H 2 ) are the ratio of excited atoms with excited level p=3 to the total density of all excited states. Therefore, the coefficients are denoted as R p n p n q , where n s (p) and n s (q) are the population density in the excited level p and q, respectively. The population densities and the population coefficients are calculated in the CR model in the same manner as in [79,80] from the rate equations, Figure 9. (a) 2D spatial profile of the H atom (H 0 ) production rate by H 2 -dissociation in the X-Y plane in the driver region (at Z=90 mm), (b) 2D spatial profile of atomic density in the driver region (at Z=90 mm) (cited from [30]).
where C s (p, q), F s (p, q) and S s (p, q) are the rate coefficients for excitation, de-excitation and ionization transition from excited level p to q, respectively. The main reactions taken into account in this analysis are in table 3. The rate coefficients are calculated by substituting the non-Maxwellian EEDF from the electron transport analysis by the KEIO-MARC code in the previous subsection. The loss rate of n H2 (p) due to the dissociation is denoted as D H H2 (p) and, conversely, the production rate of n H (q) due to dissociation is denoted as P H H2 (q). From the above transport model (Boltzmann equation) and the CR model (rate equation), the spatial distribution of population density for excited atoms with principal quantum number p=3 [excited H atom: H(3)], which contributes to the H α line emission, is numerically obtained.
It has been reported in [52] that the H α line intensity also shows a strong correlation with the spatial distribution of the negative ion density in the extraction region (distance from the PG to be Z 30  mm) via mutual neutralization (MN) process in this region [83]. However, in the present analysis, we focus on the driver region (Z>80 mm) for the reason that the localization of the high energy electrons in the driver region and the resultant localization of H 2 dissociation are the origin of the non-uniformity in the case of Cs seeded surface H − production. In the driver region, the negative ion density is generally very low because of high electron temperature as discussed above in section 2.2. Therefore, the production of H α line emission in the driver region seems to be mainly due to excitation, de-excitation, ionization and spontaneous radiation of atoms and dissociative excitation of molecules via electron impact (see table 3). Figure 10(a) shows 2D spatial profiles of the H α emission rate in the X-Y plane of the driver region (Z=80 mm). Figure 10(b) shows the H α line intensity obtained in the spectrometry and in the calculation. The spatial (longitudinal) profile of the H 0 density averaged along the line of sight is also shown in the figure.
The spatial profile H α line intensity (blue circles) obtained in the calculation shows a good agreement with that obtained in the spectrometry (red triangles).
The spatial profile of the H α line intensity calculated from the H atom density profile is in the same order as that of the of the H α line intensity obtained from the spectrometry (red triangles). Also, the agreement is good in Table 3. Main reactions taken into account in the collisional radiative (CR) models. The excited levels for atoms and molecules are denoted by p and q (cited from [30]).

No.
Reaction References the quantitative sense, because the absolute values of the line intensity obtained in the calculation and in the spectrometry are both in the order of 10 13 photons s −1 . It should be noted that spatial profile of the H 0 density is different from that of the H α line intensity in the upper region. Namely, the non-uniformity of the H 0 density is weaker than that of the H α line intensity in the negative ion extraction region (Y=−170 to 170 mm). This difference is due to the enhancement of the population ratio for p=3 level by the non-Maxwellian component of EEDF as was discussed in [30].
From the direct comparisons above of the modeling results including the non-equilibrium characteristics of the EEDF with the spectroscopic measurements, now the understanding of the multi-cusp arc-discharge becomes more robust. Namely, it has been confirmed that the local enhancement of H 2 dissociation in the driver region mainly by the high energy component of the EEDF (E>20 eV), and the resultant non-uniformity of the H atom flux onto the PG is a possible origin of the H − beam non-uniformity observed in the QST 10A source.

Modeling of QST JT-60SA kamaboko source
Being based on above the basic physics understanding of the H − beam non-uniformity in the QST 10A source, a super-giant H − source by multi-cusp arc-discharge has been developed for the N-NBI system of plasma heating and non-inductive current drive for the JT-60SA tokamak [24], which is currently one of the largest tokamaks with superconducting magnets in the world and now under construction as an EU-JAPAN international collaboration project.
The final goal of the JT-60SA N-NBI system is to produce 22 A negative ion beams for the long pulse operation of up to 100 s [22,23,84]. The KEIO-MARC code has been also applied to the systematic simulation study of the EEDF in JT-60SA H − Source. It plays a significant role in understanding the plasma non-uniformity in JT-60SA H − Source and to succeed in achieving such a large H − current (22 A) extracted from the source. As shown in figure 11(a), the arc-discharge chamber consists of the main chamber wall with the semicylindrical shape, two side walls and the PG on the bottom. Due to its semi-cylindrical shape, it is sometimes called the Kamaboko source originally from the shape of the Japanese food 'Kamaboko'. The arc chamber size is significantly large with dimensions: 340 mm in radius, 1200 mm in length. The PG has a large extraction area (450 mm 1100 mḿ ), which is further divided into 5 segments. Each segment has a large number of extraction apertures (240) and totally 1080 apertures exist for the whole extraction area. Totally 48 filaments are inserted from the side part of the main wall. In order to confine the plasma and to obtain a high plasma density and temperature in the driver region, a multi-cusp magnetic configuration has been produced by 26 rods of the permanent magnet installed along the longitudinal direction (Y-direction) on the main chamber wall, and also the 10 rods on the side walls along the transverse direction (X-direction) in figure 11(a). In addition to these permanent magnets for the multi-cusp magnetic field in the driver region, it is necessary to produce the magnetic filter configuration in the extraction region close to the PG to reduce the electron temperature and to suppress co-extraction of electrons from the extraction aperture. For this purpose, the PG magnetic filter configuration (PG-filter) has been adopted in the JT-60SA H − source [22,23]. Namely, the electric current I PG is forced to flow inside the PG in the longitudinal direction in order to produce and strengthen the transverse filter magnetic field in the X-direction in figure 11(a).

Modification of the magnetic field configuration to improve the uniformity
In the original PG-filter magnetic configuration shown in figures 11(a) and (b) a strong non-uniformity of the extracted current has been observed as will be shown later (see figure 17). This non-uniformity made it impossible to fulfill the requirement of the extracted beam current (I 22 A H ext > -) for JT-60SA. In the PG-filter configuration above, magnetic drifts are expected to take place along the chamber wall in the same manner as QST 10A source discussed in section 2.3.1. The direction of the electron magnetic drift ( B B -´ ) in the bottom of the driver region is schematically drawn in figure 11(b). In this region, the PG magnetic filter field is relatively small, but still exists. The PG-filter magnetic field B directs towards the positive X-direction, while its gradient directs B  towards the PG in the negative Z-direction. As a result, the electron magnetic drifts (− B B  ) are expected to be in the longitudinal (−Y) direction in figure 11(a).
From the above basic consideration and also the numerical simulation by the KEIO-MARC code, modification/improvement of the magnetic configuration has been done by QST in order to improve the uniformity observed in the experiments and to fulfill the requirement of the extracted current. These modifications have been done under the close collaboration between the experimental study at QST and the modeling study at Keio university.
The conventional PG-filter configuration had been modified as a tent magnetic filter (tent-filter) configuration shown in figure 12(a). In this tent-filter configuration, a pair of the permanent magnets is installed in parallel on the top of the main chamber in the longitudinal direction. These two magnets put their surface with the N-pole down towards the PG as shown in figure 12(b). In the same way, the other pair of magnets are installed on the bottom part of the arc chamber, but they put their surface with the S-pole up towards the top of the main chamber. As a result, the magnetic field line starting from the top magnets ends at the bottom magnets as shown in figure 12(b).
The shape formed by magnetic field lines from the top magnet to the bottom magnet looks like a tent shape pitched on the PG as a ground. By these modifications, as has been shown in [22], it is experimentally confirmed that (1) beam uniformity is significantly improved (from 68% to 83%) over an area of the whole extraction area of 450×1100 mm 2 , and (2) the improvement of the beam uniformity leads to the production of 32 A H − ion beams with the whole extraction area. Although the improvement for the long pulse (100 s) and high energy H − beam acceleration (500 keV) will be still needed, the obtained value of the extracted current 32 A is significantly larger than the target value (22 A).

Modeling study of JT-60SA H − source
In this section, we discuss the role/contribution of the modeling to the above experimental achievement of JT-60SA H − source.
The KEIO-MARC code has been applied step by step to each magnetic configuration shown in figures 11 and 12 above by taking into account the real 3D geometry and dimensions of the source. The following discharge parameters are used for a typical discharge condition in the experiments; arc power P arc =120 kW, arc voltage V arc =100 V, arc current I arc =1.2 kA and the initial gas pressure P gas =0.3 Pa. It should be noted the arc power is relatively large in comparison with those used in the simulation studies of QST 10A H − source discussed in section 2.3.1, but the input power density per unit volume is not so different, because the source volume of the JT-60SA H − source is larger.
The collision species taken into account for this study are almost same as those in section 2.1 (see table 1). Also, the background plasma (H + , H , 2 + H 3 + ) and neutral (H 2 , H) density and temperature have been fixed and the same values as in section 2.2 have been used. The simulation domain is divided into 456280 cells, each with the size Δx=Δy=Δz=10 mm. The time is set to be Δt=1.0×10 −10 s to integrate the equation of motion, which is much smaller than the cyclotron frequency. The collision time for the Coulomb collision is set to be Δt=1.0×10 −8 s, which is smaller than the slowing down time of the Coulomb collision. Figure 13 compares 2D profiles of the field strength (color plot) and the structure of the magnetic field lines in the X-Z cross-section with Y=0. Figures (a) and (b) correspond to the (a) conventional PG magnetic filter configuration (PG-filter: see figure 11) and (b) tent magnetic filter magnetic configuration (tent-filter: see figure 12), respectively. Figure 14 shows simulation results of the 2D electron density profile in the X-Z plane at Y=0 mm for (a) conventional PG-filter configuration and (b) tent-filter configuration, respectively. The deference in the density profiles shown in figures 14(a) and (b) is clear from the difference in the magnetic configuration shown in figures 13(a) and (b). As seen from the electron density profile in figure 14(b), the tent-shaped low-density region has been clearly formed. Figures 15(a) and (b) compare the average direction of velocity vector for test particles in the X-Y plane at Z=160 mm. The average velocity vector has been obtained by averaging over the velocities of test particles in each spatial cell and the direction of each average velocity is shown at each position in figure 15.
In the (a) PG-filter configuration, it is confirmed that the electrons drift towards the lower part (−Ydirection) along the longitudinal direction in the central part of the X-Y plane in the deriver region as shown in figure 15(a). Due to this drift in the X-Y plane, electrons are localized in the lower part of the chamber. In the same manner as in the 10A H − source discussed in section 2.3.1, this localization of the electrons, especially electrons with relatively high energy, leads to the enhancement of the H 2 dissociation in the driver region.
On the other hand, (b) in the tent-filter configuration, the rotation of the electrons has been promoted in the central part of the X-Y plane in the deriver region as seen from figure 15(b), because in this magnetic configuration the magnetic drifts along the longitudinal direction are in the opposite direction as have been already shown in figure 12(b).
In the same way as in the modeling of the 10A source, the spatial profiles of the H atom production rate by H 2 -dissociation for (a) the PG-filter and (b) the tent-filter cases have been calculated from the EEDFs obtained by the KEIO-MARC code. Subsequently, H atom density has been calculated from the Boltzmann equation for H atom as in section 2.3.1 (see equation (6)). The spatial profiles of the H atom flux (first flight atomic flux: FFAF), which is important for the surface H − production as has been pointed out in [25] and discussed in section 2.3.1, have been also calculated. The calculated H atom (FFAF) fluxes are shown in figure 16(a) along the longitudinal direction. In figure 16(a), the results for the two cases, i.e. for the conventional PG-filter (blue broken line) case and the new tent filter (red solid line) case are compared. As seen from figure 16(a), the spatial non-uniformity of the H atom FFAF on the PG is largely improved for the tent-filter configuration in comparison with the PG-filter configuration.
Using these H atom fluxes onto the PG in figure 16(a), we have also calculated the 2D spatial of the surface H − production amount on the PG in the following manner, Here, N , Hin f E H in ( ) is the number of H atoms reaching the PG, and the energy distribution function of the incoming atoms. In addition, the symbol E b ( ) in equation (11) is the H atoms' transfer rate to negative ions calculated by the equation from Rasser [75] and Seidl [76]. The Hproduction rate calculated from equation (11) Figure 13. 2D magnetic field profiles in the X-Z cross-section with Y=0: (a) conventional PG magnetic filter configuration (see figure 11) and (b) tent magnetic filter configuration (see figure 12). The field strength is shown in the color plot and the magnetic field lines are shown with white lines. for the PG-filter and the tent-filter are also compared in figure 16(b). As seen from the calculation results, the spatial uniformity of the H − SP is significantly improved for the tent-filter configuration.
The experimental results of spatial profiles of the extracted H − beam current are also shown in figure 17 [22]. The open circles are the results for the PG-filter configuration, while the closed circles are those for the tentfilter configuration. In the experiments, it is confirmed that the spatial uniformity of the extracted H − beam is largely improved by using the tent-filter configuration as in the simulation.
As mentioned at the beginning of this section, more quantitative evaluation of the improvement has been done in [22]. The beam uniformity is significantly improved (from 68% to 83%) over an area of the whole extraction area of 450×1100 mm 2 of JT-60SA H − ion source, where the uniformity [%] is defined by the fraction of the measured beam current integrated over the whole extraction area to the total beam current estimated by assuming the peak beam intensity to be constant on the whole extraction area.
Finally, it should be noted the KEIO-MARC code uses various simplifications in order to perform full 3D simulations with real 3D size and real 3D magnetic configuration of multi-cusp arc-discharge H − sources. Especially, assumptions of the electric field and the sheath potential drop implemented in the KAIO-MARC  code are rather crude. Therefore, it is necessary to validate carefully the code by comparisons with experiments. As shown in sections 2.1-2.3 above, the modeling results have been compared with the experimental results, not only in one discharge/source condition, but also in various discharge/source conditions up to now.
As seen from the results in sections 2.1-2.3, the modeling results at least give a useful physics insight to help us understand characteristic features of multi-cusp arc-discharge plasmas and to help improve the source performances. For example, the modeling results have largely contributed to improving the spatial uniformity of the H − production and the resultant increase in the H − beam current in giant H − sources for the magnetic fusion application as was discussed in section 2.3.
Taking into account the effect of the electric field and the sheath potential drop more precisely, it is necessary to perform full kinetic simulations (e.g. electrostatic PIC model with Poisson's eq). Due to the large plasma volume of such giant H − ion sources, however, it is almost impossible to perform full-kinetic 3D simulations with realistic size of H − sources, because we need a huge amount of computational cost. As will be shown later in section 4, even for the very small plasma volume close to the PG, it is necessary for most of the present 3D models to use a reduced-size model due to the massive computational cost.
On the other hand, it is easy for fluid models to perform such large-scale simulations and to take into account the effect of the electric field by current continuity equation. However, the fluid model cannot take into account the non-equilibrium features of the EEDF, especially the effect of the high energy tail of the EEDF. As is shown in sections 2.1-2.3 by the test particle Monte-Carlo modeling based on the KEIO-MARC code, the high energy tail of the EEDF (e.g. see figures 4 and 8) plays a key role for both the volume of H − production and the surface H − production.

Modeling of RF-ICP sources
As mentioned in section 1, here, we focus on a series of the modeling studies in [33][34][35] of the RF-ICP H − source for the CERN Linac4 accelerator [8]. In RF-ICP plasmas, understanding the E-H transition is one of the key issues especially for the impedance matching between the RF-system and plasma load, because the plasma characteristics and the resultant plasma loading impedance significantly change before and after the transition [10,13,85].
However, full kinetic-modeling studies of the E-H transition have been very limited up to now. Here, we briefly review our full kinetic modeling studies of the E-H transition by 2.5D EM-PIC model. From the EEDF obtained from the EM-PIC modeling, the H α intensity in the Linac4 H − source has been also calculated by the transport model and the CR model in the same manner as in section 2.3.2. The direct comparison of the modeling results with those of the experimental results is discussed. [33][34][35] First, we briefly summarize the characteristic feature of the Linac4 H − source [8,34]. Figure 18   coil is installed externally around the plasma chamber. It is operated at 2 MHz with an available peak power of 100 kW (maximum 300 A coil current). The plasma is ignited for 500 μs with a repetition rate up to 2 Hz. Three optical view ports are mounted on the back side to capture the light emitted by the plasma, pointing to: (i) the plasma chamber axis, (ii) the center of the solenoid, and (iii) the plasma electrode.

Linac4 RF-ICP H − source and EM-PIC model
The EM-PIC model and main simulation parameters used in [34] are summarized in figure 19 and table 4. Detailed descriptions of the model and calculation conditions have been given in [34]. Here, we briefly summarize basic equations, main assumptions and numerical techniques.
The EM-PIC model is a 2.5D code in the cylindrical coordinate system, in which the electromagnetic fields generated by the plasma (E pl , B pl ) are solved in the 2D space by assuming axisymmetry ( 0 q ¶ ¶ = / ), while the particle motions are followed in the 3D space. Specifically, E pl , B pl are solved with the finite difference time domain method [14,15] with absorbing boundary conditions: Mur [86] at the axial boundary, while Bayliss and Turkel [87] at the radial boundary. The dimensions of the simulation domain are twice the plasma chamber size to avoid large reflections at the boundary. The RF electromagnetic fields in vacuum (E RF , B RF ) are simulated in the frequency domain by the software package (Ansys HFSS), starting from the ion source 3D model with its respective material properties. Details of the simulations are described in [88]. From Ansys HFSS we export a field map of E RF , B RF that we average in the θ direction.
In the simulation, the trajectories of the following charged particle species, electron, H + , H 2 + are computed by solving their equations of motion. The Boris-Benemann version of the Leap-Frog method [15] has been used and weighting to the grid is performed using a volume (r 2 − z) bilinear interpolation function with volumes treated following Verboncoeur [89]. The plasma current j required for the calculation of Epl, Bpl, is the average Figure 18. Longitudinal and transversal cross-section of the Linac4 H − ion source plasma generator. The plasma chamber is made of Al 2 O 3 . The RF-coil is surrounded by 6 ferrites and embedded in epoxy to avoid RF breakdown. Gas is injected in pulsed mode by a piezo-valve. The optical view port has a capture angle of 3°. The magnetic Halbach octupole is formed by alternating magnets with clockwise and counter-clockwise magnetization. The filter field is a dipole magnet with vertical magnetization. The Cesium inlet is used to evaporate Cs (cited from [34]). Figure 19. Computational cycle performed at each time step. The plasma current j is the sum of electrons, protons and H 2 + contribution (index s). Weighting from the particle to the grid (and vice versa) as required by the PIC algorithm is implied. The transition between the 3D to the 2D is performed by averaging in the θ direction (cited from [34]). value in the θ direction. Cell size and time step are chosen to satisfy the stability conditions of the EM-PIC scheme [14,15].
Electron-neutral collisions are taken into account by a null-collision Monte-Carlo method [62] in almost the same manner as sections 2.2 and 2.3. Coulomb collisions, however, are not taken into account in this analysis, because the e-neutral collision frequency is much larger than the Coulomb collision one for the H 2 pressures (1-10 Pa) and plasma densities (<10 16 m −3 ) simulated [90]. Figure 20 shows spatial profiles of electron density n e and the axial component of the electric field at a phase of the maximum RF E-field inside the chamber for a typical low plasma density E-mode phase. Figure 21 shows the same plots, but at the opposite phase of the maximum RF E-field. As seen from the comparison of n e between figures 20 and 21 together with that of the axial electric field, the axial motion of the electrons is one of the important characteristic features for this E-mode phase. Depending on the E z -field polarity, electrons are pushed back and forth from the coil region towards the end walls at both sides in the axial z-direction.

Modeling results of plasma characteristics in E-mode and H-mode
On the other hand, in a relatively high n e regime after the E-H transition, the azimuthal electric field becomes dominant as shown in figure 22. The high-density region is localized in the region where the 5-turn RFcoil exists.  Figure 20. Distribution of the plasma parameters in the plasma chamber. From top to bottom: electron density (n e ), RF axial E-field (E z,RF ), plasma axial E-field (E z , pl ), total axial E-field (E z ); which are typical for a low density plasma driven by the capacitive field at a phase of maximum E-field. The plasma responds by generating E z,pl opposite to E z,RF resulting in a total E z that is reduced in the plasma region (cited from [34]). [34,35] In order to validate the plasma characteristics in section 3.2 obtained by electromagnetic particle in cell with Monte-Carlo method for collisions (EM-PIC-MCC) plasma modeling, comparisons have been made between the photometry and spectroscopic measurements [34,35]. Photons incoming to the plasma viewing port (see  , total axial E-field (E z ), RF azimuthal E-field (E θ , RF ), plasma azimuthal E-field (E θ , pl ), total azimuthal E-field (E θ ), which are representatives of the density ramp-up in H-mode phase. While the axial field is shielded by the plasma, the azimuthal component fully penetrates into the plasma region and drives the discharge (cited from [34]).

Model validation: comparisons with experiments
figure 18 left) from the plasma, is focused by quartz lens in the viewing port and transmitted to 3 photomultiplier tubes (PMT) by optical fiber. At each PMT, optical filters are set at the entrance to cut off the visible lights except Balmer H α (λ=656.28 nm), H β (λ=486.13 nm) and Hγ (λ=434.05 nm) lines from hydrogen atoms. The line integrated H α line intensity has been measured. Further details of the experimental set up are given in [34,35]. In order to compare the modeling results with the photometry and spectroscopic measurements, the H atom transport and collisional radiative (CR) model [30] developed in section 2.2 have been coupled with the results by the EM-PIC-MCC plasma model. In section 2.2, the EEDF calculated by the KEIO-MARC code has been used to calculate the rate coefficients of the H atom excitation and the resultant H α -line emission rate. Here, the EEDF calculated by the EM-PIC-MCC model in each spatial cell has been used to evaluate the spatial profile of the H α -line emission in the same manner as in section 2.2. Figure 23 shows the experimental result of the time evolution of the H α line intensity during the low-density plasma ignition phase, namely in the E-mode phase corresponding to figures 20 and 21. In figure 23, the time evolution of RF-coil current has been also shown.
It is interesting to note that two peaks with the different intensity, namely one large peak and one small peak, have been observed in each one RF-cycle. Figure 24 shows the modeling results. The same tendency, i.e., asymmetric peak strength in one RF-period can be clearly seen.
As is discussed in [34,35], the higher peak results from the E-field pushing the electrons towards the view port (see figure 18 left: the optical view port is located at the left end of the plasma chamber), while the smaller peak from the plasma centered in the plasma chamber where the RF-coil is located (see also figure 18 left: RF-coil located a little more right from the center of the plasma chamber).
As shown in figures 20 and 21, the main capacitive component of the electric field in the low-density E-mode phase is the axial field E z . Depending on the RF-phase, electrons are accelerated and decelerated in the zdirection. In other words, the axial electric field pushes electrons forward to the view port, and backward  towards the coil region depending on the polarity of the RF E z -field, which in turn leads to the localization of the electron density in the z-direction depending on the RF-phase.
If the electron density close to the viewing point is larger, then the H α -intensity becomes larger, because the light emission roughly scales with electron density. On the other hand, if the electron density close to the viewing point becomes smaller, then the H α -intensity also becomes smaller.
These processes lead to the asymmetric time evolution of H α -intensity within one RF-cycle observed in the experiments as shown in figure 23.
Not shown here, but both the measured and calculated results of the H α intensity become symmetric in one RF-cycle as the discharge becomes dominated by the inductive component. This result is possibly explained by the modeling result shown in figure 22. In this H-mode phase, the plasma remains centered in the coil region as E z,RF is shielded by the plasma. As electrons localized in the center region are periodically accelerated and decelerated by E θ , RF , the emission increases and decreases respectively giving rise to symmetric property.
It is generally pointed out that the E-H mode transition depends on the impedance matching network system (e.g., in [84]). In fact, in order to optimize the power transfer and impedance matching between the RFsystem and the RF-ICP in Linac4 H − source [8], an equivalent circuit model has been developed in [91,92]. This model describes not only the RF-ICP inside the ion source, but also includes the matching circuit system. All the elements have been modeled as lumped elements, which means this model is based on the 0D model. It has been confirmed that the control of the RF signal frequency is useful in increasing the efficiency of the power transfer to the RF-ICP, which has empirically been found in the Linac4 H − source experiments [93,94].
Most of the 'modeling' studies of the E-H mode transition so far have been based on the fluid or equivalent circuit model, while only a few 'full-kinetic' modeling studies of the E-H mode transition can be found up to now except for the results explained above in section 3. Recently, based on almost the same EM-PIC model in section 3.1, kinetic aspects (e.g. the change in the energy distribution function (EDF) of the electrons and ions before and after the E-H transition) have been studied more in detail in [95,96]. The EDFs deviate from the Maxwell distribution. This indicates the importance of the kinetic perspective to investigate the RF-ICP discharge.
Although recently the numerical code has been improved and speeded up with full implicit numerical scheme [97], it is still too massive for such full-kinetic 2.5D EM-PIC-MCC modeling to take into account directly the effect of the matching circuits or to couple directly with the model of the matching network system mainly from the viewpoint of the computational cost at least at present.
However, it is expected that the developed full-kinetic 2.5D EM-PIC model enables us to calculate the impedance of RF-ICP with the non-equilibrium EDFs and to implement the resultant plasma impedance into the 0D equivalent circuit model above including the effects of the network circuit system. This model will be the basis for further and more reliable modeling studies of the impedance matching between RF-system and RF-ICP. The model includes not only the matching network system but also the kinetic effects of the RF-ICP.

Extraction physics
As mentioned in section 1, an understanding of the H − extraction process from the H − sources is one of most common and critical issues of H − source development.
In various Cs-seeded negative ion sources, an 'ion-ion plasma' layer which mainly consists of positive and negative ions has been observed in front of the PG [51,52,98,99]. Especially, recent experiments on the ion-ion plasma in [51] triggered a series of 2D ES-PIC modeling studies [39][40][41][42] in order to find the answer for the questions: (i) 'why such a thick ion-ion plasma layer is produced under the H − SP?', and 'from where H − ions are extracted? Directly from the surface of the PG/extraction hole? Or from the bulk of the ion-ion plasma?
2D ES-PIC studies above pointed out that (i) electron loss along the magnetic field line in front of the PG is a crucial role in the production of the ion-ion plasma which almost consists of H + and H − without electrons [39][40][41][42], (ii) meniscus (H − emitting surface) in such an ion-ion plasma which has a relatively large curvature [41,42], especially in the peripheral region of the meniscus, (iii) H − ions extracted through the peripheral part of the meniscus with large curvature become a beam halo component [41,42], and (iv) most of the beam halo components are directly from the surface of the PG/extraction aperture, while the beam bulk component mainly consists of those H − ions through the bulk plasma [41,42].
Recently, a study of the ion-ion plasma has been done with the 3D ES-PIC model [56,57] to confirm the results obtained by the 2D ES-PIC modeling studies and to give more quantitative discussion. Following [56], we first discuss the effects of the electron loss along the magnetic field line on the ion-ion plasma formation in section 4.1. Next, in section 4.2, a comparison between 2D and 3D PIC modeling will be done. Especially, the plasma meniscus shape and the resultant beam optics will be compared between the 2D and 3D model under almost the same calculation conditions. In section 4.3, integrated 3D modeling [57] from the H − extraction to the beam acceleration and a more quantitative comparison with the experiments will be presented.   (12) and (13) are solved in 3D model geometry (see figure 25) of the extraction region in the negative ion source. The Boris-Benemann version of the leap-frog method [15] is used to solve equation (12), while the bi-conjugate gradient stabilized method (BiCGSTAB) [100] is used for equation (13). The simulation domain is shown in real-size with 17 mm 19 mm 19 mm.´The x-axis is taken to be the direction of the H − with 17 mm 19 mm 19 mm.´The x-axis is taken to be the direction of the H − extraction, while the z-axis is parallel to the direction of the PG magnetic filter field. The strength of the PG magnetic filter field is specified from the experimental data in [22]. The magnetic field by the electron suppression magnet (ESM) is not taken into account for simplicity.
As mentioned in section 1, the reduced size scaling used in the previous model [36,[38][39][40][41][42] has also been adopted in this 3D modeling. The system length in the simulation (L sim ) is given as L sim =sL real , with s being the scaling factor 3.7 10 2 =´-( ) and L real being the characteristic scale length in the real system. A more detailed description of the reduced size scaling can be found in [38][39][40][41][42]. The numerical grid has 137 153 153´in the x, y, and z directions with the following dimensions of the mesh:  is equivalent to 20-30 particles in each numerical cell. The initial ratio of the number of particles for each species is assumed to be N H+ :N e :N H-=10:9:1.
The surface produced negative ions are launched at the PG surface by a given number of particles per each time step. The initial velocity of the surface produced negative ions has been sampled from the half-Maxwellian with the temperature T b shown in table 5. Here, the number of surface produced negative ions per each time step is determined by the emission rate (mA cm −2 ) of the surface produced negative ions, j b . Here, j b is set to be ∼80 mA cm −2 which is almost the same values as used in [39][40][41][42].
The diffusion of electrons across the PG magnetic filter field (see figure 26) due to Coulomb collisions with positive ions (H + ) and elastic collisions with hydrogen molecules have been taken into account by the simple random-walk model [38].

Modeling of electron loss along the magnetic field line
As was pointed out by the series of 2D ES-PIC modeling [39][40][41][42], among various physical parameters, electron loss time along the magnetic field line is one of the key parameters to understanding ion-ion plasma formation. In most of PIC modeling studies of the multi-aperture system, only the single aperture has been modeled, even in full 3D ES-PIC simulations [53][54][55]. Careful modeling of the non-local electron loss from the whole system is very important. In other words, not only the magnetic configuration close to the aperture, but also large scale and non-local magnetic configuration/topology in the entire source is important, because electrons are easily lost along the field line.
Here, we take into account the electron loss along the field line in the same manner as in [41,42], i.e., 't t/ // -model' is used. In the t t/ // -model, the step length x D D of the random-walk process is given by,    (14) is as follows. Electrons lost along the magnetic field line cannot contribute to the diffusion across the magnetic field. A more detailed explanation of the model and derivation of equation (14) is given in [41,101].
Based on [101], t t/ // is roughly estimated to be in the order of 10 −2 in a large scale negative ion source, e.g., the negative ion source for JT-60U. However, as the strength of the magnetic filter field increases towards the chamber wall due to the external magnet, some of the electrons will be reflected from the wall side to the inner source plasma region due to the magnetic mirror trap. In order to take into account this mirror effect, we assume that 0.08 t t = / // below. The numerical results will be compared with the effect ( 0.08 t t = / // ) and without the effect ( 1.0 t t = / // ) of the electron loss along the field line. Figure 27 shows the time-averaged 2D profiles of electronegativity ( n n e H a = -/ ) in the x-z mid-plane (see figure 25) which is parallel to the PG-filter field. Figures (a) and (b) correspond to the case (a) 0.08

Effect of the electron loss along the magnetic field line in the ion-ion plasma formation
respectively. In the experiment in [51], the V-I prove characteristic curve in front of the PG becomes quite symmetric under the ion-ion plasma case, i.e., the obtained positive ion saturation current I + is equal to the negative ion saturation current I − . Therefore, from the ratio of the I − to the electron saturation current I e , here we define the ion-ion plasma region where the following inequality is satisfied; m m 40. respectively. In figure 27, x=9 mm corresponds to the PG surface (see figure 25).
According to the above definition of the ion-ion plasma layer ( 40 a > ), the thickness of the ion-ion plasma is estimated to be 3.2 mm for case (a) and 1.2 mm for case (b), respectively. From these results, it is confirmed that the ion-ion plasma is formed with relatively large thickness for case (a). On the other hand, the ion-ion plasma layer is formed just in front of the PG in case (b). These comparisons clearly show that the thickness of From these 3D results, our previous results using the 2D ES-PIC model of ion-ion plasma formation [39][40][41][42]101] have been confirmed by the 3D ES-PIC model [52]. The electron loss along the magnetic filter field is shown to be very important for ion-ion plasma formation. The thickness of the ion-ion plasma depends on the electron loss along the field line. These results, however, are based on the relatively simple modeling of the electron transport with t t/ // -model. In order to quantitatively validate the effect of the electron loss on the ion-ion plasma formation, the 'non-local' magnetic field structure or topology all over the ion source should be self-consistently taken into account as well as the local magnetic field structure near the PG in the single aperture model. Generally, electrons move easily along the field line and the electron loss time along the field line is much smaller than the diffusion time across the magnetic field. Therefore, the non-local magnetic field structure directly affects the electron confinement/loss time in the extraction region. For the development such a 3D ES-PIC model including global magnetic field structure or topology, however, it is necessary to overcome the problem concerning the huge increase in the computational cost.

2D versus 3D
PIC modeling-comparison of meniscus shape and beam halo formation with a relatively simple model geometry As mentioned in section 1, 2D ES-PIC modeling is quite useful for the basic understanding of plasma dynamics and the potential structure in the extraction region. It is, however, far from the perfect. It is necessary to develop the 3D model for more quantitative discussion and also from the viewpoint of the design oriented study of the extraction and acceleration system.
In this section, comparisons will be made between the 2D and 3D modeling. Focus is placed on the meniscus shape calculated by the 2D and 3D model, and the resulting beam optics. It has been reported that the negative ion beam consists of two Gaussian parts, namely, the beam core with good beam optics and a beam halo with poor beam optics. It is necessary to suppress the beam halo, because the beam halo reduces the negative ion beam transmission and causes the heat loads on the extraction and acceleration grids. Especially, the heat load by the beam halo is one of the critical issues to develop the high-power N-NBI system for fusion application as was discussed in section 2.3.
Here in section 4.2, the basic comparisons of the meniscus shape and the resultant beam halo will be made with a relatively simple geometrical model. In section 4.3, more quantitative discussion will be given with more realistic 3D modeling. 2D and 3D model geometry used in this comparison are shown in figures 29 and 30, respectively. As seen from the comparison between figures 29 and 30, the extraction aperture is modeled as a slit with the infinite length in 2D model, while it is modeled by more realistic cylindrical hole with a finite area in the 3D model.
For these simple 2D and 3D geometry, the basic equations (12) and (13) have been solved under almost the same conditions except for the spatial dimension and the modeling of the extraction aperture. A more detailed model description has already been given in the previous section. In both the 2D and 3D model, the main assumptions are the same. They are summarized as follows: (i) the reduced size scaling with s=3.   Here, we define this contour as a plasma meniscus, because most of the H + ions cannot go forward beyond this contour and they return towards the lefthand side i.e. the source region. On the other hand, H − ions start accelerating towards the EG, because the electric potential becomes positive and is getting larger and larger from this contour towards the right boundary, i.e., towards the extraction grid.
As seen from the comparison between figures 31(a) and (b), the plasma meniscus in the 2D model penetrates more deeply into the plasma source region. As a result, the curvature of the meniscus becomes larger in the peripheral part of the meniscus. In the 3D model, the aperture is modeled as a 'realistic hole', while it is modeled as a 'infinite slit' in the 2D model as mentioned above. This is the main reason why the penetration of the meniscus is larger in the 2D model.
Due to this deep penetration of the extraction field, the curvature of the equipotential contour at the peripheral part of the meniscus close to the PG becomes large compared with the 3D result (see figure 31). Surface  produced H − ions are much affected by the electric field close to the PG and they are accelerated to the direction perpendicular to the equipotential contour. As will be shown below, this large curvature of the equipotential contour at the peripheral part of the meniscus close to the PG leads to the large fraction of over-focused beam component, i.e., beam halo component in the 2D model.
In It should be noted that in figure 32 the H − ion beam is convergent when the sign of the beam divergence angle θ is negative (θ<0) for y>10 mm, or positive (θ>0) for y<10 mm in figure 32. On the other hand, the H − ion beam is divergent when the sign of the beam divergence angle θ is positive (θ>0) for y>10 mm, or negative (θ<0) for y<10 mm.
From figures 32(a) and (b) one can see the H − beam profile consists of two components. The one is the beam component in the central region. This beam component is convergent and corresponds to the beam core. The other is the beam component in the peripheral region. This beam component is divergent and corresponds to the beam halo. From the comparison between figures 32(a) and (b) the beam halo in the 2D model is larger than in the 3D model. This is because the curvature of the plasma meniscus in the 2D model is larger than in the 3D model.
From these comparisons, the beam halo mechanism proposed by the 2D ES-PIC model is confirmed and essentially unchanged even with the 3D model. Namely, the negative ions extracted from the periphery of the meniscus are over-focused in the extractor due to the electrostatic lens effect, and consequently become the beam halo. In the quantitative aspects, however, there exists large differences. The beam halo fraction to the total beam current is calculated to be 60% in the 2D model from figure 32, while the fraction of the beam halo is significantly reduced to 4.5% in the 3D modeling. As will be discussed in the next section, this 3D result reasonably agrees with the experimental result, while more careful modeling of the extraction voltage as in [41,101] is needed for the 2D modeling to obtain more quantitative comparison with the experiments.

Comparison with the experiments with 3D ES-PIC model [57]
As mentioned in sections 1 and 4.2 above, to suppress the beam halo component is one of the crucial issues for the high-power H − beams acceleration because the beam halo causes the heat loads on the extraction and acceleration grids.
In order to compare the simulation results and experimental results of the heat loads due to the beam halo component, 2D PIC integrated modeling from plasma meniscus formation to the beam acceleration has been first done in [42] for the model geometry of the JT-60U N-NBI (negative-neutral beam injection) system (see figure 33). The JT-60U N-NBI system has been designed to deliver the 500 keV, 22 A hydrogen negative ion beams and a rated pulse length of 10 s [20]. To achieve high energy acceleration up to 500 keV, JT-60U N-NBI beam system consists of mainly two parts: (1) multi-grid extractor and (2) multi-grid accelerator as shown in   (1) consists of three grids, i.e., the PG, EXtraction grid (EXG) and electro suppression grid (ESG), while the latter (2) consists of A1G, A2G and GRG grids (A1G: acceleration grid 1, A2G: acceleration grid 2, and GRG: grounded grid, respectively). The 2D integrated model in [42] has reproduced the basic trend of the beam halo experimentally observed in [50], e.g., (i) the beam halo is intercepted on the grids in the accelerator, while the beam core passes through the grounded grid (GRG: see figure 33) aperture without interception, and (ii) the heat loads in the accelerator due to the intercepted negative ions become larger in the downstream grid (from A1G to GRG in figure 33).
For more quantitative comparison, however, the 2D model is not enough. As can be expected from the comparison between the 2D and 3D model of the plasma meniscus shape in the previous section, 2D model overestimates the amount of the beam halo component than the 3D model. As a result, the heat load due to the beam halo on the extractor, e.g., especially the heat load on the extraction grid (EXG) becomes too large in 2D model, because the deep penetration of the plasma meniscus and the resultant large curvature of meniscus lead to a large amount of the over-focusing beam halo component. On the other hand, the heat loads on the first acceleration grid (A1G) and second acceleration grid (A2G) due to the beam halo are much smaller than the experimental result in [50].
In order to make more quantitative comparisons with the experiments, 3D integrated modeling from plasma meniscus formation to the beam acceleration has been recently started in [57]. In the following, we briefly summarize this integrated 3D model and its results.
In the model geometry for the 3D model in figure 33, the cross-sectional view of the x-y plane is shown. As in the previous section, the extraction aperture has been modeled as a cylindrical shape in the present 3D model instead of the infinite slit in the previous 2D model. In this model geometry, the negative ion trajectories from the plasma meniscus to the extractor are calculated by using the 3D PIC code in which the plasma meniscus is calculated self-consistently without any assumptions. The 3D PIC model and numerical method are the same as the model described in section 4.1 except for the model geometry. Namely, the present 3D model is mainly based on the following assumptions: (i) reduced size scaling with the same scaling factor, (ii) t t/ // -model is applied to reproduce the ion-ion plasma formation under the Cs-seeded H − SP and (iii) the surface produced negative ions are modeled to be launched from the surface of PG uniformly. The initial velocity distribution is assumed to be the half-Maxwellian, while the initial temperature is taken to be smaller 0.1 eV in this calculation from the typical emittance of H − ion beams.
The negative ion trajectories in the accelerator are calculated by using the commercial software (Omnitrak, Advanced Science Laboratory, Inc.) with the boundary condition that the potential at the location of 20 mm away from the exit of the extractor is given to be 52.4 kV. In figure 33, an example of the calculation result of contour maps of the electrostatic potentials are also shown. The x-axis is taken to be the direction of negative ion extraction and acceleration. The magnetic filter and the ESM field are taken into account. The y-axis is parallel to the direction of the magnetic filter, while the z-axis is parallel to the direction of the electron suppression magnetic field. As mentioned above, the extractor consists of the extraction grid (EXG) and the electron suppression grid (ESG). The permanent magnets are equipped in the EXG to suppress the electrons. Model geometry of the multi-stage accelerator is taken from the JT-60U N-NBI system in [102,103]. The accelerator consists of the acceleration grid 1 (A1G), the acceleration grid 2 (A2G), and the grounded grid (GRG). Each voltage is taken to be 174.7 kV, 341.3 kV, and 508.5 kV, respectively. Figure 34 shows a typical example of the negative ion trajectories of the beam core and halo in the extractor. The contour map of the electrostatic potential is also shown in figure 34. Moreover, the source region around the PG (4 mmx12 mm) is enlarged. The beam core is extracted around the center of the meniscus where the curvature of the equipotential surface is small. This negative ion is focused by the electric field around the ESG and then is convergent downstream the exit of the ESG. On the other hand, the beam halo is extracted from the periphery of the meniscus where the curvature of the equipotential surface is large. This negative ion is strongly focused due to this large curvature of the equipotential surface. The negative ion is over-focused in the extractor, and then, divergent. Therefore, the difference of negative ion extraction location results in a geometrical aberration, and then the beam halo. Thus, the physical mechanism of the halo formation reported with the 2D PIC code [41,42] is verified with the 3D PIC modeling. Figure 35 shows the emittance diagrams of the H − ions in the extractor along the lines of the exits of (a) the PG, (b) the EXG, and (c) the ESG, respectively. From figure 35(a), it is indicated that there are the deviations (the curves in green) from the ellipse (the curve in blue). These deviations correspond to the beam halo, that is, the negative ions which are extracted from the periphery of the meniscus and have very large convergent angles. At the exit of the ESG, one can observe two components of the beam core and halo, which correspond to the ellipse and the deviations in figure 35(c).
From negative ion beam trajectories in the accelerator, the fractions of heat load due to the intercepted negative ions to the total beam power are listed in table 8. The fractions of the heat loads on the acceleration grids to the total negative ion beam power are evaluated from where k corresponds to the negative ion beamlets intercepted at each acceleration grid, and N is a total negative ion beamlet. For the comparison, the result with the 2D PIC [42] and the experimental result are also shown in table 8. It should be noted that in the present PIC simulations, the stripping loss of H − ions via the collision of H − ion and the residual H 2 gas molecules is not taken into account. However, the heat loads in the experiments are composed of not only the beam halo but also the secondary electrons [102]. Thus, the heat loads in the experimental result in table 8 are obtained by extrapolating the source gas pressure to 0 Pa, in which there is no contribution from the stripped electrons.
As seen from the comparison in table 8, both 2D and 3D model reproduce the basics experimental trend, i.e., the heat loads due to the beam halo become larger in the downstream grid. However, the heat loads on the A1G and the A2G are quantitatively better for the 3D PIC simulation result than those for the 2D PIC simulation result. As mentioned above, the penetration of the electric field into the source plasma is deeper for the infinite slit aperture in the 2D PIC simulation than more realistic cylindrical aperture in the 3D PIC simulation. The negative ions are focused largely due to the deeper penetration of the electric field. Therefore, in the 2D PIC simulation, the negative ions are likely to be intercepted at the EXG or the ESG instead of the A1G.

Summary and open problems
An overview of modeling studies of H − ion sources is given. Due to the recent progress of computer technologies, kinetic modeling of even large volume H − sources with the particle model is available and is indispensable in understanding H − source plasmas. Under close collaboration with the experimental studies, modeling studies contribute not only to the basic physics understanding, but also the improvement of the source performances.
In section 2, modeling studies of multi-cusp arc-discharge H − sources are reviewed. The modeling of the compact tandem type H − sources for medical applications [17,18,64,65] is first discussed with the main attention to the H − VP. Next, the H − SP under the Cs-seeded condition in the so-called 'giant H − ' sources for fusion applications are discussed. In such giant sources with multi-extraction apertures, the H − beam nonuniformity is one of the most critical issues to obtain high-power H − beams. A series of modeling studies [25,[27][28][29][30][31] largely contribute to the understanding of its origin and improvement of the beam uniformity. The following are the summary of the above modeling studies: • 3D analysis of the EEDF in the 3D real geometry and real magnetic field configuration of the arc-discharge H − sources is now possible by the KEIO-MARC code [26][27][28].
• The EEDF in arc-discharge plasmas has a non-equilibrium feature. The EEDF mainly consists of two components: the one is a low energy thermal component with Maxwellian distribution, the other is a high energy tail component.
• The non-equilibrium feature of the EEDF is very important to evaluate/predict the H − production amount in both the H − VP and H − SP.
• By 3D simulations of the tandem type H − source with the KEIO-MARC code successfully reproduce the magnetic filter effect.
• The high energy tail component of the EEDF in the driver region plays a key role in producing H 2 (v) with high vibrational excited number v, while the high energy component in the extraction region is drastically reduced. It should also be noted that recent progress in the modeling of H atom transport and collisional radiative (CR) model with the effect of the non-equilibrium feature of the EEDF contributes to the understanding of the origin of the beam non-uniformity in H − giant sources for fusion application; • In the H − SP, the localization of the high energy tail component of the EEDF by the magnetic drift plays a crucial role to understand the origin of the spatial non-uniformity of the H − SP in giant arc-discharge H − sources.
• More specifically, the local enhancement of H 2 dissociation in the driver region mainly by the high energy component of the EEDF (E>20 eV), and the resultant non-uniformity of the Franck-Condon H atom flux onto the PG is a possible origin of the H − beam non-uniformity. Direct comparisons of the calculated profiles of H α intensities with those in the QST (National Institute for Quantum and Radiological Science and Technology) 10A H − source experiments confirm above understanding, which was originally proposed in [25]. Recent progress in the H atom transport and CR model [30,35] in section 2.3 including the non-equilibrium features of the EEDF has made it possible to carry out such direct comparisons. The progress largely contributes to the source improvement. Namely, the following significant achievement has been obtained in the experiments by improving the magnetic configuration in H − source for JT-60SA N-NBI system at QST [22,23]; (1) beam uniformity is significantly improved (from 68% to 83%) over an area of the whole extraction area of 450×1100 mm 2 , and (2) the improvement of the beam uniformity leads to the production of 32 A H − ion beams with the whole extraction area.
In section 3, recent 2.5D EM-PIC-MCC (electromagnetic particle in cell with Monte-Carlo method for collisions) modeling [33][34][35] of the Linac4 RF-ICP source [8,93,94] has been reviewed. Direct comparisons of the time evolution of the Hα intensity in the capacitive coupling (E-mode) phase and the inductive coupling (Hmode) phase have been done [30,31]. H atom transport and CR model developed in section 2 is used with the EEDF obtained from the 2.5D EM-PIC-MCC modeling; • EM-PIC modeling reproduces the E-H transition observed in the experiments.
• Interesting behavior of the time evolution of the H α intensity observed in the E-mode phase, i.e., the temporal asymmetry of the H α intensity during one RF-cycle has been reproduced by the combined analysis of the 2.5D EM-PIC, H atom transport and CR model.
• The non-equilibrium feature of the EEDF in the E-mode phase plays a key role to understanding the temporal asymmetry of the H α intensity.
Although the model validation is limited to the Linac4 H − source at present, and further model validation and model improvement will be needed in the future, the 2.5D EM-PIC-MCC modeling developed here is useful in evaluating the change in the plasma loading impedance including the non-equilibrium feature of the EEDF before and after the E-H transition [95,96]. The predictive model of the time evolution of the plasma loading impedance is useful for the efficient coupling of the RF-system and plasma loading. Such a predictive model of the RF circuit equation and plasma loading impedance [91,92] with the plasma transport coefficients evaluated by the kinetic 2.5D EM-PIC-MCC modeling [95,96] is now started.
As for the extraction physics in section 4, 3D ES-PIC model [56,57] has recently been developed on the basis of a series of previous 2D ES-PIC modeling [36-41, 101, 104]. The present status and the important understanding obtained by these 2D and 3D modeling studies are summarized as follows; • The 2D PIC model is still useful for giving a basic insight into the extraction physics and to play a key role as pilot boats for the massive 3D PIC modeling.
• The formation mechanism of the ion-ion plasma proposed by the 2D PIC molding studies has been confirmed by the recent 3D modeling study in [56], i.e., the electron loss along the magnetic filter field is shown to be very important for ion-ion plasma formation. The thickness of the ion-ion plasma depends on the electron loss along the field line.
• Therefore, taking account of the large scale and non-local structure/topology of the magnetic field line is important for the 3D modeling with the single aperture model of the large H − ion source.
• Also, the beam halo mechanism proposed by the 2D ES-PIC model [41,42] is confirmed and essentially unchanged even with the 3D model [56,57]. Namely, the negative ions extracted from the periphery of the meniscus are over-focused in the extractor due to the electrostatic lens effect, and consequently become the beam halo. However, for a more quantitative comparison with the experiments and for a more design/engineering oriented study, 2D PIC modeling is far from the purpose, e.g., the fraction of the beam halo component is 60% by the 2D model, while 4.5% by the 3D model, as was shown in section 4.2. The former is too large, but the latter reasonably agrees with the experimental observation.
• Recently, 3D integrated modeling from the plasma meniscus formation to the beam acceleration [57] has been developed for more quantitative comparisons with experiments.
• The 3D integrated model not only reproduces the basic experimental trends of heat loads on the acceleration grids due to the beam halo, but also gives better agreement with the experiments in the quantitative sense.
Owing to the recent progress of integrated 3D modeling [56] from plasma meniscus formation to the beam acceleration, it is now possible to make more quantitative comparisons of the modeling results for the beam halo component and the beam quality, e.g., the beam emittance.
It is however still necessary to validate the assumptions and simplifications not only by one particular experiments, but also by various H − ion sources/experiments. For example, (i) effects of the reduced sizescaling on the beam quality, (ii) simple model for the electron loss along the magnetic field line, (iii) the production/destruction and the energy relaxation processes of the surface produced H − ions are not included in the model, (iv) the effect of the secondary electrons in the accelerator have been neglected and so on.
These validations and improvements are now underway step by step. For example, it is pointed out that the Coulomb collision play an important role in the energy relaxation process of H − ions and the resultant beam emittance, especially for low temperature extraction region and the low H 2 gas pressure case typical for the arcdischarge H − ion sources applied to the fusion plasma heating by the 2D PIC mode [105]. Recently, the effect of the Coulomb collision has been included in the 3D PIC model [106]. It has been confirmed by the 3D PIC model that the effect plays a significant role in the extraction mechanism for surface produced H − ions. Namely, surface produced H − ions are shown to be extracted not directly from the PG, but mostly from the bulk plasma region far from the PG due to the effect of Coulomb collisions. As was discussed in sections 1 and 4, a similar tendency has also been observed in the experiments in relation to the ion-ion plasma formation [51,52]. H − density far from the PG has significantly affected and reduced by applying extraction voltage. The effect of Coulomb collision is shown to be important for more quantitative and realistic evaluation of the total extracted H − current and the resultant H − beam optics.