Modelling RF-plasma interaction in ECR ion sources

This paper describes three-dimensional self-consistent numerical simulations of wave propagation in magnetoplasmas of Electron cyclotron resonance ion sources (ECRIS). Numerical results can give useful information on the distribution of the absorbed RF power and/or efficiency of RF heating, especially in the case of alternative schemes such as mode-conversion based heating scenarios. Ray-tracing approximation is allowed only for small wavelength compared to the system scale lengths: as a consequence, full-wave solutions of Maxwell-Vlasov equation must be taken into account in compact and strongly inhomogeneous ECRIS plasmas. This contribution presents a multi-scale temporal domains approach for simultaneously including RF dynamics and plasma kinetics in a “cold-plasma”, and some perspectives for “hot-plasma” implementation. The presented results rely with the attempt to establish a modal-conversion scenario of OXB-type in double frequency heating inside an ECRIS testbench.


Introduction
Three-dimensional self-consistent numerical simulation of full-wave RF fields in inhomogeneous non-uniformly magnetized ECRIS plasma is a big challenge.Simulations are necessary to improve the comprehension of microwave-to-plasma coupling.Numerical results about electron heating, ionization, ion confinement and ion beam formation in an ECRIS can be very useful in order to increase the performances in terms of output beams currents, charge states, beam quality (emittance minimization, beam halos suppression, etc.) .Several full-wave numerical simulation codes have been developed to describe waves in fusion toroidal devices [1][2][3], exploiting symmetry or ray-tracing technique; at the same time, only few methods have been developed to obtain selfconsistent solutions for nonaxisymmetric ECRIS plasma [4][5][6].
In ECRIS devices (having magnetic field of the order of few Tesla, electron density in the order of 10 18 m -3 and frequencies in the range 3-40 GHz) the Wentzel, Kramers and Brillouin (WKB) method is not applicable because the wavelength is comparable with the plasma gradient and dielectric tensor scale lengths.Only fullwave methods are suitable to capture waves' propagation features, strongly affected by multiple reflections at the cavity-plasma chamber.
During the last years, several numerical simulations of the microwave heating process in ECR plasmas has been carried out at INFN-LNS and LNL [7][8][9].They allowed to give a reasonable explanation of the "Frequency Tuning Effect" affecting the ECRIS plasmas as depending on the excited modes into the plasma chamber.
However, a self-consistent calculation of wave fields and plasma distribution requires a closed loop joining three physics models: Wave-solver, Kinetic solver and collisions.Previous approaches to simulate RF fields in ECRIS using the conductivity kernel of 1-D mirror magnetic field are described in [6].
In this paper we present numerical results of a set of 3D simulations based on the quasi-self consistent approach described in [10,11], including RF wave propagation in a "cold-plasma" solved by a finite element (FEM) method, and a single species (electrons) kinetic code.Our code-suite conjugates FEM solutions of Maxwell equations including the resonator (metallic plasma chamber) effects and the "cold" dielectric tensor of the plasma, with a "particle mover", i.e. a kinetic code developed in MATLAB for solving the single particle relativistic equation of motion of thousands of electrons according to the Boris method [12].
A criticality in microwave-heated magnetoplasmas occurs near the resonance regions, where the plasma becomes a spatial dispersive medium, and the hot magnetized plasmas response has to be considered.In this case, the dielectric tensor depends on the wave vector and the wave solution needs to be calculated by means of iterations.In our FEM code it is proposed a method to take into account the hot plasma response in this paper.

