A model of radiofrequency and plasma coupling for compact ion sources and design

Inductive coupling of radiofrequency power to plasma is a complicate process, since it depends from the density of plasma itself, because ionization is a chain reaction process, and, at low density a capacitive coupling may mix with inductive coupling (with no Faraday screen). Plasma temperature Te , density ne and vector potential are closely coupled, giving nonlinear and singular systems of Partial Differential equations, which require slow iterative solutions, motivating the consideration of a 2D model, also as rapid design and first approximation tool. Plasma conductivity and heating depend on collision rate, which includes also the so-called stochastic collisions (mainly electron collisions with walls), proportionally more important at low gas density ng . Conductivity is also affected by static and radiofrequency magnetic fields; results for skin depth and stochastic collision estimates are reported. The transport of Te and ne inside source can be controlled by a magnetic filter B f . Considering a 5 cm radius H- ion source as example, solution of ne and Te are reported as a function of filter strength B fa, applied rf power and wall status; solver convergence methods and key plasma observables are briefly discussed. Due to small dimension, a filter strength B fa in the order of 8 mT is needed to achieve electron temperature lower than 2 eV (for negative ion production) at extraction.


Introduction
The advantage of contactless heating of plasmas (ionized gases) with radiofrequency (MHz) or microwaves (GHz) [1][2][3][4] well motivates a quest for both a theoretical understanding and a reasonably fast numerical simulation technique, capable to reproduce main experimental behaviour and to help design and optimization of plasma ion sources, also for negative ions and fusion applications [5,6].
Important interrelated physical aspects to be at least approximately accounted for include: coupling between radiofrequency (rf) and plasma electrons [7], ionization and plasma generation, electron collisions [8], transport of plasma ions to walls, effects of confinement magnetic fields and of source walls.
The steady transfer of power from rf (with frequency  = /2) to electrons (called plasma heating) request collisions, so that electrons can dissipate the received energy; plasma heating rate is proportional to electron density   itself.In turn,   depends on ionization due to energetic electrons and on plasma transport to walls.A cylindrical geometry (with lateral wall covered by multipoles [9]) is effective in confining electrons; for rapid calculation a stationary 2D axisymmetric model is here developed.Effect of filter field   and of end wall multipoles is also included.
The rest of the paper is organized as follows.Results for wave equation and plasma rf conductivity  which determines rf penetration in plasma are recalled first [10], including stochastic collisions.Then balance equation for   and   (electron temperature) are added; a constraint on rf adsorbed power helps to stabilize the iterative solution of this nonlinear Partial Differential Equation (PDE) system.Profiles of   and   as function of   and boundary conditions are reported in final section, with results for the efficiency, and extracted currents.

