A new numerical description of the interaction of an ion beam with a magnetized plasma in an ECR-based charge breeding device

The ion beam-plasma interaction is a relevant topic in several fields of plasma physics, from fusion devices to modern ion sources. This paper discusses the numerical modelling of the whole beam-plus-plasma-target system in case of 1+  ions entering an ECR-based charge breeder (ECR-CB). The model is able to reproduce the ion capture and the creation of the first charge states in the selected physics case, i.e. the interaction of a 85Rb1+ ions with the plasma of the 14.5 GHz PHOENIX ECR-CB installed at the Laboratoire de Physique Subatomique et de Cosmologie (LPSC) of Grenoble. The results show that a very narrow window of physical parameters for both the beam (energy and energy spread especially) and plasma (ion temperature, density, density structural distribution, self-generated ambipolar fields) exists which is able to reproduce very well the experimental results, providing an exhaustive picture of the involved phenomena. Possible non-linear interactions and the role played by the eventual onset of instabilities are also discussed.


Introduction
Starting from the first experiment performed in the 1980s at CRC Louvain-La-Neuve [1], the interest in post-accelerated radioactive ion beams (RIBs) has been growing constantly: this led to the construction of different ISOL (Isotope Separation On Line) facilities all over the world, like CARIBU [2], ISAC [3], SPIRAL [4], TRIAC (not operational any more) [5] and ISOLDE [6]. The European road map is now going through the development of second generation facilities pointing to EURISOL [7]: presently, three main projects are ongoing in different Laboratories, HIE-ISOLDE at CERN [8], SPIRAL2 at GANIL [9] and SPES at INFN-LNL [10].
All the above mentioned facilities adopt a charge breeder: mainly two techniques have been developed so far, consisting in adapting electron beam ion sources (EBIS) [11] or electron cyclotron resonance ion sources (ECRIS) [12] to the charge breeding process. Both techniques were compared during FP6 [13] in view of the performances expected for EURISOL, and were found to be complementary: among the two, the ECRIS-based charge breeders (ECR-CB) are usually chosen for their high acceptance in terms of injected current and emittance, the reliability, the possibility to operate both in continuous and pulsed mode, and the high charge states produced. A drawback is the presence of contaminants in the extracted beam, that could possibly hide the peaks of the radioactive ions of interest: this effect can be limited by properly treating all the surfaces exposed to vacuum, and by implementing a spectrometer of adequate resolution downstream the charge breeder, as in the case of the SPES project, which will adopt this kind of device [14,15].

A new numerical description of the interaction of an ion beam with a magnetized plasma in an ECR-based charge breeding device
In an ECR-CB, a magnetically confined plasma captures the incoming RIB through the long range i-i collisions, and then extract it as a highly charged ion beam, after charge multiplication through step-by-step ionization by energetic electrons. The effectiveness of this technique was proved first by the team of Professor Geller in the 1990s the first experiments of the → + + q 1 conversion technique employing an ECR source were performed with the ISOL-MAFIOS set-up, in the framework of the PIAFE project [16]. The results validated the methods and opened the way to its development: in the late 1990s the Laboratoire de Physique Subatomique et de Cosmologie (LPSC, but ISN at that time) build a new model, the PHOENIX booster [17], followed in the subsequent years by others developed by different Laboratories [18,19].
The performances of the ECR-CB nicely evolved in the past, however reaching recently a kind of saturation: presently, an important gap exists, between charge breeding of gaseous or metallic ions. By denoting with I q+ the extracted current on charge state q+, and I 1+ the 1+ current injected into the charge breeder, the global and single breeding efficiencies are defined as: Due to wall-recycling, ε g reaches 80% and 10% ⩽ ⩽ + ε q 15% for gaseous ions, while for metallic ones ε g decreases below 50%, while 5% ⩽ ⩽ + ε q 10% [20]: this discrepancy indicates the need of a better comprehension of the parameters affecting the capture efficiency (CB-efficiency), and numerical simulations are an effective way to fulfill this task.
A code aiming at modelling ECR-CB has to be particularly suitable to reproduce the elastic Coulomb collisions between injected and plasma ions; this issue is relevant also in other fields such as Astrophysics or inertial confinement fusion [24][25][26][27]; electron dynamics in high density processing plasmas also require advanced Coulomb collisions simulations, even in those cases where 90° scattering is dominated by e-neutral collisions: e-e collisions determine in fact the electron energy distribution function, populating the high-energy tail by energy up-scattering. Of particular relevance is the application of similar simulation strategies to ECR-CB: RIBs intensities are in fact usually orders of magnitude lower than stable beams, forcing to 'blindly' tune the post-accelerator [28]. Another relevant issue consists in the response of the buffer plasma to the injection, and further thermalization, of the incoming beam. This problem resembles the well-known plasma-beam system which is prone to the development of a number of instabilities. Despite the intensity of the incoming beam is usually considered as being negligible if compared with the energetic content (density and temper ature) of the accepting plasma, signatures of a strong interplay have been experimentally found, consisting in the depletion of the highly charged ions in the background plasma [29].
The paper presents therefore a new numerical approach, developed in a MatLab environment, for the description of the interaction and following capture of a 1+ ion beam by the plasma of an ECR-CB. The code implements a formalism based on the Langevin equation [21] and makes a significant step forward with respect to previous papers adopting different approaches [22,23]: the most relevant advancements consist in a systematic analysis of the relative importance of the injected beam and plasma parameters, making the results as directly comparable with experiments, and, moreover, in the adopted fully 3D plasma model, with several details as concern its structure. In the last part of the paper, after the description of the numerical method, and the comparison with the experimental results, some candidate mechanisms related to beam-plasma instabilities will be described in details. As it will be described along the paper, simulations results not only show a very good agreement with the theoretical predictions, but reproduce also experimental evidences obtained on the LPSC test bench with the PHOENIX charge breeder [29]. The results of the numerical simulations will be useful for the future operation of SPES, which will adopt an upgraded version of the PHOENIX charge breeder [14,15].

