Opacities of Singly and Doubly Ionised Neodymium and Uranium for Kilonova Emission Modeling

Even though the electromagnetic counterpart AT2017gfo to the binary neutron star merger GW170817 is powered by the radioactive decay of r-process nuclei, only few tentative identifications of light r-process elements have been made so far. One of the major limitations for the identification of heavy nuclei is incomplete or missing atomic data. While substantial progress has been made on lanthanide atomic data over the last few years, for actinides there has been less emphasis, with the first complete set of opacity data only recently published. We perform atomic structure calculations of neodymium $(Z=60)$ as well as the corresponding actinide uranium $(Z=92)$. Using two different codes (FAC and HFR) for the calculation of the atomic data, we investigate the accuracy of the calculated data (energy levels and electric dipole transitions) and their effect on kilonova opacities. For the FAC calculations, we optimise the local central potential and the number of included configurations and use a dedicated calibration technique to improve the agreement between theoretical and available experimental atomic energy levels (AELs). For ions with vast amounts of experimental data available, the presented opacities agree quite well with previous estimations. On the other hand, the optimisation and calibration method cannot be used for ions with only few available AELs. For these cases, where no experimental nor benchmarked calculations are available, a large spread in the opacities estimated from the atomic data obtained with the various atomic structure codes is observed.We find that the opacity of uranium is almost double the neodymium opacity.


INTRODUCTION
About half of the elements heavier than iron are produced in the astrophysical rapid neutron-capture process (or r-process) (Cowan et al. 2021). In contrast to the s-process (slow neutron capture process), r-process neutron captures occur much faster than the typical beta-decay timescale of the synthesised nuclei (Burbidge et al. 1957). One of the most promising r-process production sites are binary neutron star mergers (BNS) (Symbalisty & Schramm 1982;Lattimer & Schramm 1974;Eichler et al. 1989;Freiburghaus et al. 1999), where the r-process is expected to lead to an electromagnetic transient known as kilonova (Metzger et al. 2010). Such signal was ★ E-mail: a.floers@gsi.de recently observed following the gravitational wave event GW170817 (Abbott et al. 2017a,b).
The observed signal of the kilonova AT2017gfo -the electromagnetic counterpart following GW170817 -suggests at least some heavy r-process material was produced (Kasen et al. 2017;Tanaka et al. 2017;Perego et al. 2017). The electromagnetic emission from the first and only spectroscopically observed kilonova was studied extensively in visible and near-infrared (NIR) bands observed by ground-based telescopes for about two weeks until it faded out of reach of the 8-10 m-class facilities (Shappee et al. 2017;Pian et al. 2017;Smartt et al. 2017;Cowperthwaite et al. 2017;Drout et al. 2017;Kasliwal et al. 2017;Tanvir et al. 2017). The tail of the quasibolometric light curve of AT2017gfo is in good agreement with the energy release rate from radioactive decays of r-process elements (Metzger 2017). The rapid spectral evolution from a blue and nearly featureless continuum to a red spectrum (peaking in the NIR), rich in absorption and emission lines as well as the evolutionary timescale of the light curve indicate that high-opacity elements must have been synthesised, with opacities much higher than those typical for the iron group elements (IGE) commonly seen in thermonuclear and core collapse supernovae (Smartt et al. 2017;Pian et al. 2017;Waxman et al. 2018).
So far, only a single element -strontium -has been firmly identified in the kilonova AT2017gfo (Watson et al. 2019;Gillanders et al. 2022). The location of the proposed Sr feature could also be explained by the He i 10831 Å line. However, to have a noticeable effect on the spectrum, the required helium mass exceeds what is expected to be produced in merger simulations by about one order of magnitude (Perego et al. 2022). In particular, no r-process elements of the second or third peaks have been unambiguously identified (see however Domoto et al. 2022, for tentative identifications of La iii and Ce iii). Unsuccessful searches involve cesium/tellurium (Smartt et al. 2017) and gold/platinum (Gillanders et al. 2021). Identification of spectral features of lanthanides or actinides, together with a constraint on the abundances, would settle the debate on whether BNS mergers are responsible for r-process nuclei seen in the universe, and, if so, whether they are the dominant site of production. A straightforward, albeit very challenging approach to element identification, is through radiative transfer modelling of the observed spectra of AT2017gfo, as was done with strontium (Watson et al. 2019).
Radiative transfer models in their simplest form (1D, local thermodynamic equilibrium [LTE], no time dependence) still require precise knowledge of level energies and bound-bound atomic transitions, which make up the bulk of the photon opacity in the r-process enriched ejecta. It is expected that lanthanide and actinide ions each have of order 10 6 relevant transitions -a factor of 10-100 more than IGE ions (Kasen et al. 2013) -while only a tiny fraction has been measured for a few selected ions (see e.g. experimental data on the NIST ASD; Kramida et al. 2021). The ejecta contain only neutral to ≈ 4 times ionised atoms at phases beyond 1 day. It is not experimentally feasible to measure a full set of opacities for all quasi-neutral ions from the IGE to the actinides (≈ 250 relevant ions), and thus the only way of obtaining complete atomic data is through theoretical atomic structure calculations. Since the observation of AT2017gfo, several calculations of weakly ionised r-process opacities have been published (Gaigalas et al. 2019(Gaigalas et al. , 2022Radžiūtė et al. 2020;Gillanders et al. 2021;Fontes et al. 2020;Tanaka et al. 2020;Pognan et al. 2022a;Silva et al. 2022;Carvajal Gallego et al. 2021), focusing mainly on lanthanides. However, if material with sufficiently low electron fraction (≈ 0.15 or lower) is ejected in the merging process, nucleosynthesis can proceed to the actinides, which are expected to have photon opacities similar or even higher than lanthanides. For actinides, only few atomic structure calculations have been performed (Even et al. 2020;Fontes et al. 2023).
In this paper we present new calculations of neodymium and uranium, which act as case studies of the lanthanides and actinides. We focus on the singly and doubly ionised ions due to their importance in the line-forming regions, based on kilonova models produced with radiative transfer codes. To investigate variations in the calculated atomic properties we use two atomic structure codes: the Flexible Atomic Code (FAC, Gu 2008) and the Hartree-Fock-Relativistic code (HFR, Cowan 1981), described in Sections 2.1 and 2.2, respectively. While we compute ab-initio atomic data with the HFR code, we calibrate the local central potential and the calculated level energies in the FAC calculations to experimental data (Sections 3.1 and 3.2). Finally, in Section 3.3 we discuss and compare the resulting atomic opacities to published data within the expansion opacity (data from Gaigalas et al. 2019;Tanaka et al. 2020), line binned opacity (data from Fontes et al. 2020Fontes et al. , 2023Olsen et al. 2020), and Planck mean opacity frameworks. Conclusions are drawn in Section 4.

FAC calculations
Part of the calculations for this work were performed using the open source and freely available Flexible Atomic Code (FAC) (Gu 2008) relativistic atomic structure package based on the diagonalization of the Dirac-Coulomb Hamiltonian given, in atomic units, as where and are the 4 × 4 Dirac matrices and accounts for potential due to the nuclear charge. Recoil and retardation effects are included in the Breit interaction in the zero-energy limit for the exchanged photon, whereas vacuum polarisation and self-energy corrections are treated in the screened hydrogenic approximation.
Relativistic configuration interaction (CI) calculations are performed based on a set of basis states, configuration state functions (CSFs), which consist of linear combinations of antisymmetrised products of one-electron Dirac spinors, where represents the spin-angular function and and are the radial functions of the large and small components, respectively. Successive shells are coupled using a coupling scheme. Atomic state functions (ASFs) Ψ are then constructed from a superposition of = 1, · · · , CSF configuration state functions (CSFs) with the same and parity symmetry, where the stands for the complete relevant information to define each CSF (configuration and coupling tree quantum numbers).
The mixing coefficients { } are obtained by solving the eigenvalue problem Hc = c, with c = ( 1 , 2 , . . . , CSF ) . The eigenvalues obtained from the diagonalization of the Hamiltonian matrix H are, therefore, the best approximation for the energies in the space described by the basis of the selected CSFs. Increasing the number of CSFs used would improve the wave functions, and hence, the expected accuracy of the AELs, but the improvement is not expected to be significant enough to outweigh the increasing computational cost. Therefore, we search for the optimal set of CSFs by examining how the level energies converge as the number of configurations increases.
FAC uses a variant of the conventional Dirac-Fock-Slater method to compute one-electron radial functions. The small and large components are determined by solving self-consistently the coupled Dirac equation for a local central potential ( ) where is the fine structure constant and are the one-electron orbital energies. Although not as accurate as the multiconfiguration Dirac-Fock method used in the GRASP2K (Jönsson et al. 2013) and MCDFGME (Desclaux 1975;Indelicato 1995) structure codes, this approach has many advantages. As the same potential is felt by all the electrons of the system, all orbitals are automatically orthogonal. Moreover, the secular equation necessary to determine the eigenvalues and the mixing coefficients has to be solved only once, making the calculations much faster and computationally efficient when compared to the other codes mentioned. FAC has been shown to provide particularly good results when compared to experiments and other structure codes for calculations on lighter and/or highly excited ions (Gu 2008).
Following a similar approach to the one first described by Sampson et al. (1989) and Zhang et al. (1989), a single fictitious mean configuration (FMC) with fractional occupation numbers is adopted. The local central potential is then derived from this mean configuration using a self-consistent Dirac-Fock-Slater iteration. This unique potential is then used in Equation 4 in order to determine the remaining one-electron radial orbitals. Although this reduces the computation time and avoids convergence issues, as the orthogonality of the different orbitals with the same -value is automatically ensured with this method, the potential is not optimised for a single configuration. To accommodate multiple configurations, the FMC approach offers a good compromise for reducing the overall error of the calculation and getting a satisfactory accuracy of individual level energies. Typically, for the construction of the FMC, the occupation of the active electrons is split equally between a set of configurations. However, the weight of each configuration can be changed in order to raise (or reduce) its contribution in the construction of the mean configuration.
A summary of the calculations achieved, for each ion, with FAC is given in Table 1. Multiple models were computed, with an increasing number of configurations in order to test the sensitivity of the level energy to a higher degree of correlation and of the opacity with the inclusion of higher lying states. The configurations were sorted by the lowest level's excitation energy and then successively included in the various models considered. The largest model considered included 30 configurations, for which they were used in relativistic configuration interaction calculations; however, due to computational limits, only oscillator strengths and wavelengths of radiative transitions between the lowest 22 configurations were included.
In all calculations using FAC, the potential was adjusted to reduce the difference between our computed values and the published data in the Atomic Spectra Database (ASD) of the National Institute of Standards and Technology (NIST) database (Kramida et al. 2021;Kramida 2022) for the case of Nd (Martin et al. 1978), while for U the Selected Constants Energy Levels and Atomic Spectra of Actinides (from now on abbreviated as SCASA), available as an online database (see Blaise & Wyart 1994), was used. This was done by manually changing the contribution of the configurations used in the construction of the FMC in order to better reproduce the experimental energy of the lowest levels for each and parity values. For each ion, roughly 100 values were tested before choosing the one that best matched the experimental data available. The final set of weights included in the calculations is also presented in Table 1.
Semi-empirical corrections to the energies can be added subsequently to correct for errors induced by the use of a mean configuration. However, when this FAC functionality is applied, we have noticed significant disparities between the energy levels computed and the reported experimental NIST values. Similar findings have also been reported in Lu et al. (2021) and McCann et al. (2022). As a final step, when possible, we use a calibration technique that moves all the levels of a given symmetry by the same energy shift, estimated by the difference of the calculated and measured excitation energies of the lowest level of the considered ( -) block. That calibration process, used for all computed ions, does not affect neither the orbitals, nor the ASF wave function compositions, but corrects the transition data through the transition energies and the partition functions through the excitation energies.

HFR calculations
The Hartree-Fock Relativistic code (HFR) was developed by Cowan (1981). In this computational approach, a set of orbitals is obtained for each configuration by solving the Hartree-Fock (HF) equations, which arise from a variational principle applied to the configuration average energy. Some relativistic corrections are also included in a perturbative way, namely the Blume-Watson spin-orbit (including the one-body Breit interaction operator), mass-variation and one-body Darwin terms. The self-consistent field method is used to solve the coupled HF equations.
In the Slater-Condon approach, the atomic wavefunctions (eigenfunctions of the Hamiltonian) are built as a superposition of basis wavefunctions in the representation, i.e.
The construction and diagonalization of the multiconfiguration Hamiltonian matrix is carried out within the framework of the Slater-Condon theory. Each matrix element is computed as a sum of products of Racah angular coefficients and radial Slater and spin-orbit integrals: In the present computations, scaling factors of 0.85 are applied to the Slater integrals, as recommended by Cowan (1981). It was recently shown that the choice of scaling factors between 0.8 and 0.95 virtually does not affect the computed expansion opacities . The eigenvalues and eigenstates obtained in this way can then be used to compute the radiative wavelengths and oscillator strengths for each possible transition.
All the configurations included in the multiple models used in our HFR computations are listed in Table 2. Unlike the strategy followed with FAC, a more conventional approach was used for our HFR calculations, in which single and double electron substitutions from reference configurations were included. In the case of U ii and U iii, several models were considered and tested, with an increasing number of configurations obtained by considering single and/or double electron excitations from reference configurations to higher orbitals with an increasing principal quantum number (each model includes Table 1. Configurations included in the FAC calculations for the different models tested. The presented configurations are sorted by energy. Boldface configurations are used for the radial optimisation of the potential. The weights associated with each configuration are also shown, ordered in the same way as the boldface configurations in the Configurations column. Consecutive rows for the same ion show the additional configurations included with respect to the base model.

Nd atomic data
The lowest energy levels for each parity and for the largest calculations performed with FAC and HFR for Nd ii and Nd iii are shown in Tables 3, 4, 5, and 6. The relative difference from the data available in the NIST ASD (Kramida et al. 2021) is also evaluated as Δ = | NIST − |/ NIST with = FAC, HFR. The HFR calculations give rise to 13 004 380 and 6 481 846 lines computed for Nd ii and Nd iii, respectively, whereas the FAC calculations yield a total number of 15 031 028 lines for Nd ii and 1 833 759 for Nd iii. Each level is labeled by the largest configuration in both coupling, which are extracted from the NIST database, when available, and coupling schemes, the latter obtained directly from the FAC output. While labels use the usual 2 +1 term symbols to represent the coupling of electrons, the label provided by FAC uses non-standard notation: a given subshell is denoted as ( + ) 2 or ( − ) 2 . Here and represent the usual principal and angular momentum quantum numbers, represents the number of equivalent electrons within the same subshell while 2 denotes 2 times the total angular momentum that the electrons in the subshell couple to. Finally, we adopt the + and − notation to denote = + 1/2 and = − 1/2, respectively. The total angular momentum of the level given by the coupling of all subshells J = J 1 + J 2 + . . . is also indicated as a subscript. As an example, the ground level of Nd iii is labeled as ((4 3 − ) 9 (4 1 + ) 7 ) 8in this case 3 equivalent electrons in a 4 5/2 subshell couple to give 2 1 = 9 while there is only one electron with = 7/2, and, therefore, 2 2 = 7. The level has a total angular momentum of 2 = 8.
Associated with the density of levels, and as reported in previous works (Gaigalas et al. 2019;Radžiūtė et al. 2020;Gaigalas et al. 2022) strong mixing between configurations is anticipated for both lanthanides and actinides. In such cases, the assignment of configurations to each individual level becomes more complicated. For this reason, levels present in the NIST ASD (Kramida et al. 2021) database were identified and matched to the theoretical values primarily based on their , parity, and their position within a -group. Although comparisons with the available experimental atomic energy levels should be reliable for the first few levels, the large gaps in the experimental data prevent direct comparisons of more excited levels, as direct matching to the calculated levels becomes unreliable for the reasons mentioned above. Using this approach, a total of 611 levels were identified for the calibrated 22 config + extra CI calculation of Nd ii in the NIST ASD. The average relative difference to experimental data is of Δ FAC % = 9.21% and Δ HFR % = 18.35% with mean energy differences of Δ FAC = 1970 cm −1 and Δ HFR = 3838 cm −1 .
For Nd iii, the 29 experimental levels available in the NIST ASD represent only 0.38% of the total number of theoretical levels (15230) computed with both codes used in this work. Therefore, it is difficult to precisely assess the accuracy of these calculations, in particular for the higher energy levels. We observe, however, a larger disparity between the two sets of calculations, with a mean accuracy of the FAC calculations similar to the Nd ii case, (Δ FAC % = 9.65%, corresponding to Δ FAC = 1975 cm −1 ), while a larger deviation from the experimental data is found for the HFR energies (Δ HFR % = 49.41% and Δ HFR = 9234 cm −1 ). Following the level identification described above, a one-to-one matching between FAC and HFR was possible for 27330 levels of Nd ii and 9668 levels of Nd iii. An average relative difference (evaluated as | FAC − HFR |/ FAC ) of 10.65% and of Table 2. Configurations included in the HFR calculations for the different ions and models tested.
A visual representation of the energy levels for the largest calculations achieved with FAC and HFR as well as their comparison to the data available in the NIST ASD can be found in Figure 1 for Nd ii and Figure 2 for Nd iii. In the case of the singly ionised ion, while some differences are found between the two codes used, particularly for 2 = 7 and 2 = 23, the dispersion of levels is similar, especially for the low-lying levels, where the level density is the lowest. Larger differences between HFR and FAC results are observed for Nd iii, particularly for the even parity. HFR calculations indeed predict a higher level density than FAC in the 60 000 -10 0000 cm −1 range: 20.23 levels/1000 cm −1 compared to a level density of 7.30 levels/1000 cm −1 from FAC. This may be explained by the presence of the three extra configurations considered in the HFR calculation with respect to the FAC one (i.e., 4 2 6 2 , 4 2 5 1 6 1 and 4 2 6 1 6 1 ), which give rise to 1200 extra levels that lie between 70 000 and 200 000 cm −1 . For the odd parity states the disparity is not as large. However, our HFR calculations give rise to levels with energies lower than those of FAC, with a constant shift of about ∼ 13 000 cm −1 . This could be explained by differences in the optimization of the orbitals in each of the codes. While in FAC the optimization of the potential is based on a restricted number of configurations chosen among all the configurations included in the physical model, the average energies are minimised for the whole set of HFR configurations. This difference between the two methods is thought to explain the lower energies computed by HFR for the levels belonging to configurations that are not optimised in the FAC calculations. The resulting effect on the opacity cannot be ignored (see Section 3.3 for a further discussion). Table 4. Nd iii lowest level energies (in cm −1 ) for -split groups and comparison with those from NIST (Kramida et al. 2021). Shown are the lowest bound states for each 2 -group, as well as the relative difference (in percent, columns Δ%) compared to NIST for the FAC Nd iii calculations involving 8 (both with and without potential optimisation), 15, 22 or 22 + extra CI configurations (see Table 1).  To assess the reliability of the relevant atomic data supplied, we have studied the convergence of the energy levels with the multiple models calculated using the FAC code. The lowest level energies for every and parity symmetries for the different calculations performed were compared to the data available on the NIST ASD. Results for singly-and doubly-ionised Nd are shown in Tables 3 and 4, respectively. In both tables we see that the potential optimization of FAC (described in Section 2.1) has the biggest impact on the energies levels, with improvements of 211% for Nd ii and of 19.5% for Nd iii.
The impact of the inclusion of a higher number of configurations is particularly noticeable for Nd iii, where the mean deviation from NIST recommended values of the lowest state of each -group decreased from 19.54% to 8.06%, with a greater impact on the odd parity states (where there was an overall improvement of about 20%).
Larger deviations (of 40-50% in some cases) were found for the lowest energy levels of each -group of Nd ii, even for the largest CI calculations based on 30 configurations (22 configurations + extra CI). The mean deviation to NIST data also increased from 16.20% to 19.30%. We should highlight, however, that the impact of the level energy precision on the opacities is not the same for all levels. Lower excited levels, in particular, should have a larger effect under the assumption of local thermodynamic equilibrium (LTE).
While direct effects on the opacity will be characterised in 3.3, a more in-depth convergence analysis of the calculated energy levels using FAC was carried out for the case of Nd ii where the average accuracy (with respect to NIST ASD experimental data) of the lowest energy levels is computed for different calculations. Figure 3 shows the effect of the different models in the resulting energy levels of the 7 lowest configurations of Nd ii. In particular, we highlight the large effect of the optimization of the local central potential, described in Section 2.1. This optimization has, by far, the biggest impact on the atomic data, and illustrates how sensitive the FAC calculations are to the chosen local potential and to the choice of the FMC. In any case, it is important to ensure the convergence of the energy levels with Table 6. Excitation energies (in cm −1 ) for the first 20 energy levels of Nd iii calculated for the largest models computed with the FAC and HFR codes, with matching experimental values available in the NIST database (Kramida et al. 2021) . Relative differences with NIST data are also shown (columns Δ FAC % and Δ HFR %, in percent). the inclusion of higher degrees of correlation through the inclusion of more configurations in the CSF basis set used in the CI scheme. Nonetheless, the impact of the additional configurations on the energy levels is on average very minor (typically only a few per cent), as can be seen from Figure 4. While the effect on the accuracy seems to be irregular, indicating that only 30 configurations may still not be sufficient to achieve convergence, we note that the inclusion of the extra configurations in the largest model, for which transitions were not computed, still resulted in a noticeable effect mainly on the more excited levels. The effect of the calibration applied in the last model has a more noticeable effect for the first even (4f 4 6 1 ) and odd (4f 4 5 1 ) configurations since the constant shift applied is adjusted to match the energy of the lowest levels of each -block.
When compared to the FAC calculations, larger deviations from experimental values are found for the AELs of Nd ii computed with HFR (greater than 50%) This is partly due to the fact that the predicted HFR ground state differs from observation and from the FAC calculations, whatever the correlation model (see Table 5). It is also worth mentioning that configuration average energy adjustments were tested to match the Nd ii HFR predicted ground level to the observation. While it was possible to obtain the observed ground state in the HFR atomic data when proceeding this way, the impact on the opacities computed with this set of data was negligible ).

U atomic data
Just as for the Nd atomic data, the lowest energy levels for each and parity for the most complex models used with both FAC and HFR for U ii and U iii are respectively shown in Tables 7, 8, 9 and 10. As no data, besides for the ground state, is available in the NIST ASD, the calculation obtained for the considered uranium ions are compared to data available in the Selected Constants Energy Levels and Atomic Spectra of Actinides (SCASA). The same methodology used for NIST data is used here to evaluate the relative difference from the data available in the SCASA database with Δ FAC = | SCASA − FAC |/ SCASA and Δ HFR = | SCASA − HFR |/ SCASA . With the HFR code, 12 758 946 and 7 177 574 lines are obtained for U ii and U iii, respectively, whereas the FAC calculations give rise to 11 605 793 lines for U ii and 2 549 511 lines for U iii.
As explained in Section 3.1, energy levels from the SCASA database were identified and matched to our computed values based on their -value, parity , and their position within a -group.
-labels are, in this case, taken directly from the SCASA database. Level identification for U ii and U iii has been carried out from the experimental work of Palmer and Engleman Jr. , together with calculations of Blaise and co-workers (Blaise et al. 1984(Blaise et al. , 1987Blaise & Wyart 1992). Unfortunately, a number ofvalues are still ambiguous, with a few levels from the original analysis still to be confirmed, particularly for the doubly-charged ion. Theoretical interpretation of these levels is still lacking for several of the measured levels. For this reason, unidentified levels or levels which predicted leading component has a percentage of less than 25% are only labelled by their electronic configuration, i.e. omitting the (unidentified) term symbol.
FAC labelling follows the same format as the one introduced in Section 3.1. Some discrepancies in the configuration of the leading component of the eigenvector have been observed between labels. While some variation is to be expected due to the use of different coupling schemes, we notice the large uncertainties in the identification of levels for the uranium spectra, which are partly related to the strong mixing and high level density, even at lower energies. Such issues may not directly affect the calculation of opacities, but they should be considered, particularly when assessing the quality and reliability of atomic data.
As for Nd, except for the first few levels, the lack of data for more excited levels makes the matching of the latter to our calculated values less reliable. For U ii, the average relative difference to SCASA data is of Δ FAC % = 28.44% and Δ HFR % = 31.34% with mean energy differences of Δ FAC = 7988.46 cm −1 and Δ HFR = 8867.60 cm −1 . In the case of U iii, agreement between the reliably matched SCASA data seems to favour the HFR calculation, for which no calibration on experimental data was performed and for which we find an average relative difference of Δ HFR % = 25.10%, and mean energy different to the data of Δ HFR = 6041.28 cm −1 compared to the results obtained with FAC, Δ FAC % = 31.51% and Δ FAC = 8275.58 cm −1 . We find average relative differences (measured as | FAC − HFR |/ FAC ) of 12.36% for U ii, with 26 844 identified and matched levels between the calculations of the two codes, and of 15.64% for U iii, with 9427 matched levels.
Comparison of the lowest levels in each -group with the levels from the SCASA database is also shown in Tables 7 and 8 for the multiple models used in the FAC calculations. The biggest relative differences are found for the lowest states of U ii, due to the fact our FAC results predict a different ground level to the one reported in the SCASA database. Similar to Nd, the inclusion of extra configurations has a bigger impact on the energy levels for the doubly-ionised case, with agreement to SCASA data improving to 7.46% for U ii and 183.93% for U iii when going from 8 to 30 configurations in the CSF basis sets, all FAC CI calculations including the optimization of the central potential. Just looking at the lower -states of each group, we do observe an improvement of 33% for U ii despite the fact that the FAC ground state differs from SCASA. For U iii, the optimization process was necessary to fix the wrong ground state. On average, while the agreement with measured data worsens with the optimization for the doubly-ionised ion (much due to the odd 2 = 12 state), the energies tend to converge to values closer to the experimental ones with the increasing number of configurations added in the model.
The levels obtained with the biggest computations performed with both FAC and HFR are visually represented in Figures 5 and 6 for U ii and U iii, respectively. A comparison with the available data from SCASA is also shown, highlighting the lack of available data that we intend to fill in the present work. The dispersion of levels obtained in both cases are in good agreement for singly-ionised uranium, even if a few differences can be observed, especially for the highest values of considered. Nevertheless, in the case of doubly-ionised uranium, the level density predicted by HFR seems to be higher than the FAC one below 50000 cm −1 for even parity, while we find a much better agreement for the odd parity states. As for Nd iii, the higher level density predicted by HFR in the even parity can partially be explained by the three extra configurations included in the HFR calculation and not in the FAC model (i.e., 5 2 7 2 , 5 2 6 1 7 1 and 5 2 7 1 7 1 ), which lead to 1200 additional levels of energies comprised between 4 000 and 130 000 cm −1 . In addition, as discussed in Section 3.1, the fact that the average energies of all configurations are optimised in HFR (while in FAC the optimization of the potential is based on selected configurations of the model) is consistent with the higher level density predicted by HFR at lower energies. The differences between the codes and their effects on the opacity will be further discussed in Section 3.3.
Despite the differences between the two sets of results obtained with FAC and HFR respectively, it is important to note that both calculations confirm the previous results from Silva et al. (2022) where a higher level density of low-lying levels of U iii was found when compared to Nd iii. Even with an identical atomic structure, this effect can be understood by the higher diffuseness of the 5 shell when compared to 4 (Cowan 1981). The larger radii (e.g. 5 compared to 4 ) increase the 5 − {6 , 7 , 7 } overlap, relatively to the overlap of 4 with the outer shells, increasing the levels density. While this effect is not strong enough to produce noticeable effects for U ii, similar effects are expected to appear in other actinide elements, particularly for ions with a lower number of electrons in the 5 shell. As for Nd ii, a configuration average energy adjustment procedure was tried for U ii and U iii to match the HFR predicted ground state to the observed one, but the impact on the computed opacities was insignificant ).

Bound-Bound Opacities
For the ions considered in this study (Nd ii, Nd iii, U ii, U iii), we compute the bound-bound opacities for conditions expected for the ejecta of NS mergers ∼1 day after coalescence. However, we caution the reader that the presented opacities serve as illustrations and as a means to compare with data provided by other groups rather than direct input for radiative transfer models. Line-by-line opacities, as well as frequency dependent and grey opacities, are sensitive to the composition of the ejecta, and thus, the ionisation balance for the given composition. Approximating a more realistic (e.g. solar) composition with just a single element (Nd or U) would yield opacities that could be different by several orders of magnitude.
There exist a number of frequency-dependent opacity formalisms in the literature. In this study, we focus on the widely used expansion opacity (see Eastman where is the time since merger, is the ejecta density, Δ is the bin width and is the transition wavelength with Sobolev optical depth For the Sobolev optical depth, is the oscillator strength of the line with transition wavelength and the corresponding lower level number density . The number densities of the lower levels can be computed by assuming LTE holds using the Saha ionisation and Boltzmann excitation equations where the indices and indicate ionisation states and level numbers, respectively, is the multiplicity of the level (2 + 1), ( ) are the partition functions and is the level energy. For the bin width, we chose 100 Å, but emphasise that neither the absolute scale nor the general shape of the expansion opacity depends on the bin width. Choosing a larger bin width has the same effect as smoothing the curves. In the literature, a bin width of 10 Å is often assumed. For infinitesimally small Δ → 0 the expansion opacity reduces to the 4f 4 6s 4f 3 5d 2 4f 4 5d 4f 3 5d6s 4f 4 6p 4f 3 5d6p 4f 3 6s6p   Table 7. U ii lowest level energies (in cm −1 ) for -split groups and comparison with those from Selected Constants Energy Levels and Atomic Spectra of Actinides (SCASA) (Blaise & Wyart 1992). Shown are the lowest bound states for each -P group, as well as the relative difference (in percent, columns Δ%) compared to SCASA for the FAC U ii calculations involving 8 (both with and without potential optimisation), 15, 22 or 22 + extra CI configurations (see Table 1).  Table 8. U iii lowest level energies (in cm −1 ) for -split groups and comparison with those from Selected Constants Energy Levels and Atomic Spectra of Actinides (SCASA) (Blaise & Wyart 1992). Shown are the lowest bound states for each 2 -group, as well as the relative difference (in percent, columns Δ%) compared to SCASA for the FAC U iii calculations involving 8 (both with and without potential optimisation), 15, 22 or 22 + extra CI configurations (see Table 1).  Table 9. Excitation energies (in cm −1 ) for the first 20 energy levels of U ii calculated for the largest models computed with the FAC and HFR codes, with matching experimental values available in the Selected Constants Energy Levels and Atomic Spectra of Actinides (SCASA) (Blaise & Wyart 1992) . Relative differences with SCASA data are also shown ( columns Δ FAC % and Δ HFR %, in percent). Sobolev line opacity. The Sobolev approximation (Sobolev 1960) is well justified for NS merger ejecta as the velocity gradient (homologous expansion) is steep and the corresponding expansion velocities (≈ 0.1 ) greatly exceed the thermal line width (of order 1 km s −1 ), assuming no particularly strong hydrodynamic instabilities are in effect. While for arbitrary velocity gradients, the Sobolev optical depth depends on , for a homologously expanding atmosphere ( = ) it simplifies to Equation 8.

LS label FAC label
The opacity depends on the physical conditions of the ejecta: temperature, density and time after coalescence. To facilitate comparisons with other studies we adopt typical parameters for the ejecta at 1 day after the merger. We choose = 5000 K in accordance with the continuum of AT2017gfo inferred from the spectrum at 1.4 days, as well as a density of = 10 −13 g cm −3 , characteristic for an ejecta mass of ∼ 10 −2 ⊙ distributed uniformly within a sphere expanding at 0.1 . The choice of temperature is motivated by the fact that in LTE the radiation temperature equals the plasma temperature. Additionally, in Fig. 9, we present expansion opacities for temperatures of = 4000 K and = 6000 K, representative for a fully singly or doubly ionised plasma, respectively.
In this study, we limit ourselves to LTE conditions only (see Pognan et al. 2022b;Hotokezaka et al. 2021, for a discussion of non-LTE effects on kilonova opacities). For the above-stated ejecta properties we present the resulting ionisation balance as a function of temperature in Fig. 7 for the atomic data computed with FAC and HFR, as well as published data using the GRASP2K and HULLAC atomic structure codes 1 . Table 10. Excitation energies (in cm −1 ) for the first 20 energy levels of U iii calculated for the largest models computed with the FAC and HFR codes, with matching experimental values available in the Selected Constants Energy Levels and Atomic Spectra of Actinides (SCASA) (Blaise & Wyart 1992) . Relative differences with SCASA data are also shown ( columns Δ FAC % and Δ HFR %, in percent).  Ionisation energies for all ions were taken from the NIST ASD study all atomic levels are used to compute the partition functions, the studies as mentioned earlier use only the ground state (M. Tanaka, priv. communication). At the temperatures prevalent in kilonova ejecta after ≈ 1 day, this can lead to a factor of order unity in the ratio of the partition functions entering the Saha equations compared to neglecting the temperature dependence through the Boltzmann factor (For Nd ii at 5000 K, the ratio of the partition functions Nd ii / Nd iii is 154.26 / 39.44 = 3.91, while taking only the ground states Nd ii,ground / Nd iii,ground gives 8 / 9 = 0.89). As the ion density also enters the expression of the Sobolev optical depth, presented opacities are expected to differ between this work and the aforementioned studies. We emphasise that the different partition function treatment only affects the transition region between ionisation stages, for example, for II↔III, between 4000 and 5500 K, as long as the Boltzmann level populations are properly normalised. All expansion opacities shown throughout this work, including those from published (Kramida et al. 2021). For ions other than singly or doubly charged, we used measured levels from NIST, so that the only changes between calculations are due to the atomic data of the singly and doubly charged ions. For temperatures in the range between 2000 K and 10000 K, typical for kilonova ejecta within the first week after merger, only levels up to a few eV can be populated and thus contribute to the partition function. We find that for the ions considered in this study, these low-lying levels are sufficiently well known such that the partition functions computed from measured (NIST ASD or SCASA) and calculated levels agree within 5%.
For most highly charged ions only the ground state and the corresponding statistical weight are known. Unknown excited levels of highly charged ions (⪆ 5+) have a negligible effect on the computed ionisation balance for temperatures between 2000 and 10000 K. We find that for all sets of atomic data considered in this study, the ionisation balance is in good agreement for temperatures corresponding to the continuum of AT2017gfo within the first week (≈ 2500 K to 8000 K). However, we note that for the often used temperature of = 5000 K the balance for both Nd and U falls on the steep slope between singly and doubly ionised, where a small difference in the atomic data -for example from the partition functions -can have a significant effect on the ionisation balance and thus on the resulting opacity. Fig. 8 shows a comparison of the ionisation balance of Nd and U using the atomic data from our calibrated 22 config + extra CI FAC calculations. We find that going from the lanthanide Nd to the actinide U shifts the transition from singly to doubly ionised states to slightly higher temperatures (Δ ≈ 130 K). As a result, when the temperature of the kilonova ejecta drops, Nd iii will start to recombine shortly before U iii. As we did not compute atomic data for neutral and triply ionised ions we cannot make any quantitative claims about the I ←→ II and III ←→ IV ionisation transitions.  . Energy levels for U ii for even (top) and odd (bottom) parity for the models with the larger number of configurations using both FAC and the HFR codes, depicted in red and blue horizontal lines, respectively. Black horizontal lines shown experimental data from the Selected Constants Energy Levels and Atomic Spectra of Actinides (Blaise & Wyart 1992). Darker colors show a higher density of levels in that region.

Expansion Opacity
Following Equations 7 to 10, we compute expansion opacities for Nd and U (Fig. 9). To emphasise the effect of the temperature on the computed opacity we show the expansion opacity for 4000 K, 5000 K and 6000 K, corresponding to the fully singly ionised, partially singly and doubly ionised and fully doubly ionised cases, respectively (see ion fractions of Fig. 7). We find that the atomic data calculated by GRASP2K and HULLAC as well as by our FAC models typically yield lower opacities (especially in the shorter UV wavelength range) than the atomic data from our HFR models. This is evident from the level structure, where closeto-ground states in the former have higher excitation energies than the ones in the latter, leading to a reduced opacity due to the exponential Boltzmann factor exp(− / ). Additional configurations included in the HFR models and, as discussed in Sections 3.1 and 3.2, the different optimization procedures used in the FAC , GRASP2K , HULLAC and HFR codes (which are based on a restricted number of configurations in the first three and on all the configurations included in the model in the last one) could explain such differences.
Peak expansion opacities of Nd, for which published atomic data is available, differ by about 0.5 dex among the calculations considered. While there is some variation at longer wavelengths due to differences in the level density near the ground state, the total number of configurations included in the atomic structure calculations particularly affects the opacity red-wards of the ionisation edge at shorter wavelengths of about 1200Ȧngstroms. In this region the sharp drop in opacity occurs at longer wavelengths for the GRASP2K and HULLAC curves (see Fig. 12 for a clearer view of this portion of the opacity spectrum). This is due to highly energetic states not present in smaller calculations, which cannot be populated in LTE, but serve as upper levels for transitions from the populated nearground levels. These transitions only contribute opacity for photon energies − ≈ ion . Fig. 10 illustrates the convergence of the opacity with an increasing number of configurations treated in our FAC and HFR atomic structure calculations. While the opacity at the long-wavelength tail remains virtually unchanged for all ions, with only small variations from a small shift of the involved energy levels occurring, additional opacity components come into play at short wavelengths. The location and magnitude of these components depend on the ions and calculation strategies involved.
For the Nd atomic data computed with FAC we notice an increase in the opacity with increasing number of included configurations between 3000 Å (≈ 4 eV) and 1100 Å (≈ 11.5 eV) and between 2000 Å (≈ 6 eV) and 650 Å (≈ 19 eV), respectively, for Nd ii and Nd iii. The lower bounds are close to the ionisation energy for both ions (11.6 eV, and 19.8 eV, respectively). Increasing the number of included configurations beyond 22 (Nd ii) and 15 (Nd iii) only had a minor effect on the opacity. In either case, if applied to LTE kilonova radiative transfer modelling, including the 15 lowest configurations should suffice for capturing the majority of the opacity, as the radiation field even in the early kilonova evolution (≈1 day) is negligible at wavelengths shorter than about 2000 Å. A 10000 K blackbody, as observed in the earliest (0.5 day) spectrum of AT2017gfo (Gillanders et al. 2022), contains only about 3% of its flux at wavelengths shorter than 3000 Å. . Energy levels for U iii for even (top) and odd (bottom) parity for the models with the larger number of configurations using both FAC and the HFR codes, depicted in red and blue horizontal lines, respectively. Black horizontal lines shown experimental data from the Selected Constants Energy Levels and Atomic Spectra of Actinides (Blaise & Wyart 1992). Darker colours show a higher density of levels in that region.
configurations until convergence is achieved. Going beyond 15 configurations affects the opacity only marginally below 2000 Å or for photon energies above ≈ 6 eV. We note, however, that the calibration to measured levels (purple curves in Fig. 10) has a much greater effect on the U opacity than the Nd opacity. Our HFR calculations follow a different strategy, in which additional configurations are added based on the maximum shell number of the excited states. From Fig. 10 it is apparent that the inclusion of configurations with ≤ 7 is sufficient to capture the bulk of the atomic opacity, both around the peak and the tail. Transitions from near-ground-state levels to the = 7 shell dominate the opacity across the full optical and NIR wavelength range. Near the peak (1000 Å to 4000 Å),the opacity is almost unaffected by the inclusion of transitions to levels from the > 7 configurations.
The opacity tail at long wavelengths is a direct measure of the density of levels near the ground state. Again, as only levels close to the ground state can be populated at temperatures of a few 10 3 K, the upper level needs to be within ≈ 1 eV (corresponding to ≈ 1 m) to the lower level. Similar to the peak expansion opacities, we find significant differences between structure calculations. However, the opacities from the calibrated 22 config + extra CI FAC and the GRASP2K calculations agree exceptionally well both near the peak and on the tail, for both Nd ii and Nd iii.
Next, we explore whether actinides opacities are comparable to lanthanides, or possibly, even higher. Fig. 11 shows the expansion opacities of Nd and U computed at 5000 K. We find peak expansion opacities of Nd,FAC exp = 368 cm 2 g −1 and U,FAC exp = 736 cm 2 g −1 using FAC (an increase by a factor of 2) and Nd,HFR exp = 815 cm 2 g −1 and U,HFR exp = 1103 cm 2 g −1 (an increase by a factor of 1.35) for HFR. A similar behaviour can be seen in the atomic data from Fontes et al. (2020Fontes et al. ( , 2023 shown in Fig. 12. At temperatures of 4641 K (0.4 eV) the opacity of U is higher than that of Nd by about a factor of 4. While there is significant variation in the opacities of any given ion from different structure codes, the resulting opacity of U for a given code always seems to be significantly higher than that of Nd.

Line-binned opacities
Another frequency-dependent opacity formalism is the so-called linebinned opacity (Fontes et al. 2020(Fontes et al. , 2023, where is the ejecta density, Δ are the widths of the frequency bins containing lines with number densities of the lower level and oscillator strengths . This expression is obtained by replacing the line profile with a flat distribution across the corresponding bin (Fontes et al. 2020). Compared with the expansion opacity (Equation 7), it can be pre-computed as Equation 11 is independent of the expansion time . Fig. 12 shows the Nd and U line binned opacities from Fontes et al. (2020Fontes et al. ( , 2023 as well as the largest calibrated FAC and HFR calculations for temperatures equivalent to 0.4 eV and 0.5 eV (4641 K and 5802 K, respectively). For low photon energies (long wavelengths) we find good agreement between our calculations and the ones from Fontes et al. (2020Fontes et al. ( , 2023 Figure 8. Comparison of the ion fractions of neutral to triply ionised Nd and U as computed with FAC. Even though differences between the onset of ionisation stages of Nd and U appear to be small, due to the large gradient in the transition regions the opacity will be affected. opacities come closer to those of Fontes et al. (2020Fontes et al. ( , 2023 as they exhibit a higher density of levels near the ground state (lower level energies than what is reported on the NIST ASD). However, for high photon energies the published line-binned opacities from (Fontes et al. 2020) display a sharp cut-off at ≈ 10 eV, which is not seen in the FAC and HFR data. We identify the additional high-energy opacity Table 11. Planck mean opacity for a pure Nd or U plasma with = 10 −13 g cm −3 and = 1 day after the merger. is the Planck mean opacity averaged over the temperature range T = 2500 -7500 K and T = 5000 -10000 K, respectively.  Fontes et al. (2020Fontes et al. ( , 2023. In the high-temperature case, in which neodymium and uranium are doubly ionised, we find variations of about one order of magnitude between the various calculations. This is likely a result of the difficult calibration for Nd iii and U iii and the resulting uncertainty in the density of levels near the ground state (see discussion in Section 3.1).

Planck Mean Opacity
In addition to the wavelength dependent expansion opacity we also compute wavelength-independent 'grey' opacities in the form of the Planck mean opacity defined as Grey opacities are predominantly used for light curve modelling, in which the diffusion of photons due to opacity is the dominant effect. Use of the Planck mean opacity requires the radiation field to be a blackbody spectrum, which, in the case of kilonovae, is well justified during its early evolution (see blackbody fits in Pian 2021; Smartt et al. 2017;Watson et al. 2019;Gillanders et al. 2022). At the point at which transitions begin to dominate the spectrum -in AT2017gfo this happened ≈ 5 days after merger -the usefulness of the Planck mean opacity begins to break down. A comparison of the Planck mean opacities for atomic data computed in this work and from other sources (Gaigalas et al. 2019;Tanaka et al. 2020) is shown in Fig. 13. The peak of the Planck mean opacity is located for all ions considered at ≈ 4500 K. The peak location follows from the convolution of the blackbody spectrum with the frequency-dependent expansion opacity (Fig. 9). A lower temperature leads to an increased abundance of singly ionised ions which have a much larger opacity than their doubly ionised counterparts, but a blackbody that peaks farther in the red, where the opacity is decreasing with −1 (Silva et al. 2022). Even though the expansion opacity does not change significantly between 2500 K and 4500 K, the reduced overlap between the Planck spectrum and the wavelength-dependent opacity leads to a sharp drop of the Planck mean opacity below 4000 K. Above 5500 K the Planck mean opacity remains almost constant, as the doubly ionised state dominates the ionisation balance and the blackbody spectrum peaks in the optical (5000 Å to 10000 Å for temperatures at which doubly ionised Nd/U prevails). The local maximum near 8000 K followed by a sharp drop is due to the III↔IV transition near that temperature, surpassing the increasing overlap of the Planck spectrum with the peak of the opacity. ] T = 6000K Figure 9. Left panels: Expansion opacity of Nd at = 4000 K, 5000 K and 6000 K and for = 10 −13 g cm −3 and = 1 day with a bin width of 100Å. Shown is the atomic data computed with the FAC (red) and HFR (blue) codes as well as published atomic data using GRASP2K (black, Gaigalas et al. 2019) and HULLAC (green, Tanaka et al. 2020). In all panels we show the 22 config + extra CI FAC and the n=9 HFR calculations (see Tables 1 and 2). The bin width was chosen such that variations between bins in the NIR are smaller than the differences between codes. The choice of the bin width does not affect the absolute scale of the expansion opacities but only acts as a smoothing parameter. Right panels: Same as in left panels but for U instead of Nd.
We find mean variations in the Planck mean opacity be-

CONCLUSIONS
We computed energy levels and oscillator strengths for the singly and doubly ionised ions of neodymium and uranium, using the FAC and HFR atomic structure codes. In the calculations we included a maximum of 30 (FAC) and 27 (HFR) non-relativistic electronic configurations.
To obtain better agreement between computed and experimental level energies from the NIST ASD (Kramida et al. 2021 Table 1). Middle and bottom panels: Same as top panels but for U instead of Nd (middle panels) and for HFR instead of FAC (bottom panels).
FAC calculations. Employing this optimisation reduces the mean deviation from NIST data among the lowest 100 levels by a factor of two. An automatised methodology of optimisation is being tested at the time of writing this paper, which will be used to try to provide a more reliable, while still complete, set of atomic data for relevant lanthanide and actinide elements. (F. Silva 2023, In prep.). In addition to the potential optimisation, we calibrated the calculated level energies, split into groups of -, to experimental data. We also investigated the convergence behaviour of our atomic data with the number of included electronic configurations. We confirm that ≈ 15 configurations up to the = 7 (Nd) or = 8 shell (U) are sufficient to capture most the opacity that falls into the optical and NIR wavelength regions. The inclusion of additional configurations beyond those only yields minor contributions to the opacity close to the ionisation potential of the ions.
For each of the ions we computed bound-bound E1 transitions. From the calculated transitions we derive wavelength-dependent opacities assuming LTE conditions in the plasma. Included electron configurations were chosen such that internal convergence of the atomic opacities is obtained. We find good agreement between the optimised and calibrated FAC atomic data and calculations of Nd from Gaigalas et al. (2019), Tanaka et al. (2020) and Fontes et al. (2020). Our ab-initio calculations using the HFR code yield higher opacities than the corresponding FAC calculations. This larger opacity is a direct consequence of the higher density of levels near the ground state in the calculated HFR atomic data. Doubly ionised ions,  Figure 11. Comparison of the expansion opacities of Nd and U at = 5000 K, = 10 −13 g cm −3 and = 1 day computed with the FAC (red) and HFR (blue) codes and from published data from GRASP2K (black, Gaigalas et al. 2019) and HULLAC (green, Tanaka et al. 2020). for which less experimental data is available, display a higher opacity variance between codes than singly ionised ions. We identify the density of levels near the ground state as the main reason for the opacity difference between codes, while the total number of included levels and lines only has a minor effect on the opacity. The higher level density predicted by HFR is partially due to the additional configurations included in the HFR computations, but such extra configurations often give rise to levels, which, for typical temperatures, cannot be populated in LTE. Only transitions between these and near ground-state levels contribute to the opacity, predominantly at wavelengths close to the excitation energy of the upper level (≈ 1000 Å -3000 Å). However, the optimization procedures in HFR and FAC are fundamentally different, insofar as the average energies of all the configurations within the model are minimised in the former, while the potential optimization is based on a restricted number of configurations in the latter (as well as in the GRASP2K and HULLAC codes that were respectively used in Gaigalas et al. (2019) and Tanaka et al. (2020)). As a consequence, this difference can explain the lower energies of the predicted HFR levels, thus the higher density of levels near the ground state, leading to higher opacities.
In this study, we confirm that actinides, such as uranium, have an opacity that is a factor of a few higher than the opacity of the corresponding lanthanides, with ratios of 1.35 (HFR), 2 (FAC) and 4 (Los Alamos suite of atomic physics and plasma modeling codes, Fontes et al. 2020Fontes et al. , 2023 depending on the atomic structure code used. This result is of particular interest for the modelling of kilonovae, as the presence of even small amounts of actinides formed in the lowest-ejecta can have strong implications on the derived synthetic light curves and spectra. energy [eV] Line-binned opacities of Nd (upper panels) and U (lower panels) at = 10 −13 g cm −3 , = 1 day and = 4841 K (0.4 eV, left panels) and = 5802 K (0.5 eV, right panels), using re-binned published opacities from the Los Alamos suite of atomic physics and plasma modeling codes (black, Fontes et al. 2020Fontes et al. , 2023 Figure 13. Comparison of Planck's mean opacity as a function of temperature and for = 10 −13 g cm −3 and = 1 day for Nd and U computed with the FAC (red colors) and HFR (blue colors) codes and from published data from GRASP2K (black, Gaigalas et al. 2019) and HULLAC (green, Tanaka et al. 2020).