Wave equation, collisions and conductivity
A simplified geometry is shown in figure 1, with model coordinates  and , where as usual , ,  are the cylindrical coordinates and  the Cartesian one.As well known, plasmas typically develop thin sheath layers close to the walls [1], where strong charge unbalance may exist.As discussed later, -1 - we exclude these layers from our simulation domains, which so reduces to a quasi-neutral plasma, enclosed by the sheath input borders (close to the actual walls).Radiofrequency (rf) is here applied to a helical coil.Let   be the plasma chamber radius, ℓ  the coil length, and   the skin depth [1,11], comparable to /  with   ∝  1/2  the angular plasma frequency.In a 2D model, the quantities that depend from  are suitably averaged, for example with B  the magnetostatic induction.Two coupling regimes may be observed [1][2][3]: capacitive coupling (CCP) for low density plasma where   ≥   , so that an rf axial electric field  rf  (  /ℓ  )e i , due to the voltage   applied to coil ends, penetrates into the plasma; inductive coupled plasma (ICP), for dense plasma   ≪   , where the induction field    rf e i enters inside plasma and  rf  0. ICP may be also enforced with a Faraday shield [5,6].For any position x, let θ be a unit length vector in the direction of growing .In ICP, rf is thus given by vector potential A ℜ θ  (, )e i and electric potential  0 (inside plasma) where the operator ℜ (real part of) and time dependence will be understood in the following [8,10].From Maxwell's equations, where 2  is the voltage acting on the -th turn of the coil, and  depends on material features: conductivity  and relative permittivity   ; notation  , indicates the partial derivative of any quantity  respect to any variable .In detail  = − −2 − i 0  +   (/) 2 .
Let Ω  =   / be the static cyclotron frequency.Amplitude  rf of the rf magnetic induction is typically large, so that its cyclotron frequency Ω  =  rf / is much larger than .
With the kinetic energy E  =   |v| 2 /2 and the electron velocity v, let   = (2/3)⟨E  ⟩ be the electron temperature (in energy units), where ⟨⟩ indicate statistical average.An electron thermal velocity definition can be then   ℎ = (8  /  ) 1/2 , in analogy with the case of a Maxwellian distribution.As for ions, they are far less accelerated or decelerated by rf, or by effects of ionization collisions, so their temperature   ≪   is negligible in our simulation domain.
Electron collisions may be of several kinds: ionizations of H 2 molecules, which produce ions and more electrons; elastic collisions or excitations of gas molecules, which contribute to energy loss of electrons; stochastic collisions, where electrons are reflected from walls, or deviated by strong field gradients inside plasma.Momentum collision frequency (total) is estimated as [1,7]: -2 - with   due to stochastic collisions,   due to electron-ion Coulomb collision with ln Λ the Coulomb logarithm, and   due to electron-gas collisions (including elastic and inelastic collisions).The several regimes for   can be summarized by a single fit [12], as a function only of the product     and known parameters.An equivalent approximation   2 , as a function of     and , suitable to rapid numerical calculations, was later proposed [11] and is here compared to   1 computed from the model solution in figure 2; note parameter  = 1 2    2  2  /  as in ref. [8].For each position x, the current density j may depend from E in the plasma nearby, so that conductivity  = δj/ δE is an integral operator, where δ indicates the functional derivative.Magnetic fields (and collisions) limit the electron motions, so in the typical case  ≪   < Ω  with Ω 2  = Ω 2  + 1 2 Ω 2  , we have the local approximation   [8,10], given by note that for Ω  = 0,   is equivalent to  0 , which is the usual model for not magnetized plasma; in other words   is a generalization of  0 to magnetized plasma [10].Few critical remarks [13] and justifications for the use of local conductivity  0 should be mentioned: as shown in ref. [8], a stochastic collision term included in   and in a local  0 conductivity may be derived so that its rf power absorption equates the rf power absorption, duly computed in a kinetic model by integration on the whole electron distribution.This includes also electron reflection at a planar wall sheath, so  0 gives an effective (and simple) approximated representation of kinetic, non-local and stochastic effects.For a wall with curvature radius   this effect strictly applies at distances from the wall much smaller than   , as represented by the /  factor here introduced in eq.(2.2).As regards to   , the satisfaction of its condition   / ≫ 1 assumed in ref. [10] is well verified by figure 2 near the rf driver, where most rf absorption happens.In summary, stochastic collision terms [1,8,11] and eqs.(2.2) and (2.3) can fairly model thermal and kinetic effects.

