Modeling space-resolved ion dynamics in ECR plasmas for predicting in-plasma β-decay rates

Lifetimes of radioactive nuclei are known to be affected by the level configurations of their respective atomic shells. Immersing such isotopes in environments composed of energetic charged particles such as stellar plasmas can result in β-decay rates orders of magnitude different from those measured terrestrially. Accurate knowledge of the relation between plasma parameters and nuclear decay rates are essential for reducing uncertainties in present nucleosynthesis models, and this is precisely the aim of the PANDORA experiment. Currently, experimental evidence is available for fully stripped ions in storage rings alone, but the full effect of a charge state distribution (CSD) as exists in plasmas is only modeled theoretically. PANDORA aims to be the first to verify these models by measuring the β-decay rates of select isotopes embedded in electron cyclotron resonance (ECR) plasmas. For this purpose, it is necessary to consider the spatial inhomogeneity and anisotropy of plasma ion properties as well as the non-local thermodynamic equilibrium (NLTE) nature of the system. We present here a 3D ion dynamics model combining a quasi-stationary particle-in-cell (PIC) code to track the motion of macroparticles in a pre-simulated electron cloud while simultaneously using a Monte Carlo (MC) routine to check for relevant reactions describing the ion population kinetics. The simulation scheme is robust, comprehensive, makes few assumptions about the state of the plasma, and can be extended to include more detailed physics. We describe the first results on the 3D variation of CSD of ions both confined and lost from the ECR trap, as obtained from the application of the method to light nuclei. The work culminates in some perspectives and outlooks on code optimization, with a potential to be a powerful tool not only in the application of ECR plasmas but for fundamental studies of the device itself.

1 Introduction

Electron cyclotron resonance ion sources (ECRIS) have long been used as reliable devices supplying ion beams to particle accelerators with charge states and currents tuned as per experimental requirements [1]. Their working principle involves sequential ionization of atoms through collision with electrons which are energized using resonance heating and trapped long enough through magnetic confinement. By varying the frequency of the microwave radiation used for heating, as well as associated power, magnetic field profile and gas pressure, the charge state distribution (CSD) of the ions and output current can be tuned [2].

Coincidentally, the existence of a CSD in the ECR plasma interior is a doorway to performing in-plasma measurement of nuclear β-decay rates to emulate astrophysical scenarios. Decay lifetimes can vary according to the configuration of the surrounding atomic shell [3, 4] and such a modification has already been experimentally observed in fully stripped ions circulating in storage rings [5]. The PANDORA (Plasmas for Astrophysics, Nuclear Decay Observations, and Radiation for Archaeometry) facility aims to take this a step further and analyze the CSD-dependent t_1/2 of radioisotopes [6] based on the theory of perturbed decay rates in stellar environments [7]. For each standard β⁻ decay—continuum decay (cd), bound-state decay (bd), electron capture (ec), and continuum capture (cc)—there exists an associated Q-value which governs the reaction spontaneity.

