R-matrix atomic data for application in astrophysics and fusion

A selection of electron-impact excitation R-matrix calculations completed by the collision group at Queens University Belfast are discussed. The atomic data calculated are the foundation of applications in a variety of fields, a sample of which are presented here. They include an analysis of Fe ii emission in active galactic nuclei, the possible presence of Pt and Au in the spectra of binary neutron star mergers and the viability of W ions as a diagnostic tool in magnetically confined tokamak plasmas. The radiative atomic data were computed using the most recent version of the fully relativistic structure code grasp0 and the collision cross sections were evaluated using the recently developed pdarcR-matrix package.


Introduction
Detailed, accurate and complete atomic data sets to include energy levels, radiative transition probabilities, excitation/deexcitation rates and photoionization cross sections are necessary for understanding and calculating stellar structure. Both quality and quantity of these atomic data are crucial for the accurate modelling of a host of astronomical spectra. While some of the data can be obtained experimentally through high-powered laser facilities (ASTERIX IV Germany, National Ignition Facility NIF USA) and Z-pinch plasma experiments (Sandia National Laboratory SNL USA), they are usually limited to a small number of transitions. These experimental measurements can, however, test the theoretical models to see if they are missing important collision processes such as photoionization from high-lying states or two-photon processes. The research group at Queens University Belfast has a longstanding and very successful track record in the evaluation of these vital data. They have also played leading roles in the development, maintenance and testing of current and new variants of the powerful R-matrix codes, on numerous platforms from local clusters to HPC facilities worldwide. The R-Matrix approach is credited as one of the most powerful and reliable tools in calculating these atomic data. Recent and ongoing developments of the relativistic parallel pdarc codes have enabled an order of magnitude advance in the accuracy of the atomic structure and subsequent collision calculations that are now feasible for lowly ionised high Z ions. In this topical issue we discuss some of the most recent Rmatrix calculations that have had a significant impact on several choice applications in astrophysics and fusion research.

Methodology
Throughout this topical review we will concentrate on heavy species with atomic number equal to or greater than Fe (Z=26). It is therefore prudent to utilize the fully relativistic package grasp0 (General Relativistic Atomic Structure Package) for the computation of the radiative data (energy levels, transition probabilities, etc) and the fully relativistic pdarc package (Parallel Dirac Atomic R-matrix Code) for the computation of the collision cross sections. Below is a summary of the theory underpinning both methods.

grasp0
The grasp0 package, originally developed by [1] but published in 1996 by [2], solves the multi-configurational Hartree-Fock equations and determines an optimised where φ is the Dirac orbital, and E corresponds to energy eigenvalues of the Hamiltonian. Within grasp0 we predominantly made use of the extended average level (EAL) method in which we give the diagonal elements of our Dirac-Coulomb Hamiltonian, defined (in atomic units) as weights proportional to (2J + 1), as A variational procedure then optimises the trace of the weighted Hamiltonian, and in turn, allows one to determine a set of atomic orbitals which describe closely lying states with excellent accuracy. In Eq. (2) α and β are related to the set of Pauli spin matrices, I 4 is the 4 × 4 identity matrix, Z is the atomic number, c is the speed of light, p is the momentum operator defined as p = −i ∇, r i denotes the position of electron i and r ij = | r i − r j | is the inter-electronic distance.

pdarc
A brief introduction to the R-matrix method will be included here, but the reader is directed to the fully detailed description given by Burke [3]. R-matrix theory divides configuration space into two distinct regions, the inner and outer regions. These regions are separated by an R-matrix boundary at r = a, chosen to completely encapsulate the most diffuse orbital and hence the charge distribution of the N -electron target. This boundary thus acts as an interface between the two regions. The R-matrix is defined as follows, where E N +1 k are the eigenenergies of the (N +1) Hamiltonian, E is the energy of the incident electron, and ω ik are the surface amplitudes. Within the internal region, electron exchange and short range correlation effects between the incident electron and the target must be accounted for. This is achieved by expanding the wavefunction, in terms of energy-independent wavefunctions ψ k , in the form where A Γ jk (E) are energy dependent coefficients and X N +1 = (x 1 , x 2 , ..., x N +1 ) where X i = (r i , σ i ) are the set of space and spin coordinates of the (N + 1) electrons. The energy-independent basis functions are written as close-coupling expansions of the continuum orbitals u ij and square-integrable correlation functions HereÂ is an antisymmetrisation operator and Γ = JM J π are the conserved quantum numbers. The remainder of the parameters listed in Eq. (5) are defined as follows: Φ Γ i are channel functions, u ij are functions which describe the continuum of the (N + 1) electron system for each value of angular momentum J, square integrable functions χ Γ i describe short range correlation effects, and the coefficients a Γ ijk and b Γ ik are determined from the diagonalisation of the (N + 1) electron Hamiltonian over the energy independent basis. Once the external region has been reached, the electron moves only in the long range potential of the target and we may neglect electron exchange and correlation effects. The wavefunction here takes the much simpler form where Φ Γ i are identical channel functions to those defined in Eq. (5) and the F Γ ij are reduced radial functions. The external wavefunctions from Eq. (6) can be matched to asymptotic boundary conditions given by [4] to allow for the determination of the collision strengths, Ω ij , for excitation from some initial level i to some final level j. These collision strengths are related to the cross section σ i→j by the relation, where g i is the statistical weight of the initial state, k 2 i is the energy of the incident electron in Rydbergs, and a 0 is the mean radius of the orbit of an electron around the nucleus of a hydrogen atom in its ground state. Effective collision strengths (Υ ij ) can also be determined by Maxwellian averaging over a Boltzmann distribution of electron temperatures so that, where the excited electron has energy j . The terms k and T e are Boltzmann's constant and the electron temperature (in Kelvin), respectively. It is these Maxwellian averaged effective collision strengths that Energy level diagram for each of the four Fe ii theoretical models considered; [6][7][8][9] are commonly used by astrophysical and plasma modellers in their diagnostic applications.