Balance equations
Let   be the density of positive ions, namely  + ,  + 2 and  + 3 .We model the quasi-neutral plasma only, so that   +  −  =   , with boundaries at the sheath input border [1] and adequate boundary conditions -3 -(bc); CPU-time consuming meshing and solution of plasma sheath layers is so avoided.Moreover, gas ionization is effective only the regions of hotter plasma   ≫ 2 eV, where negative ion density  −  0 and     ; extraction region with colder plasma   < 2 eV is separately discussed.Assuming that gas density   is uniform, rate density of molecule ionization is       (  ), with   (  ) a known function [1,14] (weighted on relevant states), so that electron and positive ion balance gives with   the electron flow and   the positive ion flow.For this coupling, the ion diffusion is helped by the electron temperature   and the ponderomotive rf pressure on electrons, achieving an ambipolar equilibrium [1], as detailed in the following equation and its coefficients.At this ambipolar equilibrium, the ion flow is where   = 1 includes full ponderomotive force [10,15],   (  /  ) 1/2 is the Bohm speed,   an ambipolar diffusion coefficient and g  a gradient; for small gradients, a linear estimate is  1 =   /    , where   and   are the mass and the collision frequency for ions.We define the ambipolar flow velocity v  =   /  ; at domain borders, that is at sheath input,   accelerates up to about the Bohm speed (without surpassing it [10]); so inside domain,   <   (also for conservation of energy).This is obtained by a smooth decrease of   for large gradients [10], as written above.
Eqs. (3.1) and (3.2) give a second order Partial Differential equation (PDE) for   ; a PDE for   is the energy balance where E  is the energy spent by electrons per each ionization achieved [1] and the convective heat flux is neglected [10].The plasma thermal conductivity   is approximated as Boundary condition (bc) for   at sheath border is where n is the outward normal [16],   is the energy gained by an ion in the sheath and     is the ion flow density; eq.(3.5) says that this energy (which is lost on walls) is supplied by the energy flow.
As for   , we consider that wall can re-emit (as e − or perhaps H − ) a fraction   of the impinging particle pairs (e − +H + ), whose flow density is     ; so   boundary condition is On extraction holes,   = 0. Since   is defined also in open space, its boundary condition is   = 0 at large ,  and for  = 0. Few global equations [16] are added to our model: voltages   are solved for, so that a fixed rf current  rf circulate in each coil turn .Among output values, let   be the rf power dissipated -4 -into the plasma,   the power loss in the coil, and   in the wall, so   =   +   +   is the total active power; ideal efficiency   =   /  is easily calculated.A small part of the experimentally applied rf power   is also dissipated in the circuit (cables and matching box, not simulated yet), so say   0.9  .Ion collision frequency may be estimated by   = [(  /  ) 2  +  2  ] 1/2 /  , when the mean free path   is approximately constant [2,10].

Numerical simulations, methods and results
Solution code can be more rapidly written in a multiphysics environment [16]; as a first option [11], consider  rf an input parameter.Due to nonlinearities, system of PDE eqs.(2.1), (3.1) and (3.3) should be solved iteratively, which is slow; but this framework has also advantages: first, eq.(2.1) is complex, so it is solved separately from eqs. (3.1) and ( 3.3) at each iteration.
Second, since eq.( 3.3) is singular for   → 0 and   → 0, our code must verify (or enforce) the conditions   > 0 and   > 0, which is conveniently done at each iteration.
Third, the power setting   = 0.9  is more convenient input parameter then  rf ; during the mentioned iterations,  rf can be adjusted so that total power   progressively approaches   ; more than 20 iterations are then needed for converging to an equilibrium     = 0.9  , allowing direct comparison to experiments.Note that ionization is a chain reaction, where one fast electron gives two electrons as reaction products; so when  rf was used as input parameter, some sharp transition of   and other simulation results were observed for large  rf [11].Using   as input parameter helps to stabilize   and to limit average plasma density, thanks to eq. ( 3.3), which corroborates the use of this input parameter.

Magnetic field and other input data
Our example follows NIO1 [9,17] set-up with plasma chamber   = 48.5 mm,  1 = −78.7 mm,  2 =  1 +   = 133.3mm and a 7 turn coil, centered on the  = 0 plane; frequency  2 MHz.Simulated plasma extends up to  <   with a small gap from   (to remember of sheath layer); for visualization   = 48 mm.A current   = ±400 A flowing in plasma grid PG [17] allows to tune the magnetic filter   .
On the -axis, see figure 3(a), the dipole strength   includes this adjustable filter and a fixed and larger component, due to magnets in the backplate  <  1 and in the EG electrode (after the PG, so  >  2 ).A gaussian filter shape is used, that is with   = 80 mm the filter center and   = 30 mm a parameter.Azimuthal averaging of At maximum current in NIO1, we get   11 mT, so that  fa 8 mT.The   data also include strong multipoles, effective only near wall  =   .
NIO1 rf window spans the  ∈ [− 3 ,  3 ] range;  3 = 38 mm.Optionally, the sensitivity to a wall recycling coefficient   is considered, as given by the trial function with dimensionless parameters  0 and  1 , where Θ  (, ) is the smoothed Heaviside function with a  = 3 mm smoothing length.In detail   =  1 +  0 on rf window,   =  0 on metal walls, with -5 - The condition  1 0 will apply to future Faraday shields.