Modeling strategy
In an ECR-CB, particles extracted as a q+ ion beam are fed as an external 1+ beam, injected at an energy E = eV 1+ (in the order of tens of keV) into the charge breeder held at a potential As soon as the 1+ ions approach the charge breeder, their energy decreases due to V CB : by properly regulating the ∆V, 1+ ions can overcome the maximum of the magnetic field and the plasma potential thus penetrating the plasma core, where they undergo multiple long-range Coulomb collision with plasma ions. The effect of such interaction is twofold: on one hand, there is a general loss of directed velocity; on the other hand, there is diffusion in velocity space, leading to a distribution that has the same thermal spread as the background plasma [30].
The collisional slowing down under the influence of an inverse-square force was formerly treated by Chandrasekhar, describing the interaction of the stars in the gravitational field [31,32]; the mathematical treatment was then used by Spitzer to describe the interaction of charged particles with a plasma [33]. Coulomb collisions are a collective effect for which we can deduce some relaxation times [34]: the code described in this paper implements this process in a single particle approach, by solving the ions equation of motion under the influence of external electromagnetic forces, taking into account the presence of a plasma with a given structure, density and temperature. This approach allows, among the other things, to store and track single particles trajectories: this is important even from the point of view of radiation-protection, because the zones where particles are lost become points where radioactivity accumulates, whose position has to be known in case of unexpected maintenance. More, it is possible to implement various kind of collisions, predicting average behavior for particles motion and lifetime; additionally, it allows to figure out an explaination for the experimentally observed destruction of plasma high charge states due to the injected 1+ beam [29].