$Q_0=\left\{\begin{array}{cc}{m}_0\left({}_{\mathrm{Z}}{}^{\mathrm{A}}\mathrm{X}_{\mathrm{N}}\right)c^2-m_0\left({}_{\mathrm{Z}+1/\mathrm{Z}-1}{}^{\mathrm{A}}\mathrm{Y}_{\mathrm{N}-1/\mathrm{N}+1}\right)c^2\hfill & \left(\mathrm{c}\mathrm{d},\mathrm{b}\mathrm{d},\mathrm{e}\mathrm{c}\right)\hfill \\ m_0\left({}_{\mathrm{Z}}{}^{\mathrm{A}}\mathrm{X}_{\mathrm{N}}\right)c^2-m_0\left({}_{\mathrm{Z}+1/\mathrm{Z}-1}{}^{\mathrm{A}}\mathrm{Y}_{\mathrm{N}-1/\mathrm{N}+1}\right)c^2+K_{e_f^{-}}\hfill & \left(\mathrm{c}\mathrm{c}\right)\hfill \end{array}\right..(1)$

Here, m₀(X)c² and m₀(Y)c² represent, respectively, the rest mass energy of the parent and daughter system, while $K_{e_f^{-}}$ is the kinetic energy of free electron captured from the continuum. When the parent atom is ionized or excited, the decay energetics may be modified as

$Q_{i,j}=Q_0+\left(E_{X,K}^{\ast }-E_{Y,K^{\prime }}^{\ast }\right)+\left({ϵ}^{i,j}-{ϵ}^{i^{\prime },j^{\prime }}\right)+\left({\mathrm{\Delta }}_X-{\mathrm{\Delta }}_Y\right),(2)$

where $E_{X,K}^{\ast }$, $E_{Y,K^{\prime }}^{\ast }$ are the energies of the two nuclei in levels K and K′, ϵ^i,j; ϵⁱ′^,j′ are the atomic energies of the parent (daughter) system in the i (i′) charge state and j (j′) atomic/ionic level, and Δ_X, Δ_Y are the contributions of the ionization potential depression from the surrounding charges. In case of decays with small Q₀ (∼keV), the contribution of atomic levels can be strong enough to open up new decay channels or suppress existing ones. For any transition K → K′, the decay half-life t_1/2 can be related to nuclear properties through the expression.

$f_{K\to K^{\prime }\left(m\right)}{t}_{1/2}=\frac{\left(\ln 2\right)2{\pi }^3{\hslash }^7}{g^2{m}_e^5{c}^4|M_{K\to K^{\prime }}^{L\left(m\right)}{|}^2},(3)$

where |M_K→K′| is the nuclear matrix element (NME) describing the coupling between the initial and final nuclear states, f_K→K′(m) is the lepton phase volume, m is the type of decay (allowed, unique forbidden, or non-unique forbidden) related to the angular momentum of the transition L(m), and g is the weak interaction coupling constant. The RHS of Eq. 3 is atom-independent, making f_K→K′ the chief quantity of interest when investigating the change in t_1/2. In the presence of a CSD, and by extension a level population distribution (LPD), the lepton phase volume can be calculated as

$f_{K\to K^{\prime }\left(m\right)}^{\ast }=\left\{\begin{array}{cc}\sum _{i,j}{p}_{i,j}{\int }_1^{W_{\mathit{max}\left(i,j\right)}}{\left(W^2-1\right)}^{1/2}W{\left(W_{\mathit{max}\left(i,j\right)}-W\right)}^2{F}_0{S}_{\left(m\right)}\left(i,j\right)f_d\left(i,j\right)\mathrm{d}W\hfill & \left(\mathrm{c}\mathrm{d}\right)\hfill \\ \sum _{i,j}{p}_{i,j}{\int }_{W_{\mathit{min}\left(i,j\right)}}^{\infty }{\left(W^2-1\right)}^{1/2}W{\left(Q_{i,j}/m_e{c}^2\right)}^2{F}_0{S}_{\left(m\right)}\left(i,j\right)f_c\left(i,j\right)\mathrm{d}W\hfill & \left(\mathrm{c}\mathrm{c}\right)\hfill \\ \sum _{i,j}{p}_{i,j}\sum _{x\left(i,j\right)}{\sigma }_x\frac{\pi }{2}{\left[g_x\mathrm{o}\mathrm{r}{f}_x\right]}^2\left(Q_{i,j}/m_e{c}^2\right)S_{\left(m\right)x\left(i,j\right)}\hfill & \left(\mathrm{b}\mathrm{c},\mathrm{b}\mathrm{d}\right)\hfill \end{array}\right.,(4)$

and the modified decay rate λ given as

$\lambda =\frac{\ln 2}{f_{K\to K^{\prime }\left(m\right)}{t}_{1/2}}\sum \underset{K\to K^{\prime }\left(m\right)}{\overset{\ast }{f}}.(5)$

The summation in Eq. 5 accounts for the contribution of multiple decay channels to the overall rate. In Eq. 4, p_i,j stands for the occupation probability of the charge state i and level j, W_{max  / min(i,j)} are functions of Q_i,j, F₀ is the Fermi function, S_(m)(i, j) and S_(m)x(i,j) are shape factors associated to the momentum transfer, f_c/d(i, j) are related to the relativistic Fermi–Dirac distribution function describing the continuum electrons, x(i, j) represents the electronic levels within the atomic configuration j, σ_x is the probability of electron occupation/vacancy in that level and g_x/f_x is the maximum of the orbital radial wavefunction evaluated at the nuclear surface. It is quite evident then that calculating in-plasma decay rates requires inputs on the LPD and CSD of the plasma ions, as well as on the electron energy distribution function. Acquiring such data is a complicated process because ECR plasmas are inherently anisotropic and non-homogeneous, host multi-component electron distributions, and obey non local thermodynamic equilibrium. This calls for robust simulation tools capable of extracting spatially resolved information on charged particle properties, and associated comprehensive diagnostic methods to experimentally benchmark their results.

Over the past years, several research groups have tackled the aforementioned problem from different angles—a separate section has been devoted to discussing their methodologies and subsequent results. The INFN-LNS and LNL groups have also contributed to ECRIS simulations by developing a self-consistent model coupling electron dynamics with an EM field solver to obtain 3D space-resolved maps of electron density and energy and reproduced some well-known ECRIS phenomena like the frequency tuning effect and heavy ion charge breeding [8, 9]. These steady-state electron maps can now be used as a base to arrive at the complementary, stationary ion maps containing CSD and LPD of the species of interest as a function of their position in the plasma. This article is intended to provide an outline of the 3D coupled ion dynamics and population kinetics algorithm developed to achieve the aforementioned aim.

The contents of the article are divided into eight subsequent sections as follows: in Section 2, an in-depth review of ECRIS simulation tools is presented, followed by Section 3 wherein a brief description of the self-consistent MATLAB^ⓒ-COMSOL Multiphysics^ⓒ simulation scheme used to obtain 3D maps of electron density and energy is provided. This is extended to cover the ion dynamics algorithm in Section 4. The important modules of the algorithm, namely, the particle pusher code, Monte Carlo (MC) sampling scheme, density scaling methodology, and electrostatic field evaluation are covered in detail in Sections 5–7 respectively. The complete code flow is detailed in Section 8 while the results of the simulation for a specific plasma configuration and select operating conditions are shown in Section 9 and compared with experimental data for benchmarking. Finally, important takeaways from the procedure are presented in Section 10, as well as perspectives for future work.

2 ECR simulation models: A review

Research into ECRIS simulation has been underway since the early 1990s, with most groups focusing on the dominant electrons. Heinen and coworkers developed codes to track the motion of electrons under the effect of the ECR magnetostatic field with and without collisions and included resonance heating from single-mode microwaves. They obtained insights into the 3D profiles of electron density and energy and used the results to explain well-known phenomena such as the afterglow effect, and established a connection between localization of electron energies and production of highly charged ions [10, 11]. A similar approach was also followed by Maunoury et al who used the TrapCAD code to investigate the properties of lost and confined electrons and indirectly inferred ion properties [12]. This code was developed by Vámosi and Biri in the 1990s, with the objective to calculate and visualize the magnetostatic field of the ECR trap and follow the paths of charged particles [13]. Analysis of electron population through relativistic particle-in-cell (PIC) simulations was also conducted by Dougar-Jabon, Umnov and Diaz [14] who found specific zones of periodic bounce oscillations and banana trajectories as well as multiple electron populations.

The first algorithms aimed at coupling electrons and ions were developed by Edgell and group, and implemented in progressively advanced versions of the GEM code. They used the Fokker–Planck equation with ECR heating to obtain the electron distribution in phase space and gave it as input to detailed density and power balance fluid equations that described the collisional ions. The codes were initially 0D [15] but later upgraded to model 1D axial transport as well [16] and successfully predicted the extracted ion CSD. The 1D nature of the code was, however, quite limiting, and applicable to only the plasma core from which ion loss is modeled.

Cluggish, Kim, and Zhao later used the 1D GEM code to open investigations into the effect of pressure and microwave power on plasma parameters such as total electron density, energy, and current. They deduced some scaling laws which were in general agreement with experiments performed on different ion sources [17]. They realized the importance of studying the feedback of the plasma on the wave propagation and absorption and developed a quantitative 1D model and qualitative 3D model to describe the RF-plasma coupling based on the hot electron approximation [18]. The idea and implementation were robust in 1D but lacked completeness in 3D owing to the absence of a full dielectric tensor. Additionally, both works were centered on the electrons and ion properties were inferred on a macroscopic scale.

The most rigorous and complete treatment of ECR plasmas in terms of electrons, ions, and microwaves has been performed by Mironov and coworkers starting from 2009. In contrast to the methodologies utilized by other groups, they began by focusing directly on the ion dynamics using a particle-in-cell Monte Carlo collision (PIC-MCC) model, taking electron energy distribution function as isotropic and uniform [19]. They simulated the motion of macroparticles under a min-B field assuming collisions, implemented atomic processes such as ionization and recombination in real time, and took into account the constant influx of neutrals to replenish extracted charges just like in a real ECRIS. As a consequence, they obtained a great degree of match between predicted and measured steady-state CSD and were able to analyze various aspects of ECRIS operation such as neutral density gradient and isotope anomaly effect. In later works, they extended their analysis to include multi-component, localized electron energy distribution functions and investigated the ion confinement processes in the plasma [20] and then applied their algorithms to investigate the behavior of the device with different working materials [21]. The local neutrality assumption was dropped by coupling an electron simulation code NAM-ECRIS(e) with its ion counterpart NAM-ECRIS(i), effectively modeling electron and ion evolution simultaneously in the plasma [22]. The final upgrade came in the form of RF-electron coupling in NAM-ECRIS(e) using the cold electron approximation on the ECR surface [23].

The simulation scheme reported in this work is complementary to the approach of Mironov and group in that it also aims to investigate ECR plasmas as a whole by appropriately extending the electron simulation module to the ions. While their method captures temporal evolution of the system in real time, we instead focus on developing suitable models to represent steady-state conditions. It can be argued that the latter is merely a special case of the former, but it gives us the advantage to model RF-electron coupling with a high degree of precision while still capturing the essence of ion dynamics. The details of this methodology are presented in the following sections.

3 Self-consistent electron simulations under cold electron approximation

The importance of electrons in ECRIS has been well-established through the numerous works reported in Section 2. Mascali et al also realized the necessity to treat the plasma chamber as a resonant cavity to improve understanding of both electron and ion dynamics. They developed a hybrid model simulating electrons using a PIC code and ions using MC collisions, taking into account resonance heating from the standing-wave profile of microwaves in the resonant cavity. They reproduced the experimental beam shape and predicted the formation of the plasmoid-halo structure [24]. Some enhancements to the electron PIC module were later made by Neri et al who implemented Spitzer collisions and studied the effect of RF heating on electron confinement [25].

The simplest of ECRIS devices are operated with a peculiar magnetostatic field profile called the min-B structure as shown in Figure 1 which enables electron confinement along the axis. The electrons naturally gyrate about the magnetic field lines at their cyclotron frequency ω_c = e|B|/m_e where e is the electric charge, |B| is the field magnitude, and m_e is the mass of the electron. Depending on the position, the field may be such that ω_c matches the frequency of the circularly polarized microwave launched into the plasma, leading to resonance heating. On account of magnetic and diffusion transport, the distribution of electrons in the plasma quickly becomes non-homogeneous and anisotropic, directly affecting the dielectric tensor of the medium and therefore the EM field profile. As mentioned in Section 2, the feedback of the plasma to the wave propagation was already calculated in [18] but using the hot electrons which made it necessary to replace the dielectric tensor with a constant. Torrisi et al revisited the subject to perform FEM calculations to deduce the field profile in magnetized plasmas using the cold electron approximation instead [26] leading to a more complete yet still precise formulation of RF-electron coupling since cold electrons constitute a majority of the plasma.

FIGURE 1

FIGURE 1. (A) 3D profile of magnitude of the magnetostatic field superposed with streamlines showing field flow and ECR surface, (B) 1D on-axis field profile with demarcated ECR zones. The min-B structure of the field can be clearly appreciated.

Based on their algorithm and results, Mascali et al developed an electron kinetics model to solve the collisional Vlasov–Maxwell and Fokker–Planck equation using a large number of representative macroparticles [9]. The idea was to tackle the wave–plasma coupling in a self-consistent manner by iteratively calculating the EM field distribution using COMSOL Multiphysics ^ⓒ based on the cold tensor approximation in the plasma, and then pushing the electrons using a particle mover code written in MATLAB ^ⓒ according to the Lorentz force including the newly generated EM field. By constantly accumulating macroparticle traces in cells of a 3D matrix, they obtained relevant occupation maps and used them to update the dielectric tensor, repeating all the steps till self-consistency between the electron maps and the EM field was achieved. In the end, one obtained 3D space-resolved maps of electron density and energy corresponding to the ECR plasma at steady-state, in the form of 3D MATLAB ^ⓒ matrices. Figure 2 shows the density and energy maps of electrons simulated in an Ar plasma operated at frequency 14.28 GHz and power 200 W. These maps will serve as the source of electron-driven reactions included in the ion dynamics model described in the following section. The simulation model is currently being upgraded and the details can be found in [27].

FIGURE 2

FIGURE 2. (A) 3D profile of simulated electron density (in m⁻³) and, (B) electron energy (in eV).

4 Coupled ion dynamics and population kinetics

The standard approach to model ion population kinetics is to start from the density balance equation.

$\frac{\mathrm{d}{n}_i}{\mathrm{d}t}=\sum _{j\ne i}{n}_j{G}_{ji}-\sum _{i\ne j}{n}_i{L}_{ij}-\frac{n_i}{{\tau }_i},(6)$

where n_i,j are level populations of states i and j, G_ji and L_ij represent inter-level gain and loss rates, respectively, and τ_i is the characteristic confinement time. On extending Eq. 6 to include all levels and collision-radiative processes, a rate matrix R can be constructed and Eq. 6 can be reduced to $\frac{\mathrm{d}\mathbf{n}}{\mathrm{d}t}=\mathbf{R}\cdot \mathbf{n}$. Generally, Eq. 6 is further simplified in ECR plasmas by considering only three important processes—electron collision ionization, charge exchange through collision with neutrals, and ion transport [15, 28]. The resultant balance equation is then

$\begin{array}{l}\hfill \frac{\mathrm{d}{n}_i}{\mathrm{d}t}=n_e{n}_{i-1}\langle {\sigma }_{ion,i-1\to i}{v}_e\rangle -n_e{n}_i\langle {\sigma }_{ion,i\to i+1}{v}_e\rangle +\\ \hfill & n_0{n}_{i+1}\langle {\sigma }_{CEX,i+1\to i}{v}_i\rangle -n_0{n}_i\langle {\sigma }_{CEX,i\to i-1}{v}_i\rangle -\frac{n_i}{{\tau }_i},\hfill \end{array}(7)$

where n_e is the electron density, n_i/i−1/i+1 represents the ion density in charge state i/i − 1/i + 1, σ_ion is electron collision ionization cross section, v_e is the electron speed, σ_CEX is charge exchange cross section, and v_i represents the ion collision speed. Although Eq. 7 is complicated enough in 0D, it is quite impossible to solve in 3D wherein the densities and energies of the interacting species vary spatially and particle transport takes place simultaneously with reactions.

We decided to adapt Eq. 7 to the established electron the PIC code mentioned in Sec. 3 to calculate steady-state ion maps [29]. The idea was to generate N macroparticles in the i⁺ state obtained from the (i − 1)⁺ state (first term in Eq. 7) and track their evolution over a fixed time T_span discretized into a number of steps T_step. At each step, ionization to the (i + 1)⁺ state and exchange into the (i − 1)⁺ state would be evaluated according to their respective reaction probabilities (second and fourth terms in Eq. 7), resembling the rate matrix method. Charge exchange to the (i + 1)⁺ state (third term in Eq. 7) is ignored in this model. The remaining unperturbed macroparticles would be moved according to the equation of motion under the action of the Lorentz force F_L,i = eq_i(E_DL + v_i ×B) and the Langevin formula, with q_i being the charge state of the macroparticle and v_i its velocity. The quantity E_DL represents the electrostatic field arising from the double layer in the plasma whose origin and importance are discussed in detail in Section 7. During the course of their motion, the macroparticle tracks would be mapped to a 3D occupation matrix n_i in much the same way as for electrons, until their eventual loss from the plasma (last term in Eq. 7). The ionized/exchanged particle positions would instead be added to corresponding ionization/exchange matrices, n_i→i+1 and n_i→i−1. At the end of the simulation, the occupation, ionization, and exchange maps would be passed to the next charge state as inputs, and all the steps repeated. The results of each preceding simulation would be embedded into the successive one. If the reaction rates are implemented correctly, the confined macroparticle tracks and relative weight of state calculated according to reaction-mediated transfer would ultimately lead to a CSD map consistent with the fixed electron distribution. A summary of the simulation scheme is shown in Figure 3.

FIGURE 3

FIGURE 3. (A) Flowchart portraying the various inputs, outputs and steps in the simulation of a single charge state i⁺. ION, OCC, and CEX refer to ionization, occupation, and charge exchange matrices, respectively. (B) Overall network of simulation schemes showing flow of data between successive charge states. Ar⁰⁺ represents neutral Ar atoms since the simulations were run for such plasmas.

The key elements of this algorithm are the ion transport model, MC population kinetics, density scaling, and double layer field calculation. The subsequent sections contain more comprehensive discussions on these topics. The simulation scheme was configured for the ATOMKI, Debrecen ECR plasma trap for the same operating conditions as the electron maps in Figure 2 (14.25 GHz frequency, 200 W power, and 10^–6 mbar pressure) using Ar plasma and thus all discussions hereafter will take place assuming such a plasma.

5 Particle pusher algorithm: Ion transport

The motion of charged particles in a plasma is a notoriously difficult problem to tackle because unlike neutral atoms in a rarified gas which mostly travel in straight lines before stochastic collisions, plasma charges are subject to both long-range EM forces as well as short-range random collisions. The former, additionally, can arise from both internal as well as external fields, such as the magnetostatic field applied in ECR plasmas. The result of these forces is that ions and electrons in plasma gyrate along magnetic field lines while being simultaneously subject to continuous and infinitesimal perturbations in their trajectory on account of internal EM fields [30] while also getting randomly kicked due to head-on collisions with fellow constituents. The PIC simulations for electrons as described in Section 3 already addressed this complexity, allowing direct adaptation to the ion dynamics scheme. In the presence of electric and magnetic fields, the equations of motion of a charged particle can be analytically written as

$\frac{\mathrm{d}\mathbf{r}}{\mathrm{d}t}=\mathbf{v},(8)$

$m\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}=q\left(\mathbf{E}+\mathbf{v}\times \mathbf{B}\right),(9)$

where r and v are the position and velocity of the particle, respectively, q is the charge, m is the mass, E and B are the electric and magnetic field, respectively, and the RHS of Eq. 9 is the Lorentz force. Numerical implementation of the same has been traditionally performed with the Boris method [31] which solves the particle motion in a leapfrog manner and breaks down the acceleration into a number of steps.

$\frac{{\mathbf{r}}^{n+1}-{\mathbf{r}}^n}{T_{\mathit{step}}}={\mathbf{v}}^{n+1/2},(10)$

${\mathbf{u}}^{-}=\gamma {\mathbf{v}}^{n+1/2}+\frac{q}{2m}{\mathbf{E}}^{n+1}{T}_{\mathit{step}},(11)$

$\theta =\frac{q{T}_{\mathit{step}}}{m\gamma }|{\mathbf{B}}^{n+1}|,(12)$

$\mathbf{t}=\tan \frac{\theta }{2}\mathbf{b},(13)$

${\mathbf{u}}^{\prime }={\mathbf{u}}^{-}+{\mathbf{u}}^{-}\times \mathbf{t},(14)$

${\mathbf{u}}^{+}={\mathbf{u}}^{-}+\frac{1}{2+|\mathbf{t}{|}^2}\left({\mathbf{u}}^{\prime }\times \mathbf{t}\right),(15)$

${\mathbf{v}}^{n+3/2}=\frac{1}{\gamma }\left({\mathbf{u}}^{+}+\frac{q}{2m}{\mathbf{E}}^{n+1}{T}_{\mathit{step}}\right).(16)$

Here, $\gamma ={[1-(|\mathbf{v}{|}^2/c^2)]}^{-1/2}$ denotes the known Lorentz factor, n is the index tracking the leapfrog procedure, b = Bⁿ⁺¹/|Bⁿ⁺¹| is the unit vector of the magnetostatic field, and T_step is the discrete time interval. To save on computational time, most groups apply Eq. 13 in the “small θ” limit such that the tangent is replaced by its argument. Although this is a good approximation, the gyrophase error tends to accumulate if the number of iterations is large, which is precisely the case here. As such, a modified form of the Boris solver was used for the ion transport which is analytically exact and relatively time-saving [32]. According to the formalism, Eqs 13–15 were replaced with a rotation matrix R_rot such that

${\mathbf{u}}^{+}=R_{\mathit{rot}}{\mathbf{u}}^{-},(17)$

and R_rot is defined as

$R_{\mathit{rot}}=\left(\begin{array}{ccc}\hfill \cos \theta +b_x^2\left(1-\cos \theta \right)\hfill & \hfill b_x{b}_y\left(1-\cos \theta \right)+b_z\sin \theta \hfill & \hfill b_z{b}_x\left(1-\cos \theta \right)-b_y\sin \theta \hfill \\ \hfill b_x{b}_y\left(1-\cos \theta \right)-b_z\sin \theta \hfill & \hfill \cos \theta +b_y^2\left(1-\cos \theta \right)\hfill & \hfill b_y{b}_z\left(1-\cos \theta \right)+b_x\sin \theta \hfill \\ \hfill b_x{b}_z\left(1-\cos \theta \right)+b_y\sin \theta \hfill & \hfill b_y{b}_z\left(1-\cos \theta \right)-b_x\sin \theta \hfill & \hfill \cos \theta +b_z^2\left(1-\cos \theta \right)\hfill \end{array}\right),(18)$

The evolution of phase space of a test particle in a plasma due to long- and short-range collisions is usually modeled by the Fokker–Planck equation and simplified using the formalism of McDonald and Rosenbluth [33]. Further reduction is made possible by assuming the species of interest as test particles moving through a field composed of other plasma constituents. The previously mentioned procedure yields two important quantities at the end—friction frequency ν_s and diffusion coefficients D_‖, D_⊥ which, respectively, correspond to the continuous drag force motion and the head-on collisions. The PIC electron simulations adapted the Fokker–Planck equation to macroparticle formalism by using ν_s, D_‖, and D_⊥, in the Langevin equation which is identical to the former to first order accuracy in T_step. The theory has been rigorously discussed in [34]. Taking, at any point in time, a surrounding total plasma density n_s with particle mass m_s, atomic number Z′ , and most probable thermal speed as c_s, the friction and diffusion coefficients are obtained as

$D_{\Vert }=\langle {\mathbf{v}}_{\Vert }^2\rangle =\frac{A_D}{|{\mathbf{v}}^{n+1/2}|}G\left(\frac{|{\mathbf{v}}^{n+1/2}|}{c_s}\right),(19)$

$D_{\perp }=\langle {\mathbf{v}}_{\perp }^2\rangle =\frac{A_D}{|{\mathbf{v}}^{n+1/2}|}\left\{\mathrm{\Phi }\left(\frac{|{\mathbf{v}}^{n+1/2}|}{c_s}\right)-G\left(\frac{|{\mathbf{v}}^{n+1/2}|}{c_s}\right)\right\},(20)$

${\nu }_s=\left(1+\frac{m}{m_s}\right)\frac{A_D}{c_s^3}\frac{G\left(\frac{|{\mathbf{v}}^{n+1/2}|}{c_s}\right)}{\frac{|{\mathbf{v}}^{n+1/2}|}{c_s}},(21)$

where

$A_D=\frac{{\left(Z{Z}^{\prime }\right)}^2{e}^4{n}_s\ln \mathrm{\Lambda }}{2\pi {ϵ}_0^2{m}^2},(22)$

Φ(x) is the error function and G(x) = [Φ(x) − Φ′(x)]/2x². Z represents the atomic number of the test particle while ln λ is the Coulomb logarithm. Using these coefficients, contributions to velocity of the test particles were calculated as

${\mathbf{v}}_{\mathit{fric}}={\nu }_s{\mathbf{v}}^{n+1/2}{T}_{\mathit{step}},(23)$

${\mathbf{v}}_{\mathit{rand}}=P_{\Vert }N\left(0,D_{\Vert }\right)+P_{\perp }N\left(0,D_{\perp }\right),(24)$

where N(0, D_‖), N(0, D_⊥) indicate random sampling from a normal distribution with mean 0 and standard deviation as given by the diffusion coefficients, while P_‖, P_⊥ are projection matrices mapping the parallel and perpendicular components to the basis vector of the velocity. The final expression to advance the velocity was then

${\mathbf{v}}^{n+3/2}=\frac{1}{\gamma }\left({\mathbf{u}}^{+}+\frac{q}{2m}{\mathbf{E}}^{n+1}{T}_{\mathit{step}}\right)+{\mathbf{v}}_{\mathit{fric}}+{\mathbf{v}}_{\mathit{rand}}.(25)$

The motion of ions in the plasma was fully realized through Eqs. 10–12, 17, and 25.

6 Generalized Monte Carlo routine: Population kinetics

The density and temperatures of ECR plasmas are such that the various constituents of the system obey non local thermodynamic equilibrium, meaning ion LPD and CSD cannot be described analytically using the Boltzmann and Saha equations, respectively. Instead, one needs to apply the general collision–radiative model using the rate matrix formalism as described in Eq. 6, and take into account a multitude of reactions and levels to construct the differential equation. The model accuracy is directly correlated with the amount of physics considered because of the following reasons.

1) The reaction rates govern the percentage transfer of macroparticles within charge states and ionic levels thus defining the relative weights of each charge in the overall CSD

