^{1}Laboratori Nazionali del Sud, Istituto Nazionale di Fisica Nucleare, Catania, Italy^{2}Dipartimento di Fisica e Astronomia “Ettore Majorana”, Università degli Studi di Catania, Catania, Italy^{3}Laboratori Nazionali di Legnaro, Istituto Nazionale di Fisica Nucleare, Legnaro, Italy^{4}Institute for Nuclear Research (ATOMKI), Debrecen, Hungary

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.

Here, *m*_{0}(*X*)*c*^{2} and *m*_{0}(*Y*)*c*^{2} represent, respectively, the rest mass energy of the parent and daughter system, while

where *K* and *K*′, *ϵ*^{i,j}; *ϵ*^{i}′^{,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*_{0} (∼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.

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

and the modified decay rate *λ* given as

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*_{0} 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**. **(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**. **(A)** 3D profile of simulated electron density (in m^{−3}) 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.

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

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**. **(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^{0+} 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

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.

Here, *n* is the index tracking the leapfrog procedure, **b** = **B**^{n+1}/|**B**^{n+1}| 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

and *R*_{rot} is defined as

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

where

Φ(*x*) is the error function and *G*(*x*) = [Φ(*x*) − Φ′(*x*)]/2*x*^{2}. *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

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

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

where *σ*_{ion,i→i+1} is the ionization cross section in m^{2}, *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

where *σ*_{CEX} is the exchange cross section in m^{2}; *A*, *α*, and *β* are fit coefficients; and *I*_{0} 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

where *n*_{e} and *n*_{0} 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**. **(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*_{B}*T* = 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, *μ*_{1} and *μ*_{2}. The transmission probability through such an element would be simply

and thus in a given time interval *T*_{step}, the probabilities associated with no reaction and some reaction are

Furthermore, the probabilities associated with ionization and CEX were derived from *P*_{tot}(*T*_{step}) as

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

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

## 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**. **(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

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*_{1}, … , *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

where *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.

where *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*_{i}*N*_{i} whereas the total plasma density for friction and diffusion coefficients was *n*_{s} = *∑*_{i}*N*_{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^{6} 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*_{B}*T*_{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→2}*v*_{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*_{1} and ionization *n*_{1→2} maps to 0 with same dimensions as the plasma simulation domain (59 × 59 × 211, 1 mm^{−3} 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*_{1} obtained so far to appropriate density *ρ*_{1} and use the charge imbalance map *ρ*_{Δ} = *ρ*_{tot} − *ρ*_{e} = *ρ*_{1} − *ρ*_{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*_{1}. Repeat from evaluation.

• Check termination

– If all particles are lost or *T*_{span} is reached, save *n*_{1} and *n*_{1→2} and final electrostatic field map **E**_{DL}.

### 8.2 2^{+} → 3^{+}

• Initialization

– Generate *N* = 10^{6} macroparticles distributed in position according to *n*_{1→2}. Velocities remain Maxwellian with temperature *k*_{B}*T*_{i}.

– Use *n*_{1→2} and *n*_{1} 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→3}*v*_{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} = *ρ*_{0}*σ*_{CEX,2→1}*v*_{i}. Save the total reaction frequency map as *ν*_{tot} = *ν*_{ion} + *ν*_{CEX}.

– Define the magnetostatic field **B** configuration. Initialize occupation *n*_{2}, 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*_{1} and *n*_{2} to appropriate density maps *ρ*_{1} and *ρ*_{2}, respectively, and pass the charge imbalance map *ρ*_{Δ} = *ρ*_{tot} − *ρ*_{e} = *ρ*_{1} + *ρ*_{2} − *ρ*_{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*_{2}. Repeat from evaluation.

• Check termination

– If all particles are lost or *T*_{span} is reached, save *n*_{2}, *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^{6} 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^{3}, divided into cells of size 1 mm^{3} 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**. 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**. 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**. 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**. 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**. 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^{1+} correspond to *m*/*q* = 40, Ar^{2+} to *m*/*q* = 20, and so on. The current under each peak *I*_{i} until Ar^{7+} 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.

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

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, *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**. **(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**. 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).

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

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

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

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

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

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

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

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

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

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.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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.

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

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

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

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.

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

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

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

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

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

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

Keywords: ECR plasma, ion CSD, PIC code, β-decay rates, MC routines, magnetic confinement, double layer

Citation: Mishra B, Galatà A, Pidatella A, Biri S, Mauro GS, Naselli E, Rácz R, Torrisi G and Mascali D (2022) Modeling space-resolved ion dynamics in ECR plasmas for predicting in-plasma *β*-decay rates. *Front. Phys.* 10:932448. doi: 10.3389/fphy.2022.932448

Received: 29 April 2022; Accepted: 05 September 2022;

Published: 03 October 2022.

Edited by:

Paul Stevenson, University of Surrey, United KingdomReviewed by:

Theocharis S. Kosmas, University of Ioannina, GreeceShixiang Peng, Peking University, China

Copyright © 2022 Mishra, Galatà , Pidatella , Biri , Mauro , Naselli , Rácz , Torrisi and Mascali . This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Bharat Mishra , mishra@lns.infn.it