Background theory
The mechanism of Coulomb collisions can be described by considering a test particle p of mass m, charge Ze and velocity v, being part of an ensemble distributed following any , , , moving in a background of field particles of mass m s and charge ′ Z e, whose distribution function is ( , , s s s . The long-range Coulomb forces between the test particle and the field particles will cause p to experience a multiplicity of distant collisions, far more numerous than close ones, so that all the changes in direction and speed experienced by p will be small. The evolution of f is described by the Fokker-Plank equation = C f D that, following the formalism introduced by McDonald and Rosenbluth [35], can be expressed using the so called Superpotentials H and G : being Λ ln the Coulomb logarithm and /( ) µ = + mm m m s s the reduced mass. The two terms multiplying f in parenthesis in equation (3) represent the friction and diffusion coefficients: the former takes into account the average decrease of the particle velocity; the latter reflects the spread in the velocity space. In case of different kinds of field particles, a separate kinetic equation is required for each species.
The coefficients introduced above can be made explicit if one considers as field particles plasma ions with density n s and in thermal equilibrium at a given temperature KT i : in this case, considering the plasma as an infinite bath, the final effect of the interaction is the energy equilibrium between injected ions and plasma ions at the bath temperature, whatever the initial particle distribution be. Introducing the variable and v is the test particle speed, the coefficients can be expressed by the equations: x x x 2 2 : following the terminology of Chandrasekhar, equation (7) expresses the coefficient of dynamical friction, while equations (8) 1 and (8) 2 the parallel and perpendicular diffusion coefficients (parallel and perpendicular are referred to the direction of v). The function G(x) shows a maximum at x = 1, while ( ) Φ x increases very fast up to x = 2, and then starts saturating: this means that, for a given injected particle and plasma, the friction and parallel diffusion coefficients show a maximum when the particle velocity is equal to the most probable thermal speed, while the transversal diffusion always increases with the particle velocity.
From the coefficients defined above three relaxation times can be derived, the slowing down time τ s , the deflection time τ D and the energy exchange time τ E : where the second expression in the last equation follows considering equation (8) 1 and the approximation ( . When the test and the scatterers are identical and in thermal equilibrium there is no friction or energy exchange: in this case, equation (10) evaluated at = x 1.5 gives the well known Spitzer collision time [33]. In ECR ion sources, this time gives a measure of how fast charged particles (especially cold electrons) are scattered by Coulomb collisions into the velocity loss cone, and then lost from the magnetic trap [12]: a useful handy formula is˜τ where τ c is in s, A is in amu ( / = A 1 1836 e ), K T is in eV and ñ in cm −3 .

Case study
The validity of the proposed modelling strategy has been cross-checked by the results collected during an experimental campaign performed at LPSC within the EMILIE project [29,38]: around 350 nA of Rb 1+ ions were injected inside the charge breeder, obtaining maximum efficiencies ⩾ + ε 5% q for q+ = 1+, 18+ and 20+, and < ε 50% g , at an optimun ∆ ∼ − V 12 opt V. Then, by acquiring the so called ∆V curves (the CB-efficiency of a given ion as a function of ∆V) at different power levels, completely different trends were observed for low and high charge states: an example is shown in figure 1(a). While the CB-efficiency of Rb 20+ ions (and any q+ >5+) shows the typical peak at ∆V opt , the one of Rb 1+ ions sees a monotonic increase up to more or less ∆ = − V 25 V, where it starts saturating. The same trend is observed decreasing the microwave power feeding the ECR-CB, as part (b) of figure 1 shows.
A first interpretation of the observed discrepancy was that the recorded Rb 1+ ions were basically injected particles only weakly interacting with the plasma (hence, not captured). For a rough benchmark of this hypothesis, the charge breeding time τ CB of Rb 1+ ions (at the optimum ∆V) was compared to their 'flight time' through the charge breeder, switching off the plasma but keeping the magnetic field on: in both cases the measured value was around 500 μs. A similar scenario was then arranged in order to benchmark the numerical code: the interaction of a Rb 1+ ion beam with the plasma of the PHOENIX charge breeder was simulated, trying to reproduce the CB-efficiency curves in figure 1(b). Considering the measured value for τ CB , the simulation time was fixed to = T 500 span μs: as a further condition for validation we also asked for a global capture ⩾ ε 40% g at ∆V opt , as measured during the experiments.

Input parameters
The domain of the simulation was set inside the PHOENIX plasma chamber [36,37], as shown in figure 2: it is basically a cylinder with a radius a = 0.036 m, and a length extending from the maximum of the magnetic field at injection to the extraction hole (l = 0.248 m).
In this region, a minimum-B magnetic structure is created by three coils and a permanent magnet exapole: the maximum at injection is T. By setting the origin of a coordinates system in the middle of the central coil, the magnetic field on axis is very well interpolated by a 6th order polynomial b(z) (R 2 = 0.9999), allowing us to retrieve the other components of the magnetic field in any point of the domain, as in the following formulas:  giving the field in T with x, y, z in m and R = 617.28 T m −2 . Rb 1+ ions enter the simulation domain at the axial position corresponding to B inj (z = −0.148 m): in order to have realistic starting conditions, the transport up to the charge breeder of a 85 Rb 1+ beam was simulated first, using the tracking code SIMION [39] and the actual geometry of the LPSC beam line (including the fringe field of the charge breeder magnetic trap). The scheme of the beam line as implemented in SIMION is shown in figure 3: 85 Rb 1+ ions at an energy E = eV 1+ = 20 keV are injected into the charge breeder at a potential V CB , through a double Einzel lens (V V , 1 2 ). To take into account the deceleration due to the plasma potential, the area inside the charge breeder, extending from the maximum at injection till the end of the geometry, can be set to a given potential V P higher than V CB : particles positions and velocities are extracted at the beginning of the plasma area (vertical red line in figure 3), and become the starting conditions for the MatLab code. The only drawback here is that the plasma is seen as a fixed equipotential surface, without any sheath effect.
With this injection scheme, particles arrive at the plasma boundary with an energy P , so the simulated Rb 1+ ΔV curves can be expressed as a function of CB P (the plasma is at positive potential with respect to the charge breeder), this last quantity is more negative with respect to ∆V inj . To take into account this effect, simulated curves will be allowed to shift towards lower (more negative) ∆V in order to find the best agreement with the experimental ones: this could be a measure of the plasma potential with respect to the charge breeder, even if the verification of this parameter goes beyond the aims of this code.

Modelling strategy
The conventional technique to implement collisional processes in numerical simulations is the Monte Carlo approach. The long-range Coulomb collisions consist in many, almost continuous, very small angle scatterings: to treat them as a succession of random effects would require a numerically unmanageable integration steps. An alternative approach, which has been emphasized in the analytical development of plasma kinetic theory, is to represent Coulomb scattering through a Fokker-Planck equation, that in the present case takes the form of equation (3): direct numerical solutions of this equation are often performed, but there is no obvious way to combine this procedure with a PIC (particle-in-cell) or single particle simulation. In the case of a charged particle moving through a plasma in thermal equilibrium, Manheimer et al [40] showed that the problem can be overcome by employing the Langevin equation: the formalism is based on the theoretical analysis carried out by Chandrasekhar on the Brownian motion [21]. The evolution of the particle velocity can be described by a Langevin equation which, to first order accuracy in ∆t, is entirely equivalent to the Fokker-Planck equation expressed by equation (3). The variation ∆v c of the test particle velocity ) in a small time interval ∆t, due only to Coulomb collisions with a plasma in thermal equilibrium, can be then expressed as: that is, as the sum of a systematic part and a fluctuating part, which is characteristic of the Brownian motion. In the above formula (17) takes into account for the deterministic dynamical friction (ν τ = − s s 1 , see equation (9)), and is always directed along the original velocity, while v rand represents the results of the stochastic 'kicks' suffered by the particle. The 3D components of such vector are distributed according to the formula: can be expressed through equations (8) 1 and (8) 2 . Here v 3 lies along the original velocity, while the components v 1 and v 2 are oriented in two directions perpendicular between each other and to v 3 . By expressing the velocity at each time step through equations (16)÷(18), the Langevin equation ensures that the particle evolves towards a Maxwell-Boltzmann distribution.

Preliminary benchmarks
As a first step, equation (16) was benchmarked in a simplified scenario which consisted in a beam of monochromatic charged particles interacting with a 'single-specie' uniform plasma of infinite extent: in order to keep reasonably low the simulation time, we considered 10 000 ions with mass number A = 132 (representing the main radioactive element will be produced at SPES, 132 Sn), velocity components = = v v 0  To keep reasonably low the simulation time, for this particular case T span was fixed to 288 μs, still much longer than any characteristic time: at each time step, fixed to ∆ = − t 10 10 s, and for each particle, the code calculates ν s (so ∆v fric ), ⊥ D and ∥ D ; then, it individuates the three possible directions 1, 2 and 3, and construct the vector v rand through equation (18). The results revealed an effective thermalization at the expected temperature, fulfilling the characteristic times, as shown in figure 4. The velocity drops to 1/e of its initial value after about 13 μs, in very good agreement with the estimated value of τ s , then decreasing to 0 in around 50-60 μs; after that time, it only slightly fluctuates till the end of the simulation. Figure 5 shows the average energy for each degree of freedom reaches and keeps the expected value / KT 2 i , which is mainteined till the end of the simulation time. The isotropy in the velocity space was also checked: the standard deviation for each spatial direction, shown in figure 6, very soon start to fluctuate around the expected value of 851.3 m s −1 , in perfect agreement with the theoretical estimates. The statistical consistency of the injected ions velocity distribution with a Maxwell-Boltzmann distribution was checked by using a Matlab routine named normplot: consistency plots are given in figure 7 showing on the left a gaussian distribution along z as given by 10 000 random numbers extracted from a normal distribution with mean 0 and σ = 851.3, while on the right the results coming out from the simulations. For each velocity component, the crosscheck revealed a very good agreement, thereby confirming that the Maxwell-Boltzmann distribution is reached during the thermalization process. The computational accuracy was also checked by changing the temporal integration step. When reducing it from 10 −10 down to 10 −11 no significant changes where observed in the beam speed relaxation time nor in the other simulation outcomes. In order to check the robustness of the beam thermalization model, we also tried to simulate the heating of an 'under-thermal' beam of ions due to the interaction with the background plasma. This time the ions started with the 'under-thermal' velocity ∼ v 340 m s −1 ; expected characteristic times were τ = 9.31 s μs, τ = 0.41 D μs and τ = 0.21 E μs. Figure 8 shows the code is able to effectively describe the expected phenomenon, and the slow particles are heated by interacting with the plasma: the average energy increases up to the expected value. After this initial model benchmark, the ECR-CB plasma was implemented in two step of complexity, as described in the following sections.

Basic plasma-target model (BPM).
Several efforts were dedicated to the modelling of the plasma. We proceeded by steps of increasing complexity. For the first group of simulations, we implemented what we called a basic plasma model (BPM), where ambipolar potentials and ioniz ations of the incoming ions were not taken into account. Following the 'plasmoid-halo' scheme, proposed by theor etical papers [41,42] and experimentally verified, the plasma was assumed as subdivided in two zones: the core, for ⩽ B B ecr , with a density = n n core , and the halo, for > B B ecr , with / = = n n n 100 halo core . The plasma was made by oxygen excited at 14.521 GHz ( ∼ B 0.52 ecr T), with a given temperature KT i (several values where tested) and an average charge ⟨ ⟩ = z 2.5. Despite being over-simplified, this model revealed useful information, especially when attempting to reproduce the experimentally argued picture of 'not capture-not ionized' injected Rb 1+ ions.
Starting from the position z = −0.148 m, the code integrates in a for loop the equation of motion for N ions with a time step ∆t = 10 −10 s, calculating the time variation of the velocity due to the static magnetic field and the Coulomb collisions. At each time step it basically solves the equation: Lorentz is the contribution due to the Lorentz force, calculated through the Boris method [43]. As prevously said, the integration time was fixed to = T 500 span μs, for a total of 5·10 +6 iterations: each 5 μs the code saves the entire work space on a file, in order to postprocess the data. After each iteration, it also checks if particles went out the domain of the simulation, and in case it removes them from the calculation, storing their positions and velocities in a matrix. For the particular geometry that was used,  m.

Refined plasma model (RPM).
To let the simulated plasma be closer to reality, in a second stage we implemented a more realistic plasma model (RPM), including plasma self-generated ambipolar potentials and ionizations. It is widely known that the confinement of highly charged ions is improved by the presence of a negative potential dip in the plasmoid [44]. Its origin is still controversial: some theories argue it is due to the collisionless warm electrons component that, adiabatically confined by the magnetic field, forms a negatively charged cloud, then a potential φ ∆ , which in turn improves the ion trapping [12]. Interesting numerical studies, recently carried out at LNS on the plasma of the ECR source SERSE, attributed the origin of the potential dip to a so-called  'double layer' established at the edge of the plasmoid [45]; it was found also that the depth of the potential dip could be heuristically estimated as The 'potential dip' is seen by the simulated particles through its electric field components, that have to be included in the part of equation (19) due to the Lorentz force. Due to the lack of the real potential map for the PHOENIX CB, we rescaled and re-adapted the one obtained from the quasiself-consistent simulations performed for the SERSE source of INFN-LNS, expressing it as a function of the ratio / B B ecr , then rescaling to the magnetic configuration of PHOENIX. To obtain 3D maps for the electric field components ( ) E E E , , x y z , the potential values were first calculated and stored on a matrix coinciding with the spatial 3D mesh of the simulation domain (mesh-size element was 1 mm). Figure 9 shows the modulus of the electric field on the plane y = 0, assuming = KT 1 i eV. Ionizations are included in the calculation as a step-by-step process thorugh 3D maps of single ionization times → ( ) τ + i i ion 1 for charge states from 1+ to 30+, calculated using the well known semi-empirical Lotz formula [46,47]: In this RPM the warm electron component has the same spatial distribution as the cold one, with / = n n 10 warm : the warm electron temperature KT e warm was fixed to 1 keV inside the core, and 100 eV in the halo. After executing all the calculation sequences performed even for the BPM, in the RPM the code selects the given ionization time depending on the charge state of the particle (all are 1+ at the beginning), and turns it into an ionization probability Then, it checks if ionization occurred by applying the Monte Carlo technique, thus eventually increases the charge state of the particle of one unity.

BPM: results
For the first set of simulations the plasma temperature was fixed at   m, where r ext is the extraction hole. In view of the comparison with experimental results, the reported values take into account for the 80% of transmission of the LPSC beam line. This first set of simulations revealed a very low capture efficiency, almost independent from the plasma density, weakly from the injection energy: appreciable variations concern only the relative losses towards injection, radial or extraction sides. This picture is clearer when looking at part (a) of figure 10, showing the percentage of captured particles as a function of the plasma denisty for = − ∆ = E e V 12 inj opt eV: it is basically constant by varying the density, except for the lowest one for which the interaction with the plasma can be considered weaker. It is now interesting to evaluate the expected residence time of the injected particles in thermal equilibrium. In a magnetized plasma, the most probable ion dynamical confinement mechanism depends on the collisional/collisionless condition in which the ions can be found into the plasmoid: in other words, the Coulomb collision frequency ν coll has to be compared with their Larmor frequency ω c [45]. The condition ν ω > c coll holds for a collisional regime, vice versa the plasma is said to be magnetized. Thermalized ions do not exchange energy with the plasma, so that only the 90° diffusion influences ν coll (equation (10) evaluated for ). The trend of the collision and of the Larmor frequencies according with the assumed density and magnetic field profiles is shown in figure 11. It can be seen that, for any value of the plasma density and in any point along the axis, the condition ⩽ ν ω c coll holds, so the particles are magnetized: the confinement time can be estimated in this case from the formula where R is the mirror ratio, l the length of the plasma and v T is the one dimensional thermal speed. The previous relation  perfectly explains the results of the simulations: for a given magnetic trap and plasma length, the ion confinement time does not depend on the plasma density, in agreement with what can be observed in figure 10. More, injected ions finally reach thermal equilibrium (except for the lowest density) so, for a given KT i , all of them acquire the same value v T , explaining the weak dependence of the loss rate from the injection energy observed in table 1. Evaluating the confinement time for the PHOENIX charge breeder (R = 1.5-2.3, l = 0.288 m), it was found to be between 407 μs and 624 μs, in very good agreement with the loss rate obtained after 500 μs.
The following step was to investigate the influence of the ion temperature. Starting from an experimental ∆ = − V 12 opt V (see section 2.2), and supposing a plasma potential (with respect to the charge breeder) of 10 V, this translated in an optimum injection energy of 2 eV for Rb 1+ ions. By accepting the idea that the optimim injection energy corresponds to the maximum of the friction and longitudinal diffusion coefficiencts (i.e. x = v/C s = 1, see equations (7) and (8) 2 ), KT i was decreased to 0.376 eV: a higher influence of the plasma (due to the higher Coulomb collision frequency ν coll ) and a better confinement (due to the decrase of v T ) were expected. The results are summarized in table 2: the capture efficiency in general increased, showing now a different trend depending on the plasma density and the injection energy. This time, when looking to the captured particles fraction versus the plasma density (see part (b) of figure 10, again for = − ∆ = E e V 12 inj opt eV), the influence of the plasma density is evident, in addition to the positive effect of the lower ion temperature; for the highest densities, the capture reaches values close to the experimental ones. By checking again the coillisional/collisionless condition inside the plasmaoid it can be shown that, by increasing the density, the plasma passes from a magnetic (up to = * n n 0.3 co ) to a collisional regime (for the highest densities). In this last case the plasma behaves like a viscous medium: its ability to trap particles is due to the 'zig-zag' paths that are the consequence of very frequent collisions; in this case, the collision dynamics prevail on the magn etic field and ions are said 'unmagnezited'. A handy form ula to evaluate the confinement time of an ion of charge q in the collisional regime is given by the following equation [44]: where l is the plasmoid length in cm, Λ ∼ ln 15, A is the plasma mass number, ⟨ ⟩ z is its average charge, n the density in cm −3 and KT i is in eV. This formula is commonly applied for evaluating the residence times of highly charged ions in ECR plasmas. Considering as plasmoid length the distance between the two resonances on axis (l = 0.125 m), the above formula gives very high confinement times for the simulated plasma, in the order of tens of ms. For heavy elements like Rubidium, it is difficult to apply equation (21) just for only Rb 1+ ions. However, what can be expected is that after 500 μs a consistent number of particles are still inside the plasma chamber, as confirmed by the simulations in the higher densities regimes. At the lower densities, where the magnetic regime holds, the confinement time evaluated for the low ion temperature case ranges from 0.664 to 1.000 ms, thus validating the values of the simulated total capture and the differences with respect to the high temperature case.
In view of the comparison with the experimental results, this simple model gave interesting information: • In order to reach experimental values for the ion capture, the plasma ion temperature has to be low • At both the plasma temperatures, the density has to be decreased at least to = * n n 0.3 co in order to extract Rb 1+ ions Figure 12 compares the experimental 85 Rb 1+ efficiencies with the results of the simulations: the calculated data have similar In summary, we can say that the BPM in particular put in evidence the role played by the ion temperature: if the temperature is high, particles motion is controlled by the magnetic field and the confinement times are not consistent with the formation of high charge states trough step-by-step ionizations. By lowering the temperature, the overall confinement improves: for the higher densities the plasma shows a collisional regime and the total capture reaches values comparable with the experimental ones, even if no Rb 1+ ions are

RPM: results
The second stage of simulatins, in the frame of the refined plasma-target model, was initially based on a comparison with the previous model for the same input parameters, such as = KT 1 i eV and = n n co . The combined effect of potential dip and ionisations is clearly visible in figure 13 part (a), where the percentage of captured particles is shown as a function of the injection energy. The build up of weakly charged ions helps the trapping due to the ambipolar potential, whose influence is weighted by the factor ( / ) φ − ∆ q KT exp i : figure 13(b) shows the charge state distribution obtained at the end of the simulation for E inj = − ∆ = e V 12 opt eV. The formation of charge states higher than 1+ is consistent with the ioniz ation times evaluated through the Lotz formula. The gap between the two models when looking to the captured particles fraction increases with the injection energy, even if neither with the RPM acceptable values for the global capture are reached at this ion temperature; in addition, no Rb 1+ ions are extracted, meaning that the simulated plasma parameters are still not fully appropriate to reproduce the trend of figure 1(b). Evidently, the influence of the ion temperature on the capture process is even stronger than the combined effect of ionizations and ambipolar potentials.
Therefore a set of simulations was run in order to find the correct parameters able to reproduce the two curves of figure 1(b). The condition for which Rb 1+ ions are able to come out from the extraction hole-still having an acceptable global capture-is that both the plasma density and ion temperature be reasonably low. With = * n n 0.3 core co and = KT 0.376 i eV, both the global capture and the Rb 1+ CB-efficiency were still far from reality. The plasma density was then lowered to = * n n 0.1 core co . The results are shown in table 3: now, for the lowest energies, the global capture is comparable with experimental values. The simulated Rb 1+ efficiencies start to agree with the experimental curves, as shown in figure 14: the two trends in fact coincide with only a slight shift of the simulated values of around −2.5 V; anyway, for energies around ∆V opt , the calculated total efficiency is still too low.
The estimation of ⟨ ⟩ z was refined according to an experimentally obtained Oxygen spectrum acquired during Rb experiments, being ⟨ ⟩ = z 3; then, a realistic energy spread of 2 eV was considered for the injected beam. The ion temperature was again slightly decreased to = KT 0.3 i eV, simulating the densities = * n n 0.1 core co and * n 0.075 co . The outcomes of these last simulations are shown in tables 4 and 5: in bold the injection energies for which the simulations give values for the global capture and the Rb 1+ CB-efficiency in agreement with experimental results at ∆V opt . The comparison with the experimental curve at higher power in figure 15, parts (a) and (b), says that now the agreement is very good. : the first comparison shows a better agreement with the experimental curve, but the necessary shift is too low to be considered as a measure of the plasma potential; the second comparison was obtained with a shift consistent with a realistic value of the plasma potential, but exibit a small discrepancy at high injection energies, that can anyway be explained by the fact that for > Λ x ln 2 the evaluation of the diffusion  coefficient ⟨ ⟩ ∥ ∆v starts loosing accuracy. We also tried to reproduce the experimental trend at lower microwave power (see figure 1(b)): assuming a primary effect of the decreased RF power on the plasma density, we carried out another simulation keeping the same parameters except the plasma density, decreased to = * n n 0.05 core co . Figure 15 part (c) shows the model is again in very good agreement with the experiments (the curve was shifted of the same energy as the case of figure 15(b)).

Discussion
It is interesting now to analyse the information obtained when post-processing the results of the simulation with = E 7 inj eV and = * n n 0.075 core co . Considering the applied shift, this should correspond to an experimental ∆V of −11 V. Figure 16 shows the spatial distribution of the captured particles at the end of the simulations, with the surface corresponding to = B B ecr in blue: it can be noted that most of them are trapped inside the plasmoid, while a few stay in the halo; the track left by the particles during the injection process is also visible. Analyzing the velocitiy distribution of the confined particles, it was found that a Maxwell-Boltzmann distribution was reached, with ( / ) / σ = KT M Rb i 1 2 for each degree of freedom. The final charge state distribution in figure 17 shows the effect of ionizations: while the 62% of the captured ions is still charged to 1+, the 32%, 5% and 1% are, respectively, charged 2+, 3+ and 4+, in agreement with the ionization times estimated from equation (20).

Analysis of the ballistic phase
The ballistics of the process shed ligth on an important aspect of the ECRIS-based charge breeding. In the scientific community it is widely accepted that ∆V opt corresponds, net of the plasma potential, to the maximum of the function G(x) that describes the friction and parallel diffusion coefficients, that is x = 1. Considering Rubidium ions injected in an Oxygen plasma, this value for x translates in optimum injection energies of 5.3 eV, 2 eV and 1.6 eV for = KT 1 i eV, 0.376 eV and 0.3 eV. From the experimental value ∆ ∼ − V 12 opt V, it follows that the plasma potential should be, respectively, lower than 7 V, 10 V and 10.4 V. Simulations reveal instead that the optimum injection energy for Rubidium ions is up to four times higher than expected, as if an extra energy should be necessary for the ions to adequately penetrate the plasmoid, then being captured. This pushes for a lower estimate of the plasma potential (4 V in this case), that is also in agreement with more recent experimental results of charge breeding of  Rubidium obtained with the SPES charge breeder at LPSC [14]. The concept of optimum injection energy is perfectly clarified looking at figure 18, showing in the middle the xz projection of the particles density map obtained for = E 7 inj eV. Injected particles leave a 'density track' along the path toward the plasma then, once inside the plasmoid (indicated by the contour of the resonance surface in white, they are mostly confined on a narrow area around the axis. The situation is different for the other injection energies, as shown in the rest of figure 18: depending on the particular value, particles are stopped before reaching the middle of the plasmoid (with the risk of being lost back at injection due to diffusion), or some of them leak towards the extraction, loosing the possibility of being ionized.

The effect of the ion beam on the magnetized plasma
Another important aspect, mentioned in section 2, concerns the evident effects on the extracted charge state distribution of a so low (compared to the plasma density) injected intensity. It has been experimentally verified that the charge state distribution of the plasma ions is greatly modified by the injection of a 1+ beam: as an example, figure 19 shows the effect on an Oxygen plasma, injecting almost 1 μA of Xe 1+ ions inside the SPES ECR-CB. Looking at the spectrum, it is clear that the number of Oxygen ions extracted decreases by injecting the 1+ beam. A possible explaination of this effect could be given in terms of ion heating due to ion-ion col lisions, causing enhanced ions diffusion and a reduction of their lifetime. To prove this hypothesis, we evaluated the 'energy release maps', that is the spatial 3D map of the energy lost by the injected particles. In figure 20 it can be seen that basically the injected ions don't loose energy until they enter the plasmoid; once there, most of their energy is locally released on a very thin 'plume' around the axis. If this local release is faster than the capatibility of plasma ions to exchange energy between themselves, then a local heating is possible, with a consequent decrease of the ion confinement time. The average energy per unit volume of the plasma is n KT 3 2 i i , and it is exchanged in a time τ c , expressed by equation (12). For this particular case, we have = ⋅ n 6.5 10     eV·mm −3 ·s −1 . This value is much lower than H, so it seems to be excluded a direct ion heating causing the depletion of highly charged ions as depicted in figure 19. The picture changes a little bit if one restricts the analysis to only those points where the energy released is, respectively, higher than 10% and 20% of the maximum: the rates calculated in this cases are, respectively, = ⋅ H 2.9 10 9 eV·mm −3 ·s −1 and = ⋅ H 4.8 10 9 eV·mm −3 ·s −1 . It can be observed that the value of H can be locally more than one order of magnitude higher than what calculated on average, even if not enough to support the theory of direct ion heating.
According to the obtained results, the ion heating could be an indirect consequence of a plasma, temporal-growing instability triggered by the beam-plasma interaction. It has been shown by Neufeld and Wright [48] that an ion beam, propagating with any velocity inside a plasma, along the direction of a magnetic field B 0 , can excite electromagnetic waves. In particular, those waves have a circular polarization (in resonance with plasma electrons) and propagate along the magn etic field: the range of the excited ω i is wide, and can extend to ω ∼ Ω i e , where / Ω = eB m e 0 e is the electron cycloton frequency. Depending on the value of / β = v c (where v is the ion speed, while c is the speed of light), the beam can excite up to three differente waves, two of which are usually of low frequency (compared to Ω e ): in the present case we have β 1, and the theory forsees that only one wave is excited, at a frequency close to Ω e . Following the formalism of the mentioned paper, the dispersion equation of the excited wave can be written as a function of adimensional variables: In the present case, the ion beam moves in a decreasing field (inside the magnetic bottle towards the plasmoid), so the term Ω e refers to the particular position where the energy is released: looking again at figure 20, the most part is released in a volume of about 1 cm of thickness, just inside the plasmoid (part in yellow of the picture), where Ω e is between 14.52 GHz and 12.97 GHz. The generated waves let the ECR resonance condition to extend further towards the middle of the plasma chamber, with an effect similar to a magnetic field with a gentle gradient regime (Canobbio [50]). The situation in the background plasma during the instability onset thereby resembles the setup studied by Alton et al [51]. In the mentioned paper, a broadband microwave generator was used for producing 'withe noise signals', in order to broaden the ECR volume. Even in that case the output beam was shown to contain a lower amount of highly charged ions. Another similar situation is obtainable when increasing the / B B min ECR ratio, thus reducing the magnetic gradient at the resonance region too. As a consequence, the plasma starts to emit high energy x-rays and produces lower amounts of highly charged ions [52]. Recently, Tarvainen et al [53] found conditions where the beam ripple and emission of x-rays and microwaves bursts were caused by the variation of the magnetic field profile, making weaker the gradient. A plausible hypotesis is therefore that the former beam-plasma instability artificially broadens the ECR volume, and it makes the plasma prone to the development of further turbulence like the one associated to cyclotron instability [54]. High energy electrons formed by broad resonance volumes, in fact, can eventually lead to cyclotron instabilities, from which the plasma decays emitting cyclotron radiation and plasma electrons: to compensate this extra electron flux, the plasma could expel ions as well, leading finally to the effect observed in figure 19. It is necessary to point out that the proposed hypothesis consider the beam as a 'spark' to ignite the 'fire' of the instability; it then burns thanks to the plasma that makes the fuel.
By evaluating the above mentioned quantities for the choosen simulation ( , where the variation of B is like a 6th order polynomial in Z, see section 2.3). As a consequence, the instabilities excited at the end of the path are much stronger than the ones excited at the entrance of the plasmoid. A dedicated experiment is already under the design phase in order to detect plasma emitted radiation in microwave and x-ray domains, as fingerprint of the excited instability.

Summary and perspectives
A new numerical description of the interaction of an ion beam with a magnetized plasma, in an ECR-based charge breeding device, has been proposed and benchmarked by experimental results: its theoretical investigation revealed the leading role of the many small angle Coulomb collisions in the slowing down and capture phase. For those collisions, three characteristic times can be deduced, expressed by equations (9)÷ (11), that together give an idea for the time it takes thermalization to be reached.
The correct implementation of the above theory on a particle tracking code necessitated the implementation of the Langevin formalism. The physics of the simulated model has been progressively improved by adding more and more sophisticated hypothesis: (i) isotropic/homogeneous plasma with a given temperature and density; (ii) effect of the B field, (iii) implementation of the plasmoid/halo plasma model and (iv) implementation of the potential-dip structure and ionization as mechanisms for ions formation and confinement. The results obtained with the presented approach agree with the theoretical predictions on ion confinement and dynamics, and were benchmarked by experimental results obtained on the LPSC test bench. The simulations evidenced the critical input parameters, listed in the following in order of importance: (i) KT i , (ii) n core , (iii) E inj , (iv) injected beam energy spread and (v) KT e warm .
When exploring the multi-dimensional space of parameters, the simulations revealed to be a powerful tool for explaining the physics of the capture, thermalization and multi-ioniz ation of the charge bred ions: in particular, simulations showed that only a precise combination of the above mentioned parameters matches the experimental results described in [29]: in particular, ∼ ⋅ i . The value of the plasma density is also in good agreement with the estimates reported in the same publication. The simulations clarified also the dynamics of the ballistic phase of the capture process, leading to a rough estimation of the plasma potential, in line with other experimental evidences. Finally, the results obtained put the basis for a possible explaination of the influence of the injected beam on the charge breeder plasma, based on instabilities onset accompained by the excitation of electromegnetic waves inside the plasma.
In this perspectives, the activity carried out opens now to very interesting scenarios: from one side, the code is now able to run as a predictive tool for the charge breeding experiments, with a lot of advantages for the future operation within the SPES facility. Moreover, to have identified the ion temperature as a fundamental parameter in determining the charge breeder performances, provides a useful goal for an R&D program aiming at optimizing this plasma parameter. From another side, the strategy followed so far has permitted also the benchmarking of model assumptions, with a lot of consequences on plasma physics studies. Concerning this last point, the obtained results are in agreement with a view of the plasma based on a plasmoid/halo scheme, with a self-generated potential dip whose strength has been evaluated on the basis of the assumption found in [45].
The very good agreement between the experiments and the simulations encourages us to proceed towards a further optimization of the model. It is worth mentioning that several experimental data feature a critical sensitivity of ECRIS (and so charge breeder) plasma to the electromagnetic field established inside the plasma chamber: the two main effects are the 'two frequency heating' and the 'frequency tuning', whose consequences on the charge breeder performances are still controversial. The beneficial effect of the latter on ECRIS performances, and the fruitful application of the former at Argonne National Laboratories on the local ECR-based charge breeder suggest, as a natural evolution of the work presented in this paper, the necessity of the improvement of the plasma model when taking into account the self-consistent interaction with the microwave field. The work concerning preliminary self-consistent results, obtained through a Maxwell-kinetic loop as simulation strategy, and published in [49], will point the way towards obtaining a 'plasma-target model' under different excitation frequencies.