2) The ion transport and occupation map rely on the self-generated plasma density n_s (Section 5) and consequently, on the weights of the various charge states.

As such, while one may be limited by the availability of data on reaction cross sections or atomic levels, it is nevertheless a good idea to develop a generalized algorithm to model the population kinetics which can be later expanded to include as many levels and reactions as needed. This is particularly important not just for modeling the CSD and LPD in the plasma, but also for the β-decay model in [7].

In keeping with the simplified balance equation in Eq. 7, electron collision ionization and charge exchange after collision with neutral Ar atoms were the only reactions considered. The analytical formula for single ionization i → i + 1 is given as

${\sigma }_{ion,i\to i+1}=\frac{10^{-17}}{I_i{E}_e}\left[\sum _{n=1}^6{A}_n{\left(1-\frac{I_i}{E_e}\right)}^n+B\ln \left(\frac{E_e}{I_i}\right)\right],(26)$

where σ_ion,i→i+1 is the ionization cross section in m², I_i and E_e are the ionization energies of the charge state indexed by i and electron collision energy (in eV), respectively, and A_n, B are the fitting coefficients which vary with i. The expression, as well as the values of the coefficients, is detailed in [35].

The charge exchange cross section for the transition i → i − 1 was taken from the works of Müller and Salzborn [36] and given as