Recent calculations and applications
In this section we concentrate on a select few Rmatrix calculations which have an important astrophysical imperative. These examples serve to emphasise the difficulties encountered when dealing with heavy elements with atomic number greater than iron, while at the same time demonstrate the current capabilities of the atomic structure and collision codes.

Improved Fe ii emission-line models for AGNs
Emission lines of Fe ii contribute significantly to the spectra of AGNs (active galactic nulcei) from the IR to the UV wavelength regions. It is well known that the observed Fe ii spectra of the lower-luminosity narrowline Seyfert galaxies contain a 'UV bump' between the C iii λ1909 and Mg ii λ2800 emission lines. This bump is produced by a large number of blended lines due to transitions between high-lying states of Fe ii with energies above approximately 13.25eV. Recently, [5] incorporated three large Fe ii atomic data sets into the cloudy modelling code to predict AGN spectral energy distributions (SED's) and reproduce, for the first time, observed Fe ii emission with solar abundances. Fig. 1 displays an energy level diagram to pictorially represent the size of each of the theoretical models used to create the atomic data sets. The original work compiled by [6] contained 371 individual levels up to 0.94×10 5 cm −1 , however many of the collision strengths Fig. 2 A comparison between the observed Fe ii UV template of [10] and the cloudy predicted spectrum of [5] with V turb =100 km s −1 were generated from the 'g-bar' approximation. The [7] model incorporated 159 levels up to 0.93×10 5 cm −1 but unfortunately did not consider lines in the UV range (2000-3000Å) of relevance to the cloudy modelling in the [5] paper. The largest two, and the most recent, data sets were compiled by [8] (340 levels up to 1.3×10 5 cm −1 ) and [9] (716 levels up to 1.25×10 5 cm −1 containing a high density of autoionizing states). It was found that the [9] atomic data set with its higher density of states, produced significantly more emission in the FUV. Hence continuum and Lyα fluorescent excitation, known to be important for Fe ii, was better represented. The [5] cloudy modelling using the improved [9] data set was able to reproduce the Fe ii emission with solar abundances and, for the first time, the shape and strength of the Fe ii UV bump and optical emission. Figure 2 shows a comparison between the observed Fe ii UV spectra of [10] and the cloudy predicted UV spectrum. The agreement between observations and the predicted spectrum is excellent across the UV bump, indicating that a turbulent model reproduces the observations.

Pt and Au in the spectra of the Kilonova AT2017gfo
In 2017 the first gravitational wave from a binary neutron star merger (NSM) was detected as part of the LIGO [11] and VIRGO [12] projects. Neutron stars come in pairs and represent ultra-dense collapsed cores of stars, they ultimately collide at 20% of the speed of light. The ejected radioactive matter, which creates a bright glow via r-process nucleosynthesis, has been termed a kilonova and has provided us with an incredibly rich and complex set of spectra pictured in Fig. 3. The spectra consist of 10 X-shooter spectra which were flux calibrated to a compiled set of photometric measurements and are publicly available through the ENGRAVE project [13]. Clearly visible are 71 Page 4 of 9 Eur. Phys. J. D (2023) 77 :71 Fig. 3 ESO VLT X-Shooter spectra of the counterpart of GW170817 from [14] and [15], at phases indicated in days after merger time a myriad of unidentified emission and absorption lines which are believed to be due to r-process elements, elements such as the lanthanides, actinides and platinum group elements. Disentangling r-process abundances from these broad spectra of NSM is challenging and requires a high level of accuracy in the calculations of the ejecta opacity and the atomic data sets that are utilised in the modelling. A combination of the models with the observations should allow us to identify the most important r-process species present in these spectra and conclude whether these binary neutron star mergers are the forge or the origin of the heavy r-process elements with atomic number Z greater the Fe. Difficulties arise, however, due to the paucity of atomic data currently available in the literature for these heavy species. Vital atomic data, such as energy levels, Einstein A-coefficients, electron-impact excitation cross sections and excitation rates, simply do not exist in suitable quantities for inclusion in the modelling codes. This makes modelling difficult as the elements without complete data sets will be excluded from consideration. There has, however, been a flurry of activity in recent years by multiple research groups worldwide to rectify this. Platinum and gold are two of the most abundant third peak r-process elements predicted to be synthesised in binary neutron star mergers. Their origin is currently unknown and without complete atomic data sets the radiative transfer modelling for these elements is incomplete. In an attempt to predict the presence of Pt and Au in the spectra of the kilonova AT2017gfo, two calculations by [16] and [17] have recently been published.
McCann et al [16] reports on the computation of extensive relativistic atomic structure and electron-impact excitation collision rates for Au i-iii which were subsequently used in a collisional-radiative model to investigate line ratio diagnostics in NSM environments. To compliment this work, Gillanders et al [17] employed the computed atomic data for Au i-iii and Pt i-iii to generate photospheric and simple nebular phase model spectra for kilonova-like ejecta properties.
To emphasise the difficulties encountered with these high Z open 5d shell systems we present in Tables 1, 2, 3 and 4 the computed target state energies and transition probabilities for the lowest few levels of Pt iii and Au iii from a selection of theoretical works. These two ionization stages were chosen as there is no radiative data currently available in the nist database for either of these ions. For Pt iii the relativistic grasp0 energies of [17] in cm −1 are listed in Table 1 for the lowest four levels and one high lying 5d 2 6 s 2 5 D state, just a sample of this 1639 level calculation which included a total of 12 configurations. In 2022 [18] used the relativistic grasp2k code with 5 configurations to compute the same energies with significant differences noted although the ordering of the levels is similar. The 5 configuration HULLAC evaluations of [19] are listed in column 7 and it is found that these values lie between the other two theoretical predictions. In the final column the experimental values of [20] are presented for comparison. The experimental value for the J = 3 ground state split level is in best agreement with the grasp0 work, the J = 2 level is best described by the hullac calculation, and for the metastable 1 D J = 2 level, although all three agree quite well the grasp2k prediction is closest. For the high lying state shown significant discrepancies are evident and no experimental data were measured.  Table 2 we present M1  A-values for two low-lying transitions in Pt iii. For the J = 4 to J = 3 transition the grasp0 ΔE is just over 9000 cm −1 and for the grasp2k just over 8000 cm −1 . These differences will be carried through into the calculation of the A-values and scaled as a power three, so it is not surprising that the grasp0 result is 36% larger than the grasp2k value. For the second M1 transition considered from the 1 D 2 metastable level to the 3 F 2 ground state split term, the ΔE's are approximately 8000 cm −1 for grasp0 and considerably less 6600 cm −1 for grasp2k, so again the A-values will reflect these differences.
The situation for Au III, presented in Tables 3 and 4, is somewhat better. The grasp0 calculation is from [16] and comprises a 13 configuration 2170 level model. The FAC (Flexible Atomic Code) [21] model has also been computed by McCann in the same paper. Experimental measurements exist ( [22]) and on comparison a much better agreement between all three energy sets is evident for these few low lying levels. For the A-values we look at the strongest E1 dipoles from the 2 D e ground state split levels to the first 2 D o state, the 5d 8 6p. Much better agreement is evident with the theoretical values of [23]. For a plasma in LTE the level populations can be found by solving the Saha equation with the Boltzmann relationship, simple analytic formulae that depend on the energy and temperature. No fundamental collision cross sections are required in this analysis. For the non-LTE case, however, the full detailed rate equations must be solved to obtain the atomic level populations complete with collisional rates. This solution can take orders of magnitude more computing time than the LTE case and large-scale computing resources are required. We present in Fig. 4 a plot of the steady state population fractions divided by their temperature dependent LTE population fractions, as a function of electron densities ranging from 10 5 -10 15 cm −3 for the ground and first metastable levels of neutral gold, for three temperatures 5800, 11600 and 17400K. For the 6 s ground state, plotted in red, the largest deviations occur at very low densities regardless of the temperature, with the populations slowly approaching LTE values at higher densities. For the first metastable 6 s 2 J = 5/2 the devia- The computations were carried out by [16]. The forbidden line between the fine structure components of the ground state, 5d 9 2 D e 5/2 − 5d 9 2 D e 3/2 , is plotted in Fig. 5. Clearly evident are the mass of Rydberg resonances converging onto the target state thresholds in the collision strength presented in the left hand panel. The resolution of these resonances enhances the Maxwellian averaged effective collision strength for this transition, presented in the right panel, particularly for low temperatures. As a contrast an example of an E1 dipole allowed transition is presented in Fig. 6, 5d 9 5 D e 5/2 − 5d 8 4p...... Contributions from the high partial waves to the collision strength are significant for allowed lines of this type, with the collision strength rising for higher electron impact energies. This behaviour is clearly reproduced in the corresponding effective collision strength for this transition.

Non-invasive diagnostics for magnetically confined fusion
The final application discussed in this paper relates to the heavy species of relevance to magnetically confined fusion. It is globally recognised that magnetically confined tokamak fusion reactor research is crucial for future energy sustainability. The material commonly used to line the walls of the divertor region of these tokamaks is tungsten (Z=74) with its high-energy threshold for erosion and excellent thermal properties. W erosion is undesirable and must be characterised and for neutral and lowly-ionised W ions this is a non-trivial task. A 2010 publication by [24], while investigating and testing the cooling factor of W in tokamaks, found that quenching of the plasma increased by 20% for a W Fig. 5 Collision strength as a function of incident electron energy (Ryds) (left panel) and corresponding effective collision strength as a function of electron temperature (K) (right panel) for the 5d 9 2 D 5/2 − 5d 9 2 D 3/2 forbidden transition in Au iii. The R-matrix computations were carried out by [16] Fig. 6 Collision strength as a function of incident electron energy (Ryds) (left panel) and corresponding effective collision strength as a function of electron temperature (K) (right panel) for the 5d 9 2 D 5/2 − 5d 8 6p 2 F o 7/2 E1 dipole transition in Au iii. The R-matrix computations were carried out by [16] Fig. 7 Observed spectra from the CTH experiment (solid blue line) compared to the present theoretical results between 216 and 227 nm. Vertical green and red sticks are the relative PECs for W i [25] and W ii [26], respectively. The electron temperature was 92800K and electron density 10 12 cm 3 concentration of as little as 3.0×10 −5 , compared to a plasma without any impurities. The QUB development work on the pdarc R-matrix codes allows us now to address the challenges of W i-vii ions for characterisation in tokamaks. Recent theoretical calculations for W i [25], W ii [26] and W iv [27] have progressed this crucial work. Collisional radiative theoretical models using the atomic data produced in these three publications can give a predictive capability for providing new temperature and density diagnostics for tokamak plasmas. The atomic data necessary to compute the required photon emissivity coefficients (PECs) and produce a synthetic spectra to compare with experimental measurements, are energy levels (wavelengths), A-values and electron impact excitation rates (effective collision strengths). At present, existing models are either incomplete or of insufficient accuracy for diagnostic work, causing ambi- Fig. 8 Observed spectra from the CTH experiment (solid blue line) compared to the present theoretical results between 243-253 nm. Vertical green and red sticks are the relative PECs for W i [25] and W ii [26], respectively. The electron temperature was 92800K and electron density 10 12 cm 3 guity in the line identification of many levels in the tungsten spectrum. To illustrate the quality of the modelling that the new atomic data produces, we present in Figs. 7 and 8 synthetic spectra for two wavelength regions 216-227 and 243-253 nm. The blue lines are the observed spectra from the Compact Toroidal Hybrid (CTH) tokamak at Auburn University, and the vertical green and red lines are the photon emissivity coefficients (PECs) derived from the collisional radiative modelling using the W i atomic data of [25] (green) and the W ii data computed by [26] (red). The temperature of the plasma during the shot was 92800K and the density was 10 12 cm −3 . We see good correlation between the collisional radiative modelling and the experiment for several prominent lines, suggesting that these lines are excellent diagnostic candidates. On these two plots it is evident that in both wavelength regions we have identified areas where both W i and W ii are present in the same window region and more importantly are unblended. This proves very useful when investigating redeposition onto the divertor wall of the tokamak, particularly for unblended lines of two ionization species. In addition these spectral line diagnostics will help characterise the impurity influx predictions of W i and W ii in fusion relevant plasmas.

Conclusions
The calculations for high Z ions presented in this review paper have been made feasible by recent developments in the fully relativistic packages grasp0 and pdarc. These code developments include multi-level parallelism and memory management during run-time and were computed in house at QUB by the author CPB. These packages are freely available online at http://connorb.freeshell.org/ and are scaled to run on as little as one processor up to 100,000 on large scale supercomputers. The atomic data presented and discussed have significant use in a variety of areas, the three chosen as representative applications were Fe ii emission in active galactic nuclei, the possible presence of Pt and Au in the spectra of binary neutron star mergers and the viability of W ions as a diagnostic tool in magnetically confined tokamak plasmas. However, the atomic data computed will have relevance to many more research areas in astronomy and plasma physics.