Coupled Electromagnetic and Particle motion model
Assuming a quasi-stationary equilibrium, particles can be followed until they "die" (i.e., they impinge on the chamber walls), while their trajectories are stored in a 3D mesh which finally produces a charge density distribution: this distribution is then included in the following calculation step to compute the electromagnetic field distribution.In a magnetized plasma modelled as a cold magnetofluid with collisions, the field-plasma interaction is described by the tensorial relation εE.Hereby we will apply the "cold plasma" approximation, (i.e.v � >> v th , being v � the wave's phase speed and v th the electron thermal speed).In the case of an ECRIS, due to the shape of its magneto-static structure lacking of uniformity and axial symmetry, ε depends in a complex way from the magneto-static field and the local electron density n e (x,y,z).Full-wave simulations were carried out through the COMSOL Multiphysics FEM solver, considering the geometry shown in Fig. 1 with an initial electron density n e (x,y,z) and the above described "cold plasma" model.The simulations reproduce the main features of the Flexible Plasma Trap that recently came into operations at INFN-LNS [13].The Flexible Plasma Trap (FPT allows axial and radial launching of the waves and it is based on a simple-mirror configuration.The magnetic field is however flexible, allowing "flat" and "beach" magnetic configurations.Since in proximity of the resonance surface (individuated from the iso-surface B 0 =B ECR =m e /eω) the permittivity varies strongly, the discretization of such a narrow region needed a very fine mesh: to achieve this, we adopted an adaptive mesh procedure (that was allowed by a specific feature of the solver).The obtained solution for the electromagnetic fields {E(r), H(r)} was then coupled to the 3D kinetic Matlab® code, initially developed for particles-in-cell calculations of ion dynamics in a plasma.Such a code describes the motion of N macro-particles (being N much smaller than real plasma particles), calculating 3D spatial coordinates and velocities described by the functions r i (t) and v i (t) for i th particle.The calculations take into account only electrons motion, including the so called electro-static (Spitzer) 90° collisions and the Lorentz force: The kinetic code solves equation (2) for each simulated particle applying the relativistic Boris scheme [10,14]: the calculation stops as soon as all the particles hit the chamber walls, while at each time step particles positions are stored in a matrix forming a 3D density map.This map is the output of the kinetic code, and, once rescaled to real values, is used as local electron density n e (x,y,z) to compute again the electromagnetic field {E(r), H(r)}.This "self-consistent loop" runs iteratively, until a convergence is reached.
Hereby we show results up to step 1, i.e. after the evaluation of the magnetized plasma action on RF and viceversa.

Full-Wave simulation strategy
The following calculation uses FPT as simulation domain in order to explore the possibility to obtain overdense plasmas, generating electron Bernstein waves (EBW) by mode-conversion, from the incident O-mode (via OXB conversion) or directly through an X-mode.
The presence of the resonant cavity causes a polarization scrambling, thus strongly complicating the OXB modeling: this makes unpractical some straightforward solutions, already adopted for fusion plasma applications [15].We hereinafter describe our approach: 1) Preparatory virtual experiment: we assume n e and B profiles, then verifying the consistency of the assumed scenarios with modal conversion constraints; 2) Numerical simulations: they allows to check step-1, allowing to find cut-offs and resonances displacements into the plasma; the full-wave FEM code has been applied to study mode conversion mechanism features of O-Xconversion in FPT; 3) Starting from point 2) a microwave launcher has been designed, suitable for modalconversion in FPT.FPT has a cylindrical plasma chamber 260 mm long and 41 mm in radius, fed by microwaves in the range 4-7 GHz generated by a TWT, while the magnetic field is obtained by means of three coils.Fig. 2 and 3 display the electron density n e and magnetostatic field B 0 , used as input parameters for the simulations.In particular, we considered two magnetic field profiles (Fig. 2) as input parameters to compute the "cold" permittivity tensor for 3D full-wave simulations: 1. Simple mirror configuration with double ECR: Single-ECR configuration.
2. "Simple Mirror-single ECR" configuration.The virtual experiment consists in assuming the desired density profile, then selfconsistently checking that the resonance/cut-off layers displacements are compatible with the assumed profile.If the desired profile is a single peaked gaussian with maximum in the central part of the plasma chamber (because of energy transport issues, center-peaked profile is greatly preferable in order to have a hot plasma core), the compatibility check consists in verifying that the energy deposition, occurring at the ECR harmonics, is placed in proximity of the assumed density peak.The cut-off layers where the conversion may happen, as well as the resonance layers for absorption of the electromagnetic waves, can be predicted thanks to our full-wave code as it can be seen in Fig. 4 for the "Simple Mirror-single ECR" configuration: in this figure the chacteristic cutoffs and resonances frequencies  ℎ are plotted along the longitudinal z-axis considering the electron cyclotron frequency   , the plasma frequency   , the Upper Hybrid Resonance   and the L and R-cutoffs   and   .It could be controlled by proper adjustment of the magnetic induction inside the plasma during the experimental phase.The self-consistency check has been performed by means of these type of plots (like Fig. 4).It seems the reciprocal displacement of cutoffs and resonances is suitable for RF power conversion through OXB mechanism.

