3D-full wave and kinetics numerical modelling of electron cyclotron resonance ion sources plasma: steps towards self-consistency

Electron Cyclotron Resonance (ECR) ion Sources are the most performing machines for the production of intense beams of multi-charged ions in fundamental science, applied physics and industry. Investigation of plasma dynamics in ECRIS still remains a challenge. A better comprehension of electron heating, ionization and diffusion processes, ion confinement and ion beam formation is mandatory in order to increase ECRIS performances both in terms of output beams currents, charge states, beam quality (emittance minimization, beam halos suppression, etc.). Numerical solution of Vlasov equation via kinetic codes coupled to FEM solvers is ongoing at INFN-LNS, based on a PIC strategy. Preliminary results of the modeling will be shown about wave-plasma interaction and electron-ion confinement: the obtained results are very helpful to better understand the influence of the different parameters (especially RF frequency and power) on the ion beam formation mechanism.


Introduction
Among the various types of ion sources developed since 1950s, ECRIS are now the best candidates to support the growing request of intense beams of multicharged ions coming from both fundamental science (nuclear and particle Physics, especially) and applied research (neutrons spallation sources, subcritical nuclear reactors, hadrotherapy facilities). Inside an ECRIS machine [1,2] a dense (n e > 1 × 10 18 m −3 ) and hot (T > 1 keV) plasma, made of multicharged ions immersed in a dense cloud of energetic electrons, is confined by multi-Tesla magnetic fields and resonantly heated by some kWs of microwave power in the 2.45-28 GHz frequency range (2.45 GHz are used for under-tesla machines, while 28 GHz that is currently the state of the art require up to 4 T of maximum magnetic field). In the last three decades, ECR ion sources performances have been mainly improved following the semi-empirical Geller scaling laws [1]. These laws establish a direct (non-linear) relationship between the output currents and charge states, and the ion confinement time and plasma density; according to this simplified model, the formers can be increased by boosting the latters, i.e. by increasing the magnetic fields and the microwave heating frequency. However, this trend is now approaching saturation, due to technological issues and rising costs especially in the matter of superconducting magnets design and construction (further steps forward should require 5 T supporting 56 GHz of pumping frequency). Beside the path traced by Geller's laws, alternative approaches have been attempted by some groups spread over the world. The common idea was that the microwave-toplasma interaction was still not fully exploited in terms of energy deposition in the plasma core, and in the electron population playing the major role in multi-ionisation (i.e. the one lying in the 0.5-30 keV energy range). In fact, conversely to plasmas for fusion scenarios, in ECRIS the "first-pass" energy absorption in the plasma is deemed to be quite poor and cavity multireflections are needed in order to minimize the reflected RF power. Some laboratories have tried the option of multi-frequency heating [3,4], some others the way of superimposition of a discrete or continuous (broad) set of microwaves at different frequencies [5]. Another method that is deemed to further exploitation is the "frequency tuning effect" [4,6,7].
Due to the encouraging results coming from the frequency tuning method, the Catania's group developed a numerical code able to model the electron and ion dynamics in an ECR ion source using a Monte Carlo approach. This model allowed to explain for the first time the impact of the frequency tuning on the electron heating rapidity, on the plasma density distribution [8], and in turn on the shape, emittance and brightness of the output beam. Alternative models and numerical strategies have been attempted by other authors [9][10][11][12], that can be summarized in two groups: (a) assuming a "zero-dimensional" model for calculating the charge state distribution in a parametric way, or (b) considering a simplified magnetostatic scenario (only simple mirror configuration instead of minimum-B) for simulating the plasma dynamics via Fokker-Planck calculations based on a diffusion-like model for the RF heating [13]. Both these approaches do not include peculiarities of electromagnetic wave propagation into a dense, non-isotropic, non homogeneous plasma, that is indeed the condition holding in ECRIS-like systems. These properties are deemed to play a crucial role in frequency tuning effect, two frequency heating, etc. Actually, the challenging goal of predicting the electron/ion dynamics in the ECR source in a self-consistent way requires to consider the fundamental aspect of the coupling between the electromagnetic wave and the plasma. In this work we describe our approach to numerically solve the 3D Vlasov-Maxwell system of equations in a selfconsistent treatment.