${\sigma }_{CEX,i\to i-1}={Ai}^{\alpha }{\left(I_0\right)}^{\beta },(27)$

where σ_CEX is the exchange cross section in m²; A, α, and β are fit coefficients; and I₀ is the ionization energy of the 0⁺ → 1⁺ transition of Ar which acts as the “target” in this reaction. The values of the fit coefficients are taken as A = 1.43 × 10^–17, α = 1.17, and β = 2.76.

In order to reconcile the rate matrix formalism with macroparticles in a PIC algorithm, the reaction rate was replaced with the reaction frequency ν calculated for the processes under consideration as

${\nu }_{\mathit{ion}}=n_e{\sigma }_{ion,i\to i+1}{v}_e,(28)$

${\nu }_{\mathit{CEX}}=n_0{\sigma }_{CEX,i\to i-1}{v}_i,(29)$

where n_e and n₀ are the electron and neutral ion density, respectively, and v_e and v_i are the electron and ion collision velocity in CM frame of reference, respectively. Naturally, the higher the frequency, the more the probability of the reaction taking place in a given time interval. The interplay of collisional ionization and charge exchange in the context of ECR plasma ion population kinetics is shown in Figure 4. Exchange cross sections are independent of the ion energy below 25 keV/u (which is orders of magnitude higher than ECR ion energies) so the cross sections are reported in Figure 4A for each charge state following Eq. 27. The ionization cross sections are, instead, reported for a range of electron energies well into 100 keV, and it can be clearly seen that the values decrease with the increasing charge state and always remain below the exchange cross sections. On the other hand, converting σ into ν, it can be seen that high n_e and E_e initially drive up the reaction frequency for low-charge states (Figure 4B) but are soon overtaken by exchange reactions with the intersection charge state depending on the pressure of gas in the system. This is an important point because it implies that highly charged ions are shuffled into lower states, leading to an overall decrease in the mean charge of the plasma.

FIGURE 4

FIGURE 4. (A) Single ionization cross section as a function of electron collision energy for three different charge states of Ar ions (top) and charge exchange cross section for different charge states taking ion temperature k_BT = 0.5 eV (bottom) and (B) frequency of ionization under most optimum electron conditions and frequency of exchange for two different gas pressures as a function of charge state. The shaded area denotes the charge state after which CEX dominates over ionization, effectively lowering the plasma CSD.

The main task at hand was to develop a MC routine that would allow expressing CEX, ionization or the lack of any reaction thereof as unique and mutually exclusive outcomes of a single “effective” collision marked by a summed reaction frequency ν_tot = ν_ion + ν_CEX, so as to simplify and generalize the sampling procedure. The formalism was akin to modeling absorption by a material composed of two absorbers of coefficients, namely, μ₁ and μ₂. The transmission probability through such an element would be simply $T={\mathrm{e}}^{-({\mu }_1+{\mu }_2)x}$ and the corresponding probability distribution function would then be $p_{\mathit{abs}}=({\mu }_1+{\mu }_2){\mathrm{e}}^{-({\mu }_1+{\mu }_2)x}$. Hence in our case the probability distribution function is taken as

$p_{\mathit{tot}}\left(t\right)={\nu }_{\mathit{tot}}{\mathrm{e}}^{-{\nu }_{\mathit{tot}}t},(30)$

and thus in a given time interval T_step, the probabilities associated with no reaction and some reaction are

${\bar{P}}_{\mathit{tot}}\left(T_{\mathit{step}}\right)={\int }_{T_{\mathit{step}}}^{\infty }{\nu }_{\mathit{tot}}{\mathrm{e}}^{-{\nu }_{\mathit{tot}}t}\mathrm{d}t={\mathrm{e}}^{-{\nu }_{\mathit{tot}}{T}_{\mathit{step}}},(31)$

$P_{\mathit{tot}}\left(T_{\mathit{step}}\right)={\int }_0^{T_{\mathit{step}}}{\nu }_{\mathit{tot}}{\mathrm{e}}^{-{\nu }_{\mathit{tot}}t}\mathrm{d}t=1-{\mathrm{e}}^{-{\nu }_{\mathit{tot}}{T}_{\mathit{step}}}.(32)$

Furthermore, the probabilities associated with ionization and CEX were derived from P_tot(T_step) as

$\begin{array}{cc}\hfill P_{\mathit{tot}}\left(T_{\mathit{step}}\right)& =1-{\mathrm{e}}^{-{\nu }_{\mathit{tot}}{T}_{\mathit{step}}}=\frac{{\nu }_{\mathit{ion}}+{\nu }_{\mathit{CEX}}}{{\nu }_{\mathit{tot}}}\left(1-{\mathrm{e}}^{-{\nu }_{\mathit{tot}}{T}_{\mathit{step}}}\right)\hfill \\ \hfill & =\frac{{\nu }_{\mathit{ion}}}{{\nu }_{\mathit{tot}}}\left(1-{\mathrm{e}}^{-{\nu }_{\mathit{tot}}{T}_{\mathit{step}}}\right)+\frac{{\nu }_{\mathit{CEX}}}{{\nu }_{\mathit{tot}}}\left(1-{\mathrm{e}}^{-{\nu }_{\mathit{tot}}{T}_{\mathit{step}}}\right)\hfill \\ \hfill & =P_{\mathit{ion}}\left(T_{\mathit{step}}\right)+P_{\mathit{CEX}}\left(T_{\mathit{step}}\right)\hfill \end{array}.(33)$

This new formalism mapped the continuous probability distribution function to a discrete probability distribution with three outcomes which sum to 1 as expected. By sampling a single random number r ∈ [0, 1], a decision was made between ionization, CEX, and nothing as

• Ionization if 0 < = r < P_ion

• CEX if P_ion < = r < (P_ion + P_CEX)