Key observables for solution post-processing
The   and   solution gives necessary information to calculate experimental observables.First, the light from a given LOS (line of sight) parallel to  axis is collected by a telescope connected to a photomultiplier PMT, which gives a signal   .The power density of Balmer line H  emission averaged on that LOS is with   (  ) the number of ionization events per photon of the Balmer line H  (see eq. 15 and figure 2 of ref. [18]) and   the photon energy; ⟨⟩ LOS indicates average on the LOS.Even if calibration of   against  pmt is missing (and depends on the used telescope or PMT upgrades or PMT gain), the relative changes of  pmt at fixed gain reflect the   changes.Second, the extracted current   − depends on   and   near extraction (where H − are found), and dedicated modeling is well in progress.For a rapid trend estimate as   − =   ⟨    − ⟩  , we infer   − from a comparison with global models [19], in similar volume production conditions; here ⟨⟩  indicates integration on extraction region, say  > 0.9 2 .In detail for   = 1 eV we have  −   1 (  ), with the fit log 10  1 = 1.6(log 10   ) − 10.7 for the range 15.7 < log 10   < 16.9 and   in SI units.Over a threshold as   > 2 eV we consider  −  0. In the range 0.5 <   < 2 eV, a model with  1 = 0.754 eV is used (with  −  constant for   < 0.5 eV).and   = 1300 W case is also plotted.Increasing the filter strength decreases   at extraction, so that  −  there increases, thanks to eq. ( 4.4).To reach a goal   < 1.5 eV for  > 0.9 2 (extraction region),  fa 8 mT seems adequate for both   cases: this proves the usefulness of a filter field, consistently with the NIO1 filter tuning.Current  −  changes from 7 to 6.5 mA with the  fa span, while coextracted electron current estimate decreases by about 40% and  −  estimate (in extraction region) improves by about 40%. Figure 4(b) shows that by increasing   , that is, by moving filter away from the driver, the   < 2 eV region shrinks as expected, but the   < 1.5 eV stay constant; so current slightly improves: in detail   − = 6.0 mA and   = 0.73 for   = 7 cm, and   − = 6.8 mA and   = 0.75 for   = 9 cm.

Results for 𝒔 𝒆𝒆 > 0
The cases  0 ,  1 > 0 (figures 5 and 6) show more current, since the lower loss of charged particles on walls allows for a lower equilibrium   , also at extraction; moreover, when  0 is increased, the   increases for all , so more electrons can arrive at extraction  2 and   ( 2 ) can grow.When  1 is increased (see figure 5) the effect on   remains, while effect on   ( 2 ) is very small.The plasma heating efficiency   increases with  1 ; for   = 1500 W and  1 as in figure 5, we respectively have   = 0.741, and 0.746, and 0.751 and 0.757 and 0.765; current increases from 6.8 to 7.8 mA. Figure 6 shows that while  −  increases with  1 , global radiated   decreases (predicting the experimentally observed anticorrelation of   and  −  [17] at any fixed   value, in this figure   = 1.5 kW).When  0 rises from 0 to 0.3, with  1 = 0.1, the current increase is larger, from 6.2 mA to 8.4 mA; anticorrelation between   and  −  is still verified.On the other hand, both   and  −  estimates increase when   increases, at least in the case  0 = 0.In summary, this model can predict significant trends and correlations of many observed quantities.
In perspective, several model upgrades can be considered, for example: modeling the sheaths or the extraction or modeling ion species density, or 3D models, surely deserving further investigations; these models request not only huge computational resources, but also dealing with equation nonlinearity, singularities and many length scales, in a wider space.For that reason, the present 2D model solutions represent a useful reference and an initial guess.Moreover, the solutions for the shown cases of   ,   and  fa predict many experimental observables and quantify the efficacies of those parameters to control   and   profiles.

Figure 1 .
Figure 1.Geometry and objects for radiofrequency and plasma simulations; note plasma grid PG and electrode EG.

4. 3 Figures 2
Figures2 and 4show typical results for   and   , in the case   = 0, for   = 0.75 ± 0.05 Pa.Since rf absorption and ionization happen near the  = 0 plane, that is in the coil region || <   , density strongly peaks at plasma center  =  = 0; temperature peak is near to coil  = 0,  = 4 cm.On -axis,   profiles for several  fa ≤ 10 mT and   = 1500 W are compared; the typical  fa = 8 mT

Figure 6 .
Figure 6.Correlation of estimated plasma light emission   and extracted ion current  −  , for the figure 5 cases with fixed  0 = 0.1 (triangles); and varying  0 as labeled (circles);   = 1500 W and   0.7 Pa.