Particle method
Our approach is based on a numerical particles-in-cell or "particle" method to approximate the distribution as a probability distribution function, according to the Klimontovich formalism [14]. The Vlasov theory for studies of plasma waves and wave-particle interactions starts from the mean-field Vlasov equation as kinetic model of a collisionless plasma with distribution function f α : with m α and q α mass and charge of the α species. We consider here a plasma consisting of a single species only (electrons), so that m e and e denote electron mass and charge, respectively. We assumed the "cold" plasma approximation (i.e. v φ v th , being v φ the wave's phase speed and v th the electron thermal speed) for modelling the wave propagation, and we retrieved the electromagnetic field through 3D numerical simulations using the FEM solver: COMSOL Multiphysics r . We do not use static electric field to compute electron motion; the latter, in fact, depends mostly on the RF electromagnetic field while electrostatic fields are dominant only in the sheath region, that is not important for the scopes of the present paper.
In the PIC approach to Vlasov equation we can consider a plasma as a collection of N macro-particles (being N much smaller than real plasma particles), with corresponding spatial coordinates and momenta described by the functions r i (t) and p i (t). These functions can be viewed as trajectories in the six-dimensional singleparticle phase space. The phase space density is thereby described by the Klimontovich distribution function: The trajectories obey the equations: since according to the Vlasov equation the particle distribution function is constant following the particle motion. Our Particle method simulates a plasma system by following a number of particle trajectories that obey the single particle equations of motion (3), including external RF field (electric and magnetic) and the superimposed magnetostatic field ensuring plasma confinement. Ion motion is not considered because we would like to explore the behaviour of the electrons assuming the ions as a background positive (and uniform) fluid. By the point of view of RF energy absorption, in fact, ions can be ignored considering their high inertia relatively to the high frequency domain. A specific routine of the developed computer code has been dedicated to particle motion (hereinafter called "particle mover routine", or simply "Particle Mover"): we used the Boris' scheme including relativistic electron mass variation [15].

Mean field definition
ECRIS plasmas are confined by a magnetostatic field (in the ∼ T range) obtained as the superposition of an axissymmetric field produced by Helmoltz coils and a radial field produced by a hexapole; overall equations describing B field components are: where S ex = B ext R 2 is a constant related to the gradient of the hexapolar field, B 1 = B inj = 4B ECR in the injection region, B 1 = B ext = 2B ECR in the extraction region, and B 2 = 0.8B ECR are related to the solenoid ones.
In this perspective, equation (1) can be linearised, splitting the distribution function into a stationary (or equilibrium) f e0 and an oscillating f e1 part: As an approximation of f e0 , we will initialize the particles position by a random sampling weighted on a presetted initial background density, while the velocities can be taken from a Maxwellian velocity distribution given by Gaussian distribution in the three axes with standard deviation σ = kTe me , assuming initial electron temperature T e = 1 eV. Then, the particles phase space evolves according to deterministic equations of motion (3).
The 3D PIC code is coupled to Maxwell's equations that, in the case of harmonic time dependence e jωt , read as: Displacement jω 0 E(r), and conduction J (r) currents are collected in a tensorial term describing the field-plasma interaction introducing the dielectric tensor: so that substituting (6) in (7) we can obtain the wave equation for the anisotropic plasma medium: that is a partial differential equation (PDE) including tensorial constitutive relations. Among the various options of numerical codes able to solve Maxwell equations, we used COMSOL Multiphysics for its suitability in modelling tensorial constitutive relations. The overall simulation strategy is depicted in Figure 1. The method is based on a strict interplay between COM-SOL solver and the particle mover implemented in MAT-LAB. Several loops are expected to be needed in order to get convergence. Due to huge time-consuming calculations needed at each loop-step k, hereby we show results up to step 2, i.e. after evaluation of RF action on electrons and vice-versa, for two times.
Our code does not assume any rotational symmetry for the magnetostatic field, so that the real 3D external B field (B-minimum) is taken into account, thus retrieving the permittivity tensor of the cold plasma approximation. In this latter case the determination of is based on the equation of motion for single particle: Solving (10) in v and using the constitutive relation: We can obtain σ that is related to as we can see in [16].
in which off-diagonal terms cannot be neglected.
In the tensor, with m = x, y, z, n = x, y, z and is the real part of relative permettivity r , is the imaginary part, ω the angular frequency of the microwave, ω p = nee 2 me 0 the plasma oscillation angular frequency, n e the electron density, m e the electron mass, e the electron charge, i the imaginary unit and ω eff the collision frequency; the latter accounts for the collision friction (thus modeling wave damping) and resolves the singularity of some elements of (12).
The numerical code (written in MATLAB r ) uses a numerical artifice to consider a well defined statistical sample and accumulate statistical information for all the life of the sample itself: in other terms, the "Particle Mover" code solves the equation of motion (3) of each single macro-particle for its entire life, accumulating the electron density in 3D grid. In this sense, our approach can be defined as a "stationary" PIC.
Particles move inside a 3D frame entering different cells: the plasma chamber domain is subdivided into 200 × 200 × 600 = 24 × 10 6 cells (so obtaining a submm 3 precision) thanks to a uniform spatial grid spacing Δx = Δy = Δz = 0.5 mm. In a time-integration view, the fields are calculated from an initial charge and current density, then the particles move for a small distance over a short Δt, finally the fields are recalculated according to the new phase space configuration, by repeating this procedure for N time steps. Conversely, we assume the stationary case where a sample of macroparticles followed for their entire life (i.e. until they impinge on the chamber walls, meaning they have been lost by plasma confinement system) describes the stationary structure of the plasma in the phase space, through the "trick" of density accumulation in a 3D grid; the density accumulation is made as single particles paths evolve in the volume contained into the plasma chamber. Indeed, we assume a stationary electromagnetic field as well.
With the simulation parameters listed in Table 1, we can follow 300 000 electrons/5 h, with an average electron life time T lif e = 27 μs. We also consider electrostatic (Spitzer) collisions at 90 • between electrons as a stochastic Poisson process, so that the time-dependent probability between each pair of consecutive events (electron collisions) has an exponential distribution with parameter 1/τ 90 • . Evaluating the position (therefore the density) inside the chamber, we can extract a collision probability  Table 2.

Simulation results
In ECRIS, electrons are resonantly heated by hundred of watts (or even few kW) of microwave power in the GHz range. RF is normally coupled to the plasma chamber through a rectangular waveguide. A simplified cylindrical chamber geometry, with a single on-axis microwave injection rectangular waveguide is shown Figure 2, as it was modelled for calculating RF field in our simulation.

Modelling ECRIS empty plasma chamber
As a first step k = 0 the electromagnetic field is solved in a vacuum chamber: the dimensions of the rectangular waveguide are 40.386 mm × 20.193 mm in order to allow a single T E 10 mode propagation at frequency of 5 GHz; the radius and length of the cylindrical plasma chamber are 50 mm and 300 mm, respectively: both geometrical dimension and frequency (lower than values normally used in last generation ECRIS) reduce computational costs without loss of generality in the approach validity, maintaining all the significant physical properties of a conventional ECRIS. The meshing sequence is controlled, setting a very low mesh size h ECR = 1 mm on the ECR layer and increasing it up to h = λRF 5 in the remaining computational domain, following the FEM "rule of thumb" of five nodes per wavelength, that represents an optimal compromise between RAM requirement and accuracy. Figure 3 shows the electric field distribution inside the cylindrical cavity excited by the rectangular waveguide. The presence of the cavity walls, modeled via the perfect electric conductor boundary condition, makes the electromagnetic field pattern to configure in a "cavity mode like" structure.

Electromagnetic field applied to the particles
Now the RF field should be introduced into the Particle Mover code, in order to obtain the phase-space stationary electron distribution. This first step k = 0, in which particles move in presence of the "vacuum" RF field, has been obtained elsewhere [8]. These former results highlight that electrons mostly concentrate into the iso-magnetic closed surface corresponding to the ECR condition. This high density volume was called "plasmoid",  while the surrounding volume was named "halo". In reference [8] results about 2D time-dependence PIC model were also mentioned, confirming this separation between inner/outer ECR surface plasma regions: a dense plasmoid surrounded by a rarefied halo. The result -although not self-consistently obtained by evaluating plasma effect on RF -was however in qualitative agreement with experimental results of Bibinov et al. [17] and Tuske et al. [18], they being among the few ones having performed spaceresolved measurements of plasma density through optical emission spectroscopy techniques [17,18]. Therefore, in the perspectives on the present paper -where we want to calculate plasma influence on RF and vice-versa -we can reasonably assume that the starting condition of the step k = 0 is given by a "plasmoid-halo" model. We then proceeded to calculate RF field once defined a plasmoid region with a maximum density located at the center of the chamber (z = 0), a "quasi-Gaussian" profile showing a flat-top shape into the plasmoid and a rapid decrease from ECR layers towards the plasma chamber walls, according to the following equation: where n c = ω 2 0 me e 2 is the cutoff density of the ordinary wave, n max = 0.7n c is the maximum electron density and n halo is the minimum density (n halo = 10 15 ), two order of magnitude lower than the n max = 10 17 inside the plasmoid.
In Figure 4 the 1D density profile is shown along with the corresponding magnetic field profile. According to the 1 We run the kinetic code on different MATLAB istances to compute the motion of many set of particles samples.  derived 3D spatial distribution of either n e and B, the tensor (12) can be evaluated at each mesh point. Figures 5  and 6 illustrate the COMSOL outputs for either electric field distribution and power dissipation (P diss = (σ·E)·E) into the cavity (and the plasma in particular). The profile of the ECR region has been put in evidence. It can be noted that the electric field distribution [19] is partially perturbed by the plasma, while keeping a disuniformity that still puts in evidence a standing wave nature of the inner-cavity radiation. The absolute strength of the electric field inside the plasma chamber is even increased, (the maximum value of the electric field inside the chamber is increased from 1.2 × 10 5 to 1.41 × 10 5 ), when the plasma is considered, meaning that in the specific case (i.e. for the considered frequency) the waveguide-to-cavity coupling has been improved by the plasma itself. Power deposition plot (Fig. 6), in addition, shows areas where the best coupling between the wave and the particles takes place (note that the figure is given in logarithmic colormap scale): the largest fraction of energy is obviously absorbed at the ECR (as expected), moderate energy exchange occurs in the inner plasmoid volume, while becoming several orders of magnitudes lower in the halo part of the chamber. The total input RF power was 1 kW, the 45% of which was reflected.
Once determined how does the simplified step k = 0 model of the plasma density profile affect the RF field, we can go back using the computed RF field as mean field (input field) of the Particle Mover. The electromagnetic field coming out from COMSOL mesh-node points is imported in Matlab, interpolated using the "mphinterp" function (built in "COMSOL with Matlab" library) that evaluates the electromagnetic field at the coordinates specified (we specified those of the space domain of the Particle Mover). Starting positions of the electrons are chosen according to a weighted sampling following from the density profile of Figure 4.
In Figure 7 the obtained 3D distribution of the electron density in the k = 1 step is shown (on the left), along with the density distribution projection over the plasmoid surface (i.e. ECR iso-magnetic surface) and the projected 2D density distribution in the transversal direction. The nonhomogeneous distribution of the density in the plasmoidhalo scheme is preserved even in step k = 1, meaning that it is a characteristic feature of this kind of plasmas. In reference [20] this structure was explained considering the plugging effect provided by the RF acceleration of electrons passing for the first time through the ECR, which are expelled at large rates from magnetic loss cones and whose turning points in the magnetic trap coincide with the ECR layer itself. An additional effect is probably given by pseudo-ponderomotive potential, as formerly supposed by Dimonte et al. [21], even though more detailed analysis are needed to evaluate its relative weight. Other nonhomogeneities are due to the magnetostatic field structure (the "arms" on the plamoid surface). Another characteristic feature is given by the "hole" in the near axis region, due to a density depletion. In order to go more in detail in understanding plasma structure coming out from the simulations, we splitted the above pictures according to different electron energy ranges [22].

Analisys of electron energy subdomains
We followed a commonly used convention in the ECRIS field, selecting five different ranges: I 1 : E < 2 eV, I 2 : 2 < E < 100 eV, I 3 : 100 < E < 1000 eV, I 4 : 1 < E < 50 keV, I 5 : E > 50 keVs, corresponding to the so-called "bulkplasma" electrons (I 1 ), to the electrons producing low charged ions (I 2 ), the ones generating ions at mean charge states (I 3 ), at high charge states (I 4 ) and, finally, to the so-called "suprathermal" electrons (I 5 ) whose role and formation mechanism is still controversial [23][24][25] but useless for ionization and highly charged ions build-up [23][24][25]. Density distributions in 3D views and in 2D plots integrated over the axial coordinate are shown in Figure 8.
No electrons populate the I 1 domain, meaning that interaction with RF heats the totality of the electrons bringing them in the higher energy domains. A large amount of electrons have few eV of energy and populate the I 2 interval: electrons in this range reach the highest relative density (>10 5 in a.u.) among all the selected intervals, but only in some localized regions. Their concentration assumes the maximum value in the near axis region.The distribution of the particles in the I 3 interval is smother than in the second one, but once again concentrated in the quasi-axial zone. Different is the situation of I 4 and I 5 domains: in these cases, the density concentrates mainly in off-axial regions, leaving a depletion hole all around the axis. Figures in the lower row of Figure 8 allow to compare relative densities among the various energetic domains, then comparing them with experimental data. Densities of electrons in I 4 and I 5 domains in near-axis region are about 10% and 1% of the I 2 and I 3 domains. Since we assumed an absolute plasmoid density value of 0.7n c , these numbers are in good agreement with the results collected and discussed in reference [26]: in that paper, authors performed an experiment measuring X-rays energies coming from near axis plasma in the 2-30 keV domain and in the >60 keV domain, that are roughly comparable with the simulated ones. A good agreement also comes out if simulated data are compared to the experiments performed by Biri et al. [22,27], which showed an off-axis concentration  of highly charged ions and hot electrons [22,27]. As already argued in reference [20], the hole formation could explain hollow beams generation often observed in ECRIS working below 18 GHz of pumping wave frequency [6,7]. In reference [20], however, calculations were stopped to step k = 0, and still controversial was the effect of the plasma on the RF field, also deemed eventually causing standing wave structure to disappear.
Step k = 1 calculations now confirm hollow density distribution into the plasma. The reason why standing wave structure distribution is invoked to explain hollow density/hollow beam formation can be understood by looking to Figures 9 and 10.
In those figures, it can be seen that the density distribution in the I 2 interval is affected by local "bunches" of particles roughly corresponding to regions of higher electromagnetic field (providing electron plugging). Less evident is the effect on the I 3 and I 4 domains (more energetic electrons are more affected by magnetic drifts which partially smooth the density distribution), but a near axis density depletion is evident in I 4 , i.e. in a rather high energetic domain. Figure 10 shows that on the top of plasmoid surface, corresponding to axial region, the RF field is one order of magnitude weaker than in off-axis zones. This may explain a lower plugging efficiency of electrons and the following axial density depletion. An even more clear picture of density distribution concerning the different energetic intervals is given by the 1D profiles shown in the sequence of Figures 11-15. In particular, Figure 11 displays the density distribution along a transversal axis and it shows that the higher is the energy, the far the density distributes from near axis regions. These outputs are in agreement with results shown in reference [20], confirming that the self-consistent solution is not so far. Concentration of the plasma inside the plasmoid is evident also along the longitudinal direction (see Fig. 12). Peaks are also placed in proximity of the ECR layers, where the electrons spend a relevant fraction of their "life" since mirrored by the combined effect of magnetic field force plus RF plugging. Figure 13 illustrates the azimuthal distribution (on the midplane of the plasma chamber) of the electrons (once again, displaced according to the different  energetic intervals). It is clear that the plasma distribution is modulated by the multi-mirror magnetic structure: electrons concentrate in pole regions, while being more rarefied in the gaps zones. This effect is less evident on the cold population (strong collisionality causes the density to spread out regardless of magnetic structure) and on the ultra-hot population (>50 keV), where the electrons are subjected to strong drifts across the magnetic field lines. It comes out that the electrons between 1 and 50 keV are the most affected by the magnetic field. Figures 14 and 15 highlight in more details the difference between poles and gaps. In near pole positions, the density is increased showing a peak in near resonance zone.
If the electric field causes the density to displace as discussed above, the frequency tuning should change the RF distribution, then the density structure. This finally should modify the ion beam shape, as actually observed in real experiments [6]. Next steps will try to investigate plasma structure dependence on fine tunings of the pumping wave frequency.

Conclusions
The proposed results can be considered as the first outputs from the novel strategy depicted in the paper, and on which fully self-consistent strategy can be based in the next future. The most important clues coming out   from the simulations are that although vacuum field RF field distribution (that is a cavity, modal field distribution) is perturbed by the plasma medium, the non uniformity in the electric field amplitude still persists in the plasma filled cavity. This non-uniformity can be correlated with non-uniform plasma distribution, explaining a number of experimental observations. A first attempt toward self-consistent, 3D simulation of an ECRIS plasma has been done, including the resonant behaviour of the plasma vessel. Although we are still far from strict convergence, some plasma features are now clear: the plasma concentrate mostly in near resonance region: a dense plasmoid is surrounded by a rarefied halo; the resonant cavity effects cannot be neglected: resonant mode structures modify plasma density distribution; electrons at different energies distribute differently in the space: cold electrons in the core, hot ones in near ECR regions; simulation results explain quite well the experimental observations.
Effects on the ion beam formation and beam properties can be argued: ions are formed were hot electrons are placed; this implies the ion beam shape depends on electron/ion distribution on the plasmoid surface; optimization of the ion beam formation [28] and transport must begin well inside the plasma.