• Nothing if (P_ion + P_CEX) < = r < 1

The method is general in the sense that more the number of reactions considered, the more would be ν_tot and thus the smaller the probability of nothing happening. Meanwhile, the probability of individual reactions would vary depending on the weight $\frac{{\nu }_X}{{\nu }_{\mathit{tot}}}$. The method eliminated the need to implement sequential or conditional probabilities and allowed interpreting all possible outcomes, whether mutually exclusive or not originally, as “effectively” exclusive through integration into a single collision event.

The procedure is not without its demerits. As much as it simplifies the sampling, the reactions are not truly mutually exclusive. The source of each reaction is different physical systems and there is no rule that both reactions may not occur simultaneously or in rapid succession. The correct approach would be to devise a joint probability distribution function for a sequential application of the reaction probabilities which would require use of multiple sampling variables and/or more involved techniques such as Markov chains. The approximation proposed here can eliminate the problem of simultaneity as long as T_step is small enough such that only one reaction may occur in that period. Even if the actual reactions are independent, this would enforce mutual exclusivity within the chosen interval. Coincidentally, this is precisely what one performs in a rate matrix algorithm—choose a time step Δt small enough that the populations of involved species changes slowly. This reduces the chances of cross-talk and hence errors in the final results. Numerically, this accounts for the discretization error when solving a differential equation. In our model, mutual exclusivity is naturally enforced because T_step is fixed to 10^–10 s as a conservatively small percentage of the Larmor gyration time $(\sim 10^{-6}\mathrm{s})$ and characteristic reaction times in Figure 4 are 2–3 orders of magnitude higher.

7 Convergent density scaling and double layer field

As already mentioned in Sections 3, 4, the trajectories of macroparticles were used to update a 3D matrix that served as an occupation or accumulation map representative of the true physical density. For ions, devising a robust scaling procedure was of utmost importance because charge states were sequentially simulated and the relative weights of all preceding occupation maps were constantly changing. The objective was to allow the dynamic weights to converge to the real steady-state value and thus correctly reproduce the plasma CSD. For this, it was essential that the ion transport, and consequently the occupation maps, be modeled accurately.

Unlike in electron simulations, the ions moved in a self-generated field and thus the transport coefficients were only as precise as the total scaled density n_s. This interdependence between scaling and transport implied that a rapid and correct convergence of the weights of the various charge states could only occur if the occupation maps were scaled regularly enough. It was, thus, decided to extract the occupation maps at specific intervals before T_span was reached, denoted by the checkpoints τ_p in Section 4. This helped avoid over/underestimation of the contribution of the various charge states toward n_s while remaining relatively fast to execute. Initially, it was decided to set τ_p as 10% of the effective confinement time τ_d given as the inverse sum of characteristic diffusion and magnetic confinement times [37] but this greatly increased the number of extraction points. As such for the sake of this preliminary study, it was decided to set four checkpoints at τ_p as 1%, 10%, 70%, and 100%T_span for each charge simulation.

Density scaling also allowed adding the electrostatic field arising from the double layer potential into the ion transport model. At a steady state, ECR plasmas separate into an electron-rich ellipsoidal shell in the interior called the plasmoid and a rarified zone surrounding it called the halo. The boundary between them is the ECR surface containing fewer, but hotter, electrons than the plasmoid. Electrons diffusing out of the plasmoid are reflected back inside while ions continue entering the halo—this creates a space–charge separation. A layer of positive ions forms just outside the electron layer and the potential generated from the double layer balances the flow of charges to maintain global neutrality. The existence of the double layer and subsequent potential dip has been a matter of debate, with some experiments even pointing to its absence [38]. Even the first results from Mironov et al inferred the same [19] but later models confirmed their existence [22]. The double layer has been theoretically modeled by Mascali et al [37, 39] and even experimentally verified by Takahashi, Kaneko, and Hatakeyama [40]. Thus, a dedicated electrostatic model was developed in COMSOL Multiphysics ^ⓒ to calculate ϕ_DL and E_DL (as mentioned in Section 4) from the charge imbalance ρ_Δ = ρ_tot − ρ_e using the Poisson equation, where ρ_tot is the total positive charge density evaluated using the scaled ion occupation maps. The effect of the electrostatic field E_DL is shown in Figure 5 where the ratio of ion and electron accumulation within a region of interest (ROI) defined deep in the plasmoid is plotted for regular intervals. The ROI was chosen as the space between two ellipsoidal surfaces in the plasma interior. If the electrostatic field worked as expected, the ratio of scaled ion density and electron density N_i/n_e would converge to a constant value so that the charge flows would be equalized. Although this has not been fully achieved yet, the higher values in the presence of the field in Figure 5B point to the confining action.

FIGURE 5

FIGURE 5. (A) Space between two centered and skewed ellipsoids in the plasmoid where ratio N_i/n_e was calculated as and (B) N_i/n_e in said space as evaluated at regular intervals with and without double layer E_DL field.

The scaling procedure was based on global charge neutrality which posited that the total number of positive and negative charges in the plasma would be the same, even if there existed localized charge pockets. The simple and straightforward expression to convert any general occupation map n_i into the corresponding number density map N_i is

${\oint }_V{n}_e\mathrm{d}V=K_i{q}_i{\oint }_V{n}_i\mathrm{d}V,(34)$

$N_i=K_i{n}_i,(35)$

where q_i is the charge state of the ion, K_i is the scaling factor, and the integrals represent volume integration. However at any checkpoint τ_p, there would be the additional presence of the ionization map n_i→i+1 and pre-simulated occupation maps n₁, … , n_i−1. Eqs 34, 35 were thus modified to account for n_i→i+1 and n_i−1 through coefficients representing the degree of ionization and exchange, respectively, and for other preceding occupation maps through weights emulating successive ionization from one charge state to the next. These are collectively referred to as transfer coefficients hereafter. As an example, at checkpoint τ_p for 1⁺ ions, density scaling was implemented as

${\oint }_V{n}_e\mathrm{d}V=K_1\left[\left(1-k_{1\to 2}\right){\oint }_V{n}_1\mathrm{d}V+2k_{1\to 2}{\oint }_V{n}_{1\to 2}\mathrm{d}V\right],(36)$

$N_1=K_1\left(1-k_{1\to 2}\right)n_1,(37)$

where $k_{1\to 2}=\frac{{\oint }_V{n}_{1\to 2}\mathrm{d}V}{N}$ and N is the initial number of macroparticles. The term in parenthesis indicates that the contribution of 1⁺ to the overall positive charge in the plasma was reduced by k_1→2 because this amount of macroparticles was ionized to the next charge state. As another example, at τ_p for the 3⁺ simulations, density scaling resembled.

$\begin{array}{cc}\hfill {\oint }_V{n}_e\mathrm{d}V=& K_3\left[\left(1-k_{1\to 2}+2k_{1\to 2}{k}_{2\to 1}+k_{1\to 2}{k}_{2\to 3}{k}_{3\to 2}\right){\oint }_V{n}_1\mathrm{d}V+\right.\hfill \\ \hfill & 2k_{1\to 2}\left(1-k_{2\to 3}-k_{2\to 1}+k_{2\to 3}{k}_{3\to 2}\right){\oint }_V{n}_2\mathrm{d}V+\hfill \\ \hfill & 3k_{1\to 2}{k}_{2\to 3}\left(1-k_{3\to 4}-k_{3\to 2}\right){\oint }_V{n}_3\mathrm{d}V+\hfill \\ \hfill & \left.4k_{1\to 2}{k}_{2\to 3}{k}_{3\to 4}{\oint }_V{n}_{3\to 4}\mathrm{d}V\right]\hfill \end{array},(38)$

$N_1=K_3\left(1-k_{1\to 2}+2k_{1\to 2}{k}_{2\to 1}+k_{1\to 2}{k}_{2\to 3}{k}_{3\to 2}\right)n_1,(39)$

$N_2=K_3{k}_{1\to 2}\left(1-k_{2\to 3}-k_{2\to 1}+k_{2\to 3}{k}_{3\to 2}\right)n_2,(40)$

$N_3=K_3{k}_{1\to 2}{k}_{2\to 3}\left(1-k_{3\to 4}-k_{3\to 2}\right)n_3,(41)$

where $k_{1\to 2}=\frac{{\oint }_V{n}_{1\to 2}\mathrm{d}V}{N}$, $k_{2\to 3}=\frac{{\oint }_V{n}_{2\to 3}\mathrm{d}V}{N}$, $k_{2\to 1}=\frac{{\oint }_V{n}_{2\to 1}\mathrm{d}V}{N}$, $k_{3\to 4}=\frac{{\oint }_V{n}_{3\to 4}\mathrm{d}V}{N}$, $k_{3\to 2}=\frac{{\oint }_V{n}_{3\to 2}\mathrm{d}V}{N}$ and N still represents the number of macroparticles. By regularly extracting simulated occupation and ionization maps and dynamically updating the transfer coefficients in the charge neutrality expressions Eqs. 36, 38, convergence toward true CSD was ensured. Additionally, keeping track of the transfer coefficients allowed simulating the same number of macropartices for each charge state while accurately gauging their contribution to the total positive charge, avoiding errors from low statistical counts in the process. The charge density was calculated from the number density as ρ_i = eq_iN_i whereas the total plasma density for friction and diffusion coefficients was n_s = ∑_iN_i. Once ρ_i was obtained, it was used to calculate the imbalance ρ_Δ = ρ_tot − ρ_e = ∑_iρ_i and passed to the electrostatic solver to evaluate 3D maps of double layer potential ϕ_DL. The raw ϕ_DL maps were peak normalized and scaled with a reasonable potential of 5 V [39] and the field components were re-evaluated using the gradient of the scalar field to obtain E_DL. These were passed back to MATLAB ^ⓒ to perform the ion transport until the next extraction checkpoint.

8 Description of code flow

Based on the discussion in the preceding sections, the detailed algorithm is presented below for the first two charge states and can be extrapolated similarly for the remaining.

8.1 1⁺ → 2⁺

• Initialization

– Generate N = 10⁶ macroparticles distributed throughout the plasma simulation domain. The initial positions would follow a pre-generated ionization map n_0→1 while the velocities would be generated according a Maxwell distribution of temperature k_BT_i.

– Use the electron density map ρ_e and energy density map E_e to calculate the ionization frequency for each plasma cell according to the expression ν_ion = ρ_eσ_ion,1→2v_e. Calculate plasma-averaged collision time ⟨τ_ion⟩ and set the total simulation time T_span = 10⟨τ_ion⟩ during which full ionization of the particles would be achieved.

– Define the magnetostatic field B configuration. Initialize occupation n₁ and ionization n_1→2 maps to 0 with same dimensions as the plasma simulation domain (59 × 59 × 211, 1 mm⁻³ cell volume). Also, initialize the lost particle map to 0. Set iteration step time T_step = 10^–10 s to increase precision. The value can be optimized later for more efficient computation.

– Set an estimate for the checkpoint τ_p at which the double layer field E should be evaluated. The importance of τ_p and evaluation method is discussed in Section 7.

• Evaluation

– Check for occurrence of ionization at each cell according to ν_ion. If the macroparticle is ionized, map to ionization matrix n_1→2. More details are given in Sec. 6.

– Apply the equation of motion to the remaining 1⁺ particles and update the position and velocity vectors. Check if the new position exceeds the plasma dimensions–if yes, update the lost particle matrix otherwise map the positions to the occupation map. The exact transport equations can be found in Section 5.

• Check extraction

– If checkpoint τ_p is reached, scale the occupation map n₁ obtained so far to appropriate density ρ₁ and use the charge imbalance map ρ_Δ = ρ_tot − ρ_e = ρ₁ − ρ_e as input to COMSOL Multiphysics ^ⓒ to calculate the double layer potential ϕ_DL and the corresponding field E_DL. Also, calculate the total plasma density N_s and mean charge ⟨Z⟩. More details on the methodology are presented in Section 7.

– Give E_DL as input back to MATLAB ^ⓒ and re-trace the motion of the 1⁺ macroparticles in presence of all fields. Update the occupation map n₁. Repeat from evaluation.

• Check termination

– If all particles are lost or T_span is reached, save n₁ and n_1→2 and final electrostatic field map E_DL.

8.2 2⁺ → 3⁺

• Initialization

– Generate N = 10⁶ macroparticles distributed in position according to n_1→2. Velocities remain Maxwellian with temperature k_BT_i.

– Use n_1→2 and n₁ to calculate the total plasma density n_s and mean charge ⟨Z⟩. At this point, the contribution of 1⁺ would be maximum. The methodology for calculation of N_s is given in Section 7.

– Use ρ_e and E_e to calculate the ionization frequency for each plasma cell ν_ion = ρ_eσ_ion,2→3v_e. Calculate cell-averaged collision time ⟨τ_ion⟩ and set the total simulation time T_span = 10⟨τ_ion⟩. Also, calculate the frequency of a single electron exchange with neutrals in each plasma cell, using the expression ν_CEX = ρ₀σ_CEX,2→1v_i. Save the total reaction frequency map as ν_tot = ν_ion + ν_CEX.

– Define the magnetostatic field B configuration. Initialize occupation n₂, ionization n_2→3, and charge exchange (CEX) n_2→1 maps to 0. Initialize lost particle map to 0. Set iteration step time T_step = 10^–10 s.

– Set a reasonable estimate for checkpoint τ_p at which density scaling and E_DL evaluation should be performed. More details are presented in Section 7.

• Evaluation

– Check for occurrence of reaction at each cell according to ν_tot. If the macroparticle is ionized, map to ionization matrix n_1→2 and if it is exchanged, map to CEX matrix n_2→1. More details are presented in Section 6.

– Apply the equation of motion to the remaining 2⁺ particles and update the position and velocity vectors. Check if the new position exceeds the plasma dimensions–if yes, update the lost particle matrix otherwise map the positions to the occupation map. The exact transport equations can be found in Section 5.

• Check extraction

– If checkpoint τ_p is reached, scale the occupation maps n₁ and n₂ to appropriate density maps ρ₁ and ρ₂, respectively, and pass the charge imbalance map ρ_Δ = ρ_tot − ρ_e = ρ₁ + ρ₂ − ρ_e as input to COMSOL Multiphysics^ⓒ to calculate ϕ_DL and E_DL. Also, calculate updated n_s and ⟨Z⟩—the contribution of 2⁺ would gradually increase. More details on the methodology are presented in Section 7.

– Give E_DL as input back to MATLAB ^ⓒ and re-trace the motion of the 2⁺ macroparticles in presence of all fields. Update the occupation map n₂. Repeat from evaluation.

• Check termination

– If all particles are lost or T_span is reached, save n₂, n_2→3, n_2→1, and final electrostatic field map E_DL.

9 Results and discussion

The results obtained from the application of the algorithm to the ATOMKI plasma trap configuration are presented in this section. A total of N = 10⁶ macroparticles were simulated for each charge state, using the magnetostatic field profile of Figure 1 and reaction cross sections from Eqs 26, 27. To be consistent with the electron maps, the simulation domain was taken as the 3D matrix of size 59 × 59 × 211 mm³, divided into cells of size 1 mm³ each. All ionization, occupation, and exchange maps were set to the same dimensions.

Figure 6 shows the distribution probability of the initial positions of the macroparticles along the X, Y, and Z axes for three different charge states. Since the initial positions were generated according to the ionization maps of the preceding charge state as mentioned in Section 4, and the plots are representative of the local energy content in the plasma. The initial positions of high charge states are peaked on the axis because this is where the more energetic electrons with sufficient ionizing power are located. The energy distribution along the axis is more uniform, as can be evinced from the near identical distributions in Figure 6C.

FIGURE 6

FIGURE 6. Normalized distribution of initial (A) X, (B) Y, and (C) Z coordinates of macroparticles for charge states 1⁺, 4⁺ and 7⁺.

A comparison of the initial and final positions of the macroparticles along the Cartesian axes is shown in Figure 7 for two charge states, namely, 2⁺ and 5⁺. The initial distributions follow the same trend as in Figure 6 but the final positions are starkly different, particularly for the Y and Z axes. On their journey through the plasma, the macroparticles tend to exit the plasmoid along the magnetic field lines and those that are not lost remain aggregated near the ECR surface, marked by the sharp peaks in Figures 7B,C. The peaks are obscured in Figure 7A because the X-axis does not pass through any magnetic poles and the position histograms are projections from neighboring hexapole axes.

FIGURE 7

FIGURE 7. Normalized distribution of initial and final (A) X, (B) Y, and (C) Z coordinates of macroparticles for charge states 2⁺ and 5⁺.

Moving on to the ion density profiles, a 3D map of regions with highest occupancy of various charge states is shown in Figure 8. Highly charged ions such as 6⁺ and 7⁺ tend to be localized near the axis rather than in the halo or on the ECR surface, whereas 1⁺–4⁺ are the opposite. The same is shown in Figure 9 where the 1D profile of the density of 1⁺, 4⁺, and 7⁺ along the coordinate axes is shown. Although these plots cannot be directly compared with results from other ECRIS simulations since the plasma configurations differ, the general trend of Figures 6, 7, 9 can be inferred from [20, 22] where the authors show tracks of Ar 1⁺ and 8⁺ in the plasma—the former (low charge) get regularly ionized in the near-axis region and leave traces around the plasma core leading to a spaced-out distribution while the latter (high charge) remain localized axially.

FIGURE 8

FIGURE 8. Region of peak occupancy of various charge states at the end of simulations until 7⁺ in a (A) sliced view and (B) X–Y projection view. The localization of high charge states near the axis is clearly visible.

FIGURE 9

FIGURE 9. Density profiles of 1⁺, 4⁺, and 7⁺ along the (A) X-axis at Y = 0, Z = 0, (B) Y-axis at X = 0, Z = 0, and (C) Z-axis at X = 0, Y = 0.

These peak occupancy maps can be correlated with the spatial distribution maps of the lost particles in Figure 10 where more diffused ions such as 4⁺ appear to be spread out while other charges leave finer tracks. Also, as with Figure 8B, Figure 10B shows that higher charge states remain stagnated near the axis are lost from a smaller region. A similar behavior has been reported in [22].

FIGURE 10

FIGURE 10. Scatter plots of macroparticles lost from the simulation domain in an (A) isometric view and (B) X–Y projection view for charge states 1⁺, 4⁺, and 7⁺. The downward pointed tri-star corresponds to the extraction side at Z = 0.09 m and matches the profile in Figure 8B.

In order to benchmark the simulations, the charge state dependent current extracted from the plasma during an experimental run at ATOMKI [41] was compared with that from the simulated density map. Figure 11A shows the experimental beam current as a function of the charge to mass ratio of the ions— Ar¹⁺ correspond to m/q = 40, Ar²⁺ to m/q = 20, and so on. The current under each peak I_i until Ar⁷⁺ was added to obtain the total current for every charge. To calculate the simulated current distribution, the mean density inside the region of potential dip ⟨N_i⟩ (Figure 12) was considered and converted to current using the following expression.

$I_i=\kappa \frac{\left(2L\right)S}{2}\frac{⟨N_i⟩q_{i}e}{{\tau }_i},(42)$

where κ is the transmission factor, L is the semi-plasma length, S is the area of the extraction hole, and τ_i is the confinement time. Given the high plasma density and existence of the potential dip, τ_i was calculated as

$\frac{1}{{\tau }_i}=\frac{1}{{\tau }_{ES,i}}+\frac{1}{{\tau }_{d,i}},(43)$

${\tau }_{ES,i}=R\frac{\sqrt{\pi }L}{v_i}\exp \left(\frac{q_{i}e⟨{\varphi }_{DL}⟩}{k_B{T}_i}\right),(44)$

${\tau }_{d,i}=7.1\times 10^{-20}L{q}_i\ln \mathrm{\Lambda }\sqrt{A}\frac{n_e{Z}_{\mathit{eff}}}{k_B{T}_i^{3/2}E},(45)$

where τ_ES,i is the electrostatic confinement time arising from the double layer, τ_d,i is the ambipolar diffusion time connected to the space–charge field along the magnetostatic field lines [38], R = B_max/B_min is the mirror ratio, ⟨ϕ_DL⟩ is the average potential inside the dip, kT_i is the ion temperature, ln Λ is the Coulomb logarithm, $Z_{\mathit{eff}}=\sum N_i{q}_i^2/n_e$ is the effective charge density of the plasma, and E is the ambipolar diffusion field. Taking ⟨ϕ_DL⟩ slightly smaller than what is directly calculated from the potential map, E ∼ kT_i, L = 0.105 m according to the chamber design and aperture diameter 18 mm, the current extracted from the plasma was calculated. Figure 11B shows the comparison between the experimental and simulated data. It can be seen that the general relation between the current and charge state is well reproduced, as is the absolute magnitude. This match is not obtained without considering the electrostatic confinement model, and while this does not solve the debate about the existence of the double layer by itself, it at least confirms the equivalence between the particle transport strategy and the resultant current. There is, however, an inversion in trends for 2⁺ and 3⁺ which may be partly due to the ion recycling effect [42] but more probably due to an oversimplified implementation of Eq. 6. The match is expected to improve when the simulation of all other charge states is completed, and when a more rigorous charge exchange model is implemented to account for interactions between all ions and not just neutrals.

FIGURE 11

FIGURE 11. (A) Experimental spectrum of ion current as a function of the mass–charge ratio and (B) comparison of experimental vs. simulated extracted current data until 6⁺. The blue (red) dotted line and axis correspond to the experimental (simulated) data. Simulations for higher charge states are still underway and will feed into 7⁺ state, hence it is not shown.

FIGURE 12

FIGURE 12. Norm of E_DL calculated at the end of simulation runs for 2⁺, 5⁺, and 7⁺ along the (A) X-axis at Y = 0, Z = 0, (B) Y-axis at X = 0, Z = 0 , and (C) Z-axis at X = 0, Y = 0, and (D) 3D map of ϕ_DL showing the potential dip confining the ions.

As a final remark, the profile of |E_DL| calculated by COMSOL Multiphysics ^ⓒ at the end of 2⁺, 5⁺, and 7⁺ simulations is shown in Figure 12, as is the potential map obtained at the end of 7⁺ run. The formation of the potential dip is noticeable and in keeping with the results of [22]. The convergence of the field profiles at the ECR layer after successive evaluations is clear, while the introduction of smaller peaks in the near-axis region is indicative of the localization of highly charged ions. A full comprehensive analysis of the electrostatic field behavior is still underway and will be a topic of discussion in future work.

10 Conclusion and future perspectives

A novel simulation scheme for modeling ion dynamics together with population kinetics in a space-resolved manner has been reported. The novelty does not arise from any breakthrough made in ECRIS simulations—as discussed in Section 2 numerous toolkits already exist to study the system—but from the adaptability of the technique which allows direct application to PANDORA. The algorithm has been developed in such a way as to exploit a pre-existing RF-electron coupling module and extend it to emulate the steady-state distribution of ions without starting from the scratch. The method simultaneously solves the CR-model to generate 3D space-resolved maps of level population and CSD consistent with steady-state electron maps. The scheme is general, makes few assumptions about the state of the plasma and can be extended to include reactions that populate ionic levels to generate the LPD needed in Eq. 4.

The algorithm presented earlier is still quite basic and needs improvement on several counts. For one, as already mentioned in Section 6, the effective mutual exclusivity approach only works as long as T_step is small enough to neglect simultaneous reactions. Even though this is true when considering only ionization and exchange, the approximation is expected to break down when increasing the number of reactions to include collisional excitation and spontaneous emission. The sampling procedure in such a case would inevitably require use of joint probability distributions invoking complex methods such as Markov chains. The development of such a model is underway. The implementation is also deficient in that the neutrals were considered uniformly distributed in the plasma chamber according to a fixed pressure, but this is not necessarily true because under stationary conditions the plasma core is almost devoid of neutral atoms due to ionization. The algorithm can, however, be very easily corrected for this effect, and the next version of the code will be free of this assumption.

Beyond this, as reported in Sec. 9, the precise self-consistent evaluation of |E_DL| is still ongoing because while the spatial profile seems to converge after sequential runs, the magnitude of the raw ϕ_DL maps as extracted from the electrostatic solver remain too high and need to be scaled for implementing in the particle pusher code. The procedure is similar to what Mironov and group have reported in [22] but their overall approach is more physically intuitive because the double layer formation in their model is a natural artifact of equalization of electron and ion diffusion flux. In future work, we aim to implement a similar calculation to directly evaluate correct field potentials. With these modifications in place, the algorithm may turn out to be a useful tool for not only fundamental research in the working of ECR plasmas but for applied research as well, similar to those needed for PANDORA.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

Conceptualization: BM, AG, AP, GT, GM, and DM; Data curation: SB, RR, EN, and DM; Formal analysis: BM, AG, AP, and DM; Investigation: BM, AG, AP, SB, RR, GT, GM, EN, and DM; Methodology, BM, AG, AP, and DM; Resources, AG, SB, RR, and DM; Visualization: BM, AP, and RR; Writing—original draft: BM; Writing—review and editing: AG, AP, SB, RR, GT, GM, EN, and DM.

Funding

The authors gratefully acknowledge the support of 3^rd National Committee of INFN and funding from project PANDORA_Gr3.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Geller R. Electron cyclotron resonance ion source and ECR plasmas. Bristol: IOP Publishing (1996).

Google Scholar

2. Girard A, Hitz D, Melin G, Serebrennikov K. Electron cyclotron resonance plasmas and electron cyclotron resonance ion sources: Physics and technology (invited). Rev Scientific Instr (2004) 75:1381–8. doi:10.1063/1.1675926

CrossRef Full Text | Google Scholar

3. Emery G. Perturbation of nuclear decay rates. Annu Rev Nucl Sci (1972) 22:165–202. doi:10.1146/annurev.ns.22.120172.001121

CrossRef Full Text | Google Scholar

4. Bahcall J. Beta decay in stellar interiors. Phys Rev (1961) 136:1143–9. doi:10.1103/PhysRev.126.1143

CrossRef Full Text | Google Scholar

5. Bosch F, Faestermann T, Friese J, Heine F, Kienle P, Wefers E, et al. Observation of bound-state β⁻ decay of fully ionized : - cosmochronometry. Phys Rev Lett (1996) 77:5190–3. doi:10.1103/PhysRevLett.77.5190${}^{187}\mathrm{R}\mathrm{e}$${}^{187}\mathrm{R}\mathrm{e}$${}^{187}\mathrm{O}\mathrm{s}$

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Mascali D, Santonocito D, Amaducci S, Ando L, Antonuccio V, Biri S, et al. A novel approach to β-decay: Pandora, a new experimental setup for future in-plasma measurements. Universe (2022) 8(2):80–16. doi:10.3390/universe8020080

CrossRef Full Text | Google Scholar

7. Takahashi K, Yokoi K. Nuclear β-decays of highly ionised heavy atoms in stellar interiors. Nucl Phys A (1983) 404:578–98. doi:10.1016/0375-9474(83)90277-4

CrossRef Full Text | Google Scholar

8. Galatà A, Mascali D, Gallo C, Torrisi G. Self-consistent modeling of beam-plasma interaction in the charge breeding optimization process. Rev Scientific Instr (2020) 91:013506–7. doi:10.1063/1.5130704

CrossRef Full Text | Google Scholar

9. Mascali D, Torrisi G, Neri L, Sorbello G, Castro G, Celona L, et al. 3d-full wave and kinetics numerical modelling of electron cyclotron resonance ion sources plasma: Steps towards self-consistency. Eur Phys J D (2015) 69:27–7. doi:10.1140/epjd/e2014-50168-5

CrossRef Full Text | Google Scholar

10. Heinen A, Ruther M, Ducrée J, Leuker J, Mrogenda J, Ortjohann H, et al. Successful modeling, design, and test of electron cyclotron resonance ion sources. Rev Scientific Instr (1998) 69:729–31. doi:10.1063/1.1148667

CrossRef Full Text | Google Scholar

11. Heinen A, Ruther M, Ortjohann H, Vitt C, Rhode S, Andra H. Heating and trapping of electrons in ecris, from scratch to afterglow. unpublished (1999) 1–8.

Google Scholar

12. Maunoury L, Pierret C, Biri S, Pacquet J. Studies of the ecr plasma using the trapcad code. Plasma Sourc Sci Technol (2009) 18:015019–7. doi:10.1088/0963-0252/18/1/015019

CrossRef Full Text | Google Scholar

13. Vamosi J, Biri S. Trapcad - a program to model magnetic traps of charged particles. Comput Phys Commun (1996) 98:215–23. doi:10.1016/0010-4655(96)00053-7

CrossRef Full Text | Google Scholar

14. Dougar-Jabon V, Umnov A, Diaz DS. Three-dimensional simulation of an ecr plasma in a minimum-b trap. Rev Scientific Instr (2002) 73:629–31. doi:10.1063/1.1429774

CrossRef Full Text | Google Scholar

15. Edgell D, Kim J, Bogatu I, Pardo R, Vondrasek R. Model of charge-state distributions for electron cyclotron resonance ion source plasmas. Phys Rev ST Accel Beams (1999) 2:123502–6. doi:10.1103/PhysRevSTAB.2.123502

CrossRef Full Text | Google Scholar

16. Edgell D, Kim J, Bogatu I, Pardo R, Vondrasek R. Electron cyclotron resonance ion source one-dimensional fluid modeling. Rev Scientific Instr (2002) 73:641–3. doi:10.1063/1.1427031

CrossRef Full Text | Google Scholar

17. Cluggish B, Zhao L, Kim J. Simulation of parameter scaling in electron cyclotron resonance ion source plasmas using the gem code. Rev Scientific Instr (2010) 81:02A301–3. doi:10.1063/1.3259166

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Cluggish B, Kim J. Modeling of wave propagation and absorption in electron cyclotron resonance ion source plasmas. Nucl Instr Methods Phys Res Section A: Acc Spectrometers Detectors Associated Equipment (2012) 664:84–97. doi:10.1016/j.nima.2011.10.015

CrossRef Full Text | Google Scholar

19. Mironov V, Beijers J. Three-dimensional simulations of ion dynamics in the plasma of an electron cyclotron resonance ion source. Phys Rev ST Accel Beams (2009) 12:073501–10. doi:10.1103/PhysRevSTAB.12.073501

CrossRef Full Text | Google Scholar

20. Mironov V, Bogomolov S, Bondarchenko A, Efremov A, Loginov V. Numerical model of electron cyclotron resonance ion source. Phys Rev ST Accel Beams (2015) 8:123401–23. doi:10.1103/PhysRevSTAB.18.123401

CrossRef Full Text | Google Scholar

21. Mironov V, Bogomolov S, Bondarchenko A, Efremov A, Kuzmenkov K, Loginov V. Simulations of ECRIS performance for different working materials. J Instrum (2018) 13:C12002. doi:10.1088/1748-0221/13/12/C12002

CrossRef Full Text | Google Scholar

22. Mironov V, Bogomolov S, Bondarchenko A, Efremov A, Loginov V, Pugachev D. Spatial distributions of plasma potential and density in electron cyclotron resonance ion source. Plasma Sourc Sci Technol (2020) 29:065010. doi:10.1088/1361-6595/ab62dc

CrossRef Full Text | Google Scholar

23. Mironov V, Bogomolov S, Bondarchenko A, Efremov A, Loginov V, Pugachev D. Three dimensional modelling of processes in electron cyclotron resonance ion source. J Instrum (2020) 15:P10030–19. doi:10.1088/1748-0221/15/10/P10030

CrossRef Full Text | Google Scholar

24. Mascali D, Neri L, Gammino S, Celona L, Ciavola G, Gambino N, et al. Plasma ion dynamics and beam formation in electron cyclotron resonance ion sources. Rev Scientific Instr (2010) 81:02A334–4. doi:10.1063/1.3292932

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Neri L, Mascali D, Celona L, Gammino S, Ciavola G. A 3d Monte Carlo code for the modeling of plasma dynamics and beam formation mechanism in electron cyclotron resonance ion sources. Rev Scientific Instr (2012) 83:02A330–3. doi:10.1063/1.3670341

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Torrisi G, Mascali D, Sorbello G, Neri L, Celona L, Castro G, et al. Full-wave fem simulations of electromagnetic waves in strongly magnetized non-homogeneous plasma. J Electromagn Waves Appl (2014) 28:1085–99. doi:10.1080/09205071.2014.905245

CrossRef Full Text | Google Scholar

27. Galatà A, Mascali D, Mishra B, Pidatella A, Torrisi G. On the numerical determination of the density and energy spatial distributions relevant for in-plasma β-decay emission estimation. Front Phys (2022) 10:1–9. doi:10.3389/fphy.2022.947194

CrossRef Full Text | Google Scholar

28. Angot J, Luntinen M, Kalvas T, Koivisto H, Kronholm R, Maunoury L, et al. Method for estimating charge breeder ecr ion source plasma parameters with short pulse 1+ injection of metal ions. Plasma Sourc Sci Technol (2021) 30:035018. doi:10.1088/1361-6595/abe611

CrossRef Full Text | Google Scholar

29. Mishra B, Galatà A, Mengoni A, Naselli E, Pidatella A, Torrisi G, Mascali D Predicting β-Decay Rates of Radioisotopes Embedded in Anisotropic ECR Plasmas. Nuovo Cimento C (2022) 45 (5):1–4. doi:10.1393/ncc/i2022-22117-5

CrossRef Full Text | Google Scholar

30. Chandrasekhar S. Principles of stellar dynamics. Chicago: The University of Chicago Press (1942).

Google Scholar

31. Boris J. Relativistic plasma simulation—Optimization of a hybrid code. In: Proceedings of 4th Conference on Numerical Simulation of Plasmas. Washington, DC: Naval Research Laboratory (1970). p. 3–67.

Google Scholar

32. Zenitani S, Umeda T. On the boris solver in particle-in-cell simulation. Phys Plasmas (2018) 25:112110–7. doi:10.1063/1.5051077

CrossRef Full Text | Google Scholar

33. MacDonald W, Rosenbluth M, Chuck W. Relaxation of a system of particles with coulomb interactions. Phys Rev (1957) 107:350–3. doi:10.1103/PhysRev.107.350

CrossRef Full Text | Google Scholar

34. Galatà A, Mascali D, Neri L, Celona L. A new numerical description of the interaction of an ion beam with a magnetized plasma in an ecr-based charge breeding device. Plasma Sourc Sci Technol (2016) 25:045007–17. doi:10.1088/0963-0252/25/4/045007

CrossRef Full Text | Google Scholar

35. Povyshev V, Sadovoyl A, Shevelko V, Shirkov G, Vasina E, Vatulin V. Electron impact ionisation cross sections of H, He, N, O, Ar, Xe, Au, Pb atoms and their ions in the electron energy range from the threshold up to 200 keV. IAEA INIS (2001) 33:1–46.

Google Scholar

36. Muller A, Salzborn E. Scaling of cross sections for multiple electron transfer to highly charged ions colliding with atoms and molecules. Phys Lett A (1977) 62A:391–4. doi:10.1016/0375-9601(77)90672-7

CrossRef Full Text | Google Scholar

37. Mascali D, Neri L, Celona L, Castro G, Torrisi G, Gammino S, et al. A double-layer based model of ion confinement in electron cyclotron resonance ion source. Rev Sci Instrum (2014) 85:02A511–4. doi:10.1063/1.4860652

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Douysset G, Khodja H, Girard A, Briand JP. Highly charged ion densities and ion confinement properties in an electron-cyclotron-resonance ion source. Phys Rev E (2000) 61:3015–22. doi:10.1103/PhysRevE.61.3015

CrossRef Full Text | Google Scholar

39. Mascali D, Gammino S, Celona L, Ciavola G. Towards a better comprehension of plasma formation and heating in high performances electron cyclotron resonance ion sources (invited). Rev Scientific Instr (2012) 83:02A336–6. doi:10.1063/1.3672107

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Takahashi K, Kaneko T, Hatakeyama R. Double layer created by electron cyclotron resonance heating in an inhomogeneously magnetized plasma with high-speed ion flow. Phys Plasmas (2008) 15:072108. doi:10.1063/1.2951997

CrossRef Full Text | Google Scholar

41. Biri S, Rácz R, Pálinkás J. Status and special features of the atomki ecr ion source. Rev Scientific Instr (2012) 83:02A341–3. doi:10.1063/1.3673006

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Mironov V, Bogomolov S, Bondarchenko A, Efremov A, Loginov V. Simulations of charge-breeding processes in ecris. In: Lettry J, Mahner E, Marsh B, Pardo R, Scrivens R, editors. Proceedings of 17th international conference on ion sources. College Park, Maryland: AIP Publishing (2011). p. 1–3.

Google Scholar