Overview of Numerical Simulations for Calculating In-Plasma β -Decay Rates in the Framework of PANDORA Project

. β -decay rates are essential inputs in stellar nucleosynthesis models to explain observed nuclear abundances. While current models continue to use terrestrial values, experiments in storage rings indicate strong divergence be-tween decay rates of neutral and ionised atoms [1], necessitating renewed investigations into stellar decay rates. The PANDORA project aims at measuring lifetimes of specific radio-isotopes trapped in an ECR plasma (which mimics astrophysical environments to some extent) and compare them with theoretical predictions [2], consequently verifying the models and allowing decay rate estimation for any isotope in the stellar interior. We present here a simulation scheme to characterise the space-resolved charge state and level population distribution of bu ff er and radio-isotope ions in ECR ion sources in order to calculate in-plasma decay rates. The algorithm is based on a Particle-in-Cell Monte Carlo (PIC-MC) routine that simultaneously models charge transport with collision-radiative processes. Preliminary results from the simulation are also shown, along with important takeaways for code-optimisation.


Introduction
Electron cyclotron resonance ion sources (ECRIS) are devices to produce ion beams using dense and energetic plasmas. Their operation relies on the dual concepts of magnetic confinement and resonance heating wherein electrons energised through resonance interaction with microwaves are trapped in a magnetic bottle ( Fig. 1(a)) and collide with gas atoms to sequentially ionise them. Conventional ECR plasma chambers are filled with light "buffer" gases like H, He and Ar that either directly form the ions or "plasmise" injected heavy isotopes to generate radioactive ion beams. PANDORA (Plasmas for Astrophysics, Nuclear Decay Observations and Radiation for Archaeometry) is an upcoming facility at INFN-LNS designed to exploit the high density (n e ∼ 10 11 − 10 13 cm −3 ) and temperature (T e ∼ 0.01 − 100 keV) of ECR plasmas to emulate astrophysical environments and study phenomena like plasma-induced modification of β-decay rates [3]. PANDORA intends to directly measure T 1/2 of radioactive nuclei like 176 Lu and 94 Nb diffused in the ECR plasma through secondary γ-tagging and compare the values with model predicted decay rates [2]. Once verified, the model can be adapted to stellar density and temperature conditions, and decay rates evaluated even for those isotopes which do not emit secondary γ photons. * e-mail: mishra@lns.infn.it Theoretical calculation of in-plasma T 1/2 requires inputs on ion charge state distibution (CSD) and level population distribution (LPD) which are difficult to evaluate in ECRIS because of system non-homogeneity and anisotropy. Additionally, many different physical processes affect the CSD and LPD, necessitating self-consistent evaluation of plasma properties and reaction rates. The main task then is to develop suitable models to extract 3D spaceresolved CSD/LPD of buffer ions and radioisotopes, taking into account their transport in the plasma as well as reactions relevant to each. An overview of such an algorithm is reported in the following sections, whose details can be found in Ref. [5]. The algorithm is based on, and intended to be an extension of, an already existing steady state PIC code for modelling ECR plasma electron dynamics. The code is written in MATLAB © to take advantage of implicit matrix vectorisation offered by the language, and while it is mostly targeted towards use for PANDORA, its generic nature allows for expansion into fundamental ECRIS research as well.

Algorithm for Buffer Ions
The temporal evolution of population density in each charge state i is described by the standard ECRIS balance equation [4] dn where n i is the ion number density, n e is the electron density, γ represents the collisional ionisation rate, n 0 is the neutral gas density, E is the single-electron charge exchange rate and τ i is the characteristic confinement time of the charge state. In steady-state equilibrium, the time derivative of n i vanishes, and all gain/loss processes in the RHS of Eq. 1 can be calculated independently to reconstruct the distribution of n i . The simulation algorithm attempts to solve Eq. 1 using a PIC-MC code run sequentially for each i -N representative macroparticles are sampled according to n e γ i−1,i and moved in a collisional plasma containing electromagnetic fields ( Fig. 1(a)) according to Lorentz and Langevin equation. During their transport, n e γ i,i+1 and n 0 E i,i−1 are calculated and the trajectories of remaining particles mapped to accumulation maps. The accumulation maps of all charge states simulated until i are converted into density maps using charge neutrality to generate the CSD which converges to the true steady-state value as more charge states are simulated. By simultaneously modelling particle transport and atomic processes and tracking transfer between charge states, the algorithm captures the essential components of ECR ion dynamics and allows accurate reconstruction of density maps and ion CSD. The simulation scheme rests on 3 essential modules modelling, respectively, charge particle transport, MC-based interaction sampling and density scaling using transfer coefficients. Fig. 1(b) shows first results from the simulation in the form of 3D distribution of mean charge state -higher charges tend to accumulate near the axis, and all ions leave traces along the magnetic field lines on comparing with the field lines of Fig. 1(a).

Algorithm for Radioisotopes and Decay Rate
The PIC-MC scheme described in Sec. 2 can be identically applied to the radioisotope ions in order to calculate their 3D space-resolved CSD and LPD. The corresponding balance equation to be solved is modified to account for the disproportional concentration of the two ion species and a second set of equations added to model the distribution among levels i ′ within  Figure 1: (a) 3D magnetostatic field profile of typical ECRIS with resonance surface (green shell) and associated field streamlines (red curves) generated using solenoid (blue) and hexapoles (yellow) and (b) 3D map of mean charge of Ar plasma when simulated till 8 + . Operating conditions like microwave frequency, power and gas pressure are mentioned in insets. a charge state.
Here n R i represents the isotope ion density in charge state i while n R ii ′ is the same for level i ′ , n j is the buffer ion density distribution at the end of Sec. 2, and C and A respectively stand for collisional excitation and spontaneous emission rates. Eq. 2 differs from Eq. 1 in that selfcollisions between isotope atoms are neglected and charge exchange is evaluated through collisions with all charge states n i of the buffer ion. Eq. 3 is the added balance equation aiming to implement the so-called "coronal" model to describe level population distribution using electron collision excitation and spontaneous emission.
Once n R ii ′ is known, it can be converted to its equivalent probability distribution p R ii ′ and used to calculate the position-dependent decay rate λ tot as [2] where f IF(m) is the lepton phase volume and t 1/2 is the decay half-life measured on Earth, f * IF(m) is the in-plasma lepton phase volume with m denoting type of transition, W, W max and W min are energy variables connected to the Q-value Q(ii ′ ), F 0 is the Fermi function, S m (ii ′ ) is the shape factor, f d , f c are connected to the Fermi-Dirac distribution function, m e is electron rest mass, x(ii ′ ) denotes electron orbital configuration with charge state i and level i ′ , σ x is occupancy probability and f (x), g(x) represent the larger of electron wavefunction evaluated at nuclear surface. The summation in Eq. 4 indicates the contribution of multiple decay channels to total rate λ tot . Since the variation of λ tot with individual levels i ′ is small and depends more strongly on averages over several levels, grouping can be done based on type of decay, reaction energetics and kinetics to simplify Eq. 3 and gain computational time. This may be particularly useful for for heavy isotopes which have complex electronic configurations.

Future Upgrades
The simulation algorithm is being continuously upgraded from both physical and computational points of view. One of the main improvements involves updating Eq. 1 to include a more complete exchange module, one that considers collisions between charge states and modifies n 0 dynamically. Other tasks foreseen are inclusion of electron capture, collisional de-excitation and implementation of radiation transport equation (RTE) in all balance equations. There are also efforts underway to upgrade from a particle-in-cell to cloud-in-cell approach to enhance numerical precision, as well as attempts to port the algorithm to GPUs.