Numerical Results and Discussion
The following step consisted in the full-wave simulation of the field structure into the plasma chamber.Then, the full system of Maxwell's equations, including the cold magnetized plasma already described in section 2, has been solved for the geometry of the FPT plasma chamber.
Fig. 4 Cutoffs and resonances in FPT linear plasma along 1D direction, assuming the density of Fig. 3 and the case 2 magnetic profile for frequency normalized to heating wave angular frequency 4.65 GHz.
We performed a mesh refinement study, exploiting the adaptive mesh tool based on an error estimate during the solution stage.An extremely fine mesh on the ECR and UHR surface has been obtained, thanks to "functional evaluation" based on the electric field gradient.
Figures 5 and 6 illustrate the COMSOL outputs for electric field distribution in the two magnetic fields cases: for the Simple mirror configuration the electric field has a strong gradient near the ECR layer, that corresponds to resonance-induced wavelength shortening.The electric field strength inside the cavity appears enhanced by the presence of the plasma and, as expected, the largest fraction of input power is absorbed at ECR.The situation in the case 2 (single-ECR point) is somehow different since most of the power seems to be absorbed in the peripheral region of the plasma.The two different displacements of power deposition regions are consistent with visual inspection of the plasma, as it appears from images of figure 7. Experimentally measured X-ray emission is instead a factor 10 larger in case 2, probably due to the smaller volume of the hot plasma, which causes some generation of suprathermal electrons.The fields obtained for both magnetic configurations were used to calculate the electrons motion as described in the following section.

Electromagnetic field applied to the particles
The Particle Mover uses the calculated full-wave RF field in order to compute the phase-space stationary electron distribution by integrating equation of motion according to the Boris method.The simulation time was set to 40 µs, while the time step of each iteration was fixed to 10 -12 s.By dividing the simulation domain in small cells of 1 mm 3 , the code provides two important outputs: the occupation number of each cell (resembling the spatial plasma electron density n e ) and the average energy deposited in each cell, from which a local plasma temperature can be derived.Moreover, the code is able to store positions and velocities of lost electrons for postprocessing.Figure 8 shows the energy distribution of both lost (a) and survived (b) electrons for the case 1. Lost electrons reach energies in the order of tens of keV.Experimentally observed strong X-ray emission (SDD and HpGe X-ray detectors collected photons up to several tens of keVs) is evidently due to such lost energetic particles.Figure 9 shows a comparison between the two magnetic configurations.In particular, figures 9-a and 9-b show the transversal electron density distribution: case 1 is characterized by higher values if compared to case 2, with a more localized density distribution in the center of the trap, in good agreement with the higher light emission rate observed experimentally as shown in figure 7. Figure 9-c and 9-d show the longitudinal (xz plane) energy dnsity distribution: the different displacement of ECR layers (i.e., of the regions of power deposition) is clearly visible; the effect of the magnetic field on the volumetric density distribution is finally shown in figures 9-e and 9f.

"Hot-plasma" model
The next step would be the coupling of a second frequency in presence of a self-consistently established plasma.The idea is to sustain the peaked-profile plasma thanks to the axial microwave injection at frequency f 1 ; then using the radial microwave port in order to inject an O-mode wave at the second harmonic frequency f 2 =2f 1 .A dedicate launcherwith orientable emission lobe [16,17] has been developed on purpose in order to establish proper OXB scenarios.However, in terms of modelling, we need to go beyond cold-plasma approximation.A method to solve the wave-plasma interaction in the "hot plasma" scenario consists in postprocessing to the full-wave fields by means of a Continuous-Wavelet-Transform (CWT).We opted for a complex-Morlet CWT, decomposing the spectral local components along a (magnetostatic) field line.We have started the modeling of well-known systems, possibly giving the opportunity of analytical solutions.We formerly performed the 1D-Wavelet Spectrogram of the E z (z) component on the electric field E, that is shown in figure 10.It is evidentdespite the fact the method must be still improvedthe analysis was able to "capture" the displacement of the two resonances along the chamber axis.Due to the afore-mentioned lack of symmetries in ECRIS machines, 1D-CWT must be extended to a fully-3D system.This upgrade requires many computational efforts which are still ongoing.

Conclusions
Power deposition and plasma density profiles in the FPT scenario have been investigated in order to establish a proper OXB mode conversion scenario.The simulations have been performed with a numerical suite by coupling the full-wave solutions of the field together with particle kinetics.The simulation will support the doublefrequency launching inside the FPT based on a phasevariable waveguides array.

Fig. 2 Fig. 3
Fig.2Profile of magnetic field B 0 along the longitudinal z-

Fig. 5
Fig. 5 Simple mirror magnetic configuration 1 of Simulated electric field |E| [V/m] (on the Top) distribution and dissipated Power P d [W/m^-3] (on the bottom) in proper colours scales range that emphasize the strong gradient near the resonance layers.

Fig. 6
Fig. 6 Simulated electric field amplitude |E| [V/m] distribution and dissipated Power P d [W/m^-3] for the Simple Mirrorsingle ECR configuration in proper scale to emphasize the strong gradient near the resonance layer.

Fig. 7 .
Fig. 7. Imaging of the plasma of FPT in the optical range.Case 1 is depicted on the left: most of the plasma is visible in the central part of the trap.Case 2 is shown on the right: the emission is much weaker and the plasma assumes an annular shape.The rectangular shadow is the axial waveguide.