*Article* **Structure Calculations in Nd III and U III Relevant for Kilonovae Modelling**

**Ricardo F. Silva 1,2,\*, Jorge M. Sampaio 1,2, Pedro Amaro 3, Andreas Flörs 4, Gabriel Martínez-Pinedo 4,5,6 and José P. Marques 1,2**


**Abstract:** The detection of gravitational waves and electromagnetic signals from the neutron star merger GW170817 has provided evidence that these astrophysical events are sites where the *r*-process nucleosynthesis operates. The electromagnetic signal, commonly known as kilonova, is powered by the radioactive decay of freshly synthesized nuclei. However, its luminosity, colour and spectra depend on the atomic opacities of the produced elements. In particular, opacities of lanthanides and actinides elements, due to their large density of bound–bound transitions, are fundamental. The current work focuses on atomic structure calculations for lanthanide and actinide ions, which are important in kilonovae modelling of ejecta spectra. Calculations for Nd III and U III, two representative rare-earth ions, were achieved. Our aim is to provide valuable insights for future opacity calculations for all heavy elements. We noticed that the opacity of U III is about an order of magnitude greater than the opacity of Nd III due to a higher density of levels in the case of the actinide.

**Keywords:** opacity; atomic data; kilonovae; oscillator strengths; neutron stars

#### **1. Introduction**

The production mechanisms of elements heavier than iron have been studied for many decades. They involve a sequence of neutron captures and beta-decays. Depending on the neutron density reached in the astrophysical environment, one distinguishes between the *s*-process (*s* for slow) and *r*-process (*r* for rapid). While it has been known for a long time that the *s*-process operates on the asymptotic giant branch (AGB) stars, only recently have we been able to identify one of the astrophysical sites where the *r*-process operates [1]. This has been made possible by the observation of a kilonova associated with the collision of two neutron stars. This observation took place in August 2017 after the detection of gravitational waves from a neutron star merger by the LIGO-Virgo experiment, the wellknown GW170817 event [2]. Its electromagnetic counterpart, designated by AT2017gfo [3], exhibits a number of unique characteristics that set it apart from other transients, including an unusually high optical brightness in the days following the explosion and a longlived infrared emission that lasted nearly two weeks. These characteristics, which are

**Citation:** Silva, R.F.; Sampaio, J.M.; Amaro, P.; Flörs, A.; Martínez-Pinedo, G.; Marques, J.P. Structure Calculations in Nd III and U III Relevant for Kilonovae Modelling. *Atoms* **2022**, *10*, 18. https://doi.org/ 10.3390/atoms10010018

Academic Editors: Mustapha Laatiaoui and Sebastian Raeder

Received: 30 December 2021 Accepted: 2 February 2022 Published: 7 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

associated with the ejecta's rapid colour evolution, appear to be consistent with theoretical models of kilonovae, which consider them to be potential sites for the occurrence of heavy *r*-processes [4].

Numerous explanations have been put forward to account for the optical and nearinfrared spectral features observed in this so-called kilonova. Due to a lack of information about the atomic properties of lanthanide and actinide ions, the majority of radiation transport simulations are based on grey opacity schemes in which the value of the opacity is adjusted to reproduce the colour of the emission. It is expected, however, that the opacities of the former ions are roughly 10 times higher than the ones associated with ironlike elements [5,6]. Regarding this lack of data, calculations for selected *r*-process elements [7–11] and for all lanthanide elements [12–14] have been published in recent years.

In this work, atomic structure and opacity calculations were carried out for two representative *r*-process ions, Nd III and U III. The expansion opacity is compared between these two ions in order to evaluate the possible impact of actinides. This is particularly significant given the scarcity of publications incorporating actinides into their opacity models. With this in mind, we discuss the effect of level density in the computation of expansion opacities and in particular how actinides can be prone to having a higher density of low-lying levels.

#### **2. Methods**

#### *2.1. Expansion Opacities*

Pinto & Eastman's previous work on the light curves of type Ia supernovae [15] have shown bound–bound transitions to be the major source of opacity, accounting for two orders of magnitude more than contributions from electron scattering, bound–free, and free–free transitions.

In rapidly evolving environments, such as those of SN Ia and of kilonovae, the velocity gradient of the expansion of the ejecta is much greater than thermal and turbulent motions responsible for Doppler broadening. When that is the case, and by also assuming the impact of the overlap of strong lines is negligible (as demonstrated by Tanaka et al. 2020 for the case of lanthanides [14]), we can use the expansion opacity formalism to compute the bound–bound opacity for a certain wavelength grid Δ*λ* [16,17]. Assuming homologously expanding ejecta with density *ρ*, the expansion opacity at the time *texp* after the explosion is given as

$$\kappa\_{\exp}(\lambda\_k) = \frac{1}{\rho c t\_{\exp}} \sum\_{k} \frac{\lambda\_k}{\Delta \lambda} \left(1 - e^{-\tau\_k} \right), \tag{1}$$

where the sum is taken over all the lines within a wavelength interval of width Δ*λ*. *τ<sup>k</sup>* is the Sobolev's optical depth

$$
\pi\_k = \left(\frac{\pi e^2}{m\_\text{c} c}\right) n\_k t\_{\text{exp}} f\_k \lambda\_{k\text{s}} \tag{2}
$$

where *n* is the population of the lower level of the transition *k* with wavelength *λ<sup>k</sup>* and *fk* is the oscillator strength of the line.

#### *2.2. Atomic Calculations*

Most of the calculations for this work were performed using the open source and freely available Flexible Atomic Code (FAC) [18]. In particular, the forked version of the original code, cFAC (version 1.6.3a), was used [19]. The atomic structure fully takes into account relativistic effect as it is based on the diagonalization of the Dirac-Coulomb Hamiltonian:

$$H\_{\rm DC} = \sum\_{i=1}^{N} \left( c\boldsymbol{\alpha}\_{i} \cdot \boldsymbol{\mathcal{p}}\_{i} + (\beta\_{i} - 1)\boldsymbol{c}^{2} + V\_{i} \right) + \sum\_{i$$

where *α<sup>i</sup>* and *β<sup>i</sup>* are the 4 × 4 Dirac matrices and *Vi* accounts for potential due to the nuclear charge. Higher order QED effects such as self-energy and vacuum polarization effects are included in the screened hydrogenic approximation, while the Breit interaction is only included in the zero energy limit for the exchanged photon. Relativistic configuration interaction (RCI) calculations are performed based on a set of atomic state functions (ASFs), Ψ, given by a superposition of *i* = 1, ··· , *N*CSF configuration state functions (CSFs), *ψi*, with the same symmetries,

$$\Psi = \sum\_{i}^{\text{Nens}} c\_i \psi^i. \tag{4}$$

The CSFs used in the calculation are constructed from a linear combination of Slater determinants of one-electron orbitals, which are determined by solving self-consistently a set of Dirac–Fock–Slater differential equations for a local central potential which includes both nuclear field and electron–electron interactions. In this approach, the central potential is derived using a self-consistent method that minimizes the energy of a chosen mean configuration. In all the calculations provided with FAC, and as it is typically considered, the central potential was optimized for the ground configuration, [Xe] 4 *f* <sup>4</sup> in the case of Nd III and [Rn] 5 *f* <sup>4</sup> for U III.

#### **3. Results and Discussion**

*3.1. Nd III*

Starting with Nd III, a summary of the configurations used in each of the calculations, including the total number of levels and lines computed, is given in Table 1. For all calculations, a [Xe] ground configuration was used for the core of the ion. We performed two different FAC calculations (labelled accordingly as Calculations A and B) which are distinguished only by the inclusion in the basis set of the additional configuration 4 *f* <sup>3</sup> 6 *f* in Calculation A.

**Table 1.** Summary of the different set of configurations used on different calculations for Nd III, including experimental results from NIST [20]. For the GRASP2K calculations [21], only the configurations from the multireference space are shown. The full active space of configurations used in that calculation is shown in Ref. [22].


We compare our results with structure calculations performed with the GRASP2K code by Gaigalas et al. [22] as well with the NIST [20] database. As expected, however, the number of levels and lines measured experimentally is still very reduced.

For the data supplied by Gaigalas et al., MCDHF calculations were initially performed for the states of the ground configuration. The wave functions derived from these calculations were adopted as the initial ones to calculate the states of the multireference configurations (given in Table 1). Subsequent RCI calculations account for a larger number of configurations not included in the initial MCDHF self-consistent field. The Breit interaction and leading order QED effects are also included at this stage. The full details of the calculation and the construction of the active space used are given in [22]. As a result of the

increased degree of optimization, MCDHF-RCI calculations are expected to be more precise than RCI calculations based on a central-field potential. Nonetheless, the latter has been used in calculations for lanthanides, with a typical precision for the lowest energy states of each configuration of 20% for neutral atoms and of less than 10% for the first ionization stages (see, for example, [14]).

A visual comparison of our results with those of Gaigalas et al. and with the NIST data is shown in Figure 1. We can clearly see that the FAC results for Calculation B appear to match reasonably well to the ones obtained with the GRASP2K code, especially for the lowest lying configurations. A worse match is achieved for the 4 *f* <sup>2</sup> 5*d*<sup>2</sup> and 4 *f* <sup>2</sup> 5*d* 6*d*, as FAC seems to underestimate the energy values of those configurations. In order to compare our results, we have evaluated the difference in the lowest energy level of each configuration (normalized by the ionization potential) between our FAC calculations and the data from Gaigalas et al. and NIST. Agreement for Calculation A with the data from Gaigalas et al. and NIST is, on average, within 4.1% and 3.9%, respectively. Calculation B has an average agreement of 3.6% with Gaigalas et al. data and of 0.6% with NIST data. However, we also note that only data for the ground state and for 4 *f* <sup>3</sup> 5*d* are available in the NIST database. No other experimental data were found for Nd III configurations beyond 4 *f* 35*d*.

In our particular case, the disparity of results between Calculations A and B of FAC points to the fact that the energy levels have not yet converged. Hence, although RCI codes can provide time-efficient calculations when only few configurations are included, which allow for systematic calculations of many ions (as achieved in [14]), one should note that the inclusion of a much higher number of configurations can have a considerable impact in the calculations.

**Figure 1.** Energy levels for configurations of Nd III calculated with FAC and compared with the results from Gaigalas et al. [22] and with NIST data [20].

#### *3.2. U III*

As for the U III, the basis set for the FAC RCI calculations was determined by increasing the principal quantum number of the configurations used in Nd III by one. The set of configurations from Calculation B of Nd III was chosen because it produced the best results when compared to the experimental NIST data. At the time of writing, the only experimental data available for actinide ions are from Blaise et al. [23]. Furthermore, two independent groups from the *Los Alamos* institute have developed calculations for the

first ionized states of uranium with the CI-MBPT code [24] and the ATOMIC codes [25]. Only the data from Savukov et al. [26], which used a hybrid configuration interaction (CI) plus linearized coupled cluster (LCC) methods described in [27], are publicly available. Sultana N. Nahar of the Ohio State University's Astronomy Department hosts and maintains the NORAD-Atomic-Data database [28], which contains calculations for a broad range of structure calculations for important astrophysical ions, including uranium, using the SUPERSTRUCTURE [29] algorithm. In that database, however, only highly ionized elements are available.

Table 2 provides an overview of the FAC calculation achieved in the work as well as the prior calculations for U III with the CI-MBPT code from Savukov et al. and of the experimental data from Blaise et al.

**Table 2.** Summary of the different set of configurations used on the FAC calculation for U III. An overview of the experimental data from Blaise et al. [23]. The calculations produced by Savukov et al. [26] using a hybrid configuration interaction (CI) plus linearized coupled cluster (LCC) methods (described in [27]) are also shown. For the CI+LCC calculations, we only show the configurations for which levels and lines data were published.


<sup>a</sup> Only configurations for published levels and lines are shown. <sup>b</sup> Only the energies for 96 levels are published. <sup>c</sup> Only the *g f*-values for 20 lines are published.

We can observe from Figure 2 that, as with Nd III, we were able to reproduce the lowest lying levels fairly accurately when compared to both experimental and computational data. From spectroscopic studies, Blaise et al. determined 5 *f* <sup>4</sup> to be the ground configuration for U III. On the other hand, the CI+LCC calculations suggest an electron in the 6*d* shell in the ion's lowest energy state. Our FAC calculations do provide better agreement with the experimental data in this regard, giving the even 5 *f* <sup>4</sup> configuration as the ground state of U III. Comparisons for the lowest energies levels of each configuration were evaluated in the same way as for NdIII. We found relative differences of 3.4%, on average, when comparing our FAC calculation with the experimental data from Blaise et al. The agreement for the lowest energy levels with calculations from Savukov et al. was within 2.2%, on average.

One important point to keep is the larger radius of the 5 *f* shell when compared to the 4 *f* shell [30]. As a result, 5 *f* electrons tend to be less deeply buried in the core and less shielded from the effect of outer valence electrons than 4 *f* electrons. This effect, associated with the higher *Z* of actinide elements can, in theory, contribute to a smaller gap between the ground state and the first few excited levels. In the particular case of U III, it is expected that the excitation energies of 5 *f* and 6*d* are exceedingly close to each other.

The greater radius of 5 *f* shells also contributes to a significantly higher level density in U III when compared to the previously discussed calculations for Nd III, as can be seen in more detail in Figure 3, despite the fact that the same number of levels was calculated in both cases. The number of levels is particularly high in the case of the actinide at energies below 10 eV. This is especially important to consider in the opacity calculations as the population of low-lying levels is favoured in LTE conditions. The gap between the two peaks that we find at roughly 12 eV is likely due to the limited set of configurations used in the calculation, as it was based from the calculation on Nd III. Therefore, contributions from 5*g* shells, for example, are not included, and can potentially contribute with energy levels at those energies.

**Figure 2.** Energy levels for configurations of U III calculated with FAC. The excitation energies are compared with the experimental results from Blaise et al. [23] and the CI+LCC calculations from Savukov et al. [26].

**Figure 3.** Level density distribution of U III and Nd III as calculated from FAC. The configurations with greater contribution for each peak are also highlighted.

#### *3.3. Expansion Opacities*

Following previous works, the expansion opacity of Nd III and U III was calculated individually for each ion, assuming a gas composition given just by those ions and thermal population of excited states. This allows us to investigate the effect of only the computed bound–bound transitions on the opacity. They were evaluated over a time period of one day following the explosion, when the medium density was about *ρ* = 10−<sup>13</sup> g cm−<sup>3</sup> [12]. Additionally, a temperature of *T* = 10,000 K was specified based on prior estimates for doubly ionized ions [8,14]. The results are presented in Figure 4, including the expansion opacity calculated with the data provided by Gaigalas et al. [22] for Nd III (Strategy C). A good agreement was found between the opacity calculations using the data from Calculation B of Nd III and the data from Gaigalas et al., particularly at visible and infrared wavelengths (average deviations on the expansion opacity of 13% for *λ* > 5000 Å). Differences at lower energies are due to our opacity calculations taking into account a greater number of high energy levels than Gaigalas et al. The spectrum obtained with Calculation A for Nd III is overall similar to the one presented here for Calculation B (both in shape and magnitude), with only minor differences at very high energies (∼1000 Å) due to the presence of more highly energetic levels.

**Figure 4.** Expansion opacity for U III and Nd III (Calculation B) from the calculations performed with FAC. The expansion opacity of Nd III calculated using the data from Gaigalas et al. [22] is also to provide a more visual comparison between the two results. The opacity was evaluated at *t* = 1 day for the typical density and temperature values of *ρ* = 10−<sup>13</sup> g cm−3and *T* = 10,000 K assumed by other calculations. A wavelength bin of Δ*λ* = 10 Å was used. An average over a set of 10 bins is also shown in darker blue and orange lines. The black dashed line highlights the wavelength dependence of the opacity following approximately a *λ*−<sup>1</sup> power law.

As it can be seen, the opacity of U III is nearly an order of magnitude greater than that of Nd III. This results can be explained by the higher density of low-lying levels in the case of U III, which have a greater contribution to the opacity in LTE than more excited ones. The difference is most noticeable in the visible range, and while it remains significant in the IR, it is less noticeable due to large fluctuations in the opacity of Nd III.

Another interesting observation is that the number of lines of both actinide and lanthanide elements seems to vary smoothly with the wavelength. In particular, we found that after the initial peak at *λ* ∼ 1000 Å, the number of transitions decreases smoothly with *<sup>N</sup>* ∼ *<sup>λ</sup>*−2, especially at infrared wavelengths. This dependence is based on the low-energy tail seen empirically in the distribution of bound–bound transitions in Nd III and U III, that can be seen in Figure 5.

A parameterization of the opacities can, therefore, be obtained, when there is a high number of transitions with a high optical depth (1 <sup>−</sup> *<sup>e</sup>*−*<sup>τ</sup>* <sup>∼</sup> 1). In that case

$$
\kappa\_{\rm exp} \approx \frac{1}{\rho c t\_{\rm exp}} \sum\_{k} \frac{\lambda\_k}{\Delta \lambda} \approx \frac{1}{\rho c t\_{\rm exp}} N \frac{\lambda}{\Delta \lambda} \sim \lambda^{-1}.\tag{5}
$$

Equation 5 emphasizes the importance of the number of lines taken into account in these calculations. Given a high optical depth and LTE conditions, the primary source of opacity is the ion's line count per unit wavelength, rather than individual strong transitions. As a result, the precision of individual lines will be negligible in environments where the density of levels is sufficiently high enough to sustain local thermodynamical equilibrium and a high optical depth. These insights are particularly pertinent in the case of lanthanides and actinides, owing to their extremely complicated shell structure, which makes accurate computation extremely difficult to achieve within a reasonable amount of time and computational resources.

**Figure 5.** Total number of lines for Nd III and U III as a function of the wavelength. The *<sup>N</sup>* <sup>∼</sup> *<sup>λ</sup>*−<sup>2</sup> inverse power law dependence on the wavelength is highlighted. A wavelength bin of Δ*λ* = 10 Å was used.

#### **4. Conclusions**

In this work, the FAC code has been used to compute level energies, transition wavelengths and oscillator strengths for electric dipole (E1) transitions for Nd III and U III. We have noticed a reasonable agreement of the energies of low-lying levels with experimental data as well as with other theoretical calculations achieved for those ions. We predict the opacity of U III to be roughly an order of magnitude higher than of Nd III, in all wavelength ranges of interest. Moreover, the larger number of strong transitions of uranium makes the spectra extremely dependent on the number of transitions included, making the precision of individual lines negligible. Due to the higher density of levels expected for actinide elements when compared to lanthanides, we predict these differences to extend beyond Nd III and U III.

**Author Contributions:** Formal analysis, R.F.S., P.A. and A.F.; Project administration, G.M.-P. and J.P.M.; Supervision, J.M.S., G.M.-P. and J.P.M.; Validation, J.M.S., P.A., A.F. and G.M.-P.; Visualization, R.F.S.; Writing—original draft, R.F.S; Writing—review & editing, R.F.S., J.M.S., P.A., A.F., G.M.-P. and J.P.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was partially supported by the Fundação para a Ciência e Tecnologia (FCT), Portugal through the contract UIDP/50 0 07/2020 (LIP). R.F.S acknowledges the support from the ChETEC COST Action( CA16117) during his short-term scientific mission at GSI. A.F. and G.M.-P. acknowledge the support of the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (ERC Advanced Grant KILONOVA No. 885281).

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

