© ESO 2022

1 Introduction

Observations of gravitational waves (GWs) from the neutron star (NS) merger GW170817 (Abbott et al. 2017a) and its electromagnetic counterpart (AT2017gfo, Abbott et al. 2017b) paved the way for a new method with which to study the origin of heavy elements in the Universe. The observed properties of AT2017gfo at ultraviolet, optical, and infrared wavelengths are found to be consistent with theoretically expected “kilonova”, that is, thermal emission powered by radioactive decay of newly synthesized r-process elements (e.g., Li & Paczyński 1998; Metzger et al. 2010). This agreement indicates that r-process synthesis takes place in the NS merger.

Atomic data for heavy elements can be used to interpret the observed properties of kilonovae. For example, the spectra of AT2017gfo show various broad features (e.g., Chornock et al. 2017; Pian et al. 2017; Smartt et al. 2017), which are imprints of heavy elements synthesized by the NS merger. However, the elemental species have not yet been identified, except for Sr (Watson et al. 2019; Domoto et al. 2021). To firmly identify the heavy elements in kilonova spectra, accurate atomic data for selected, dominant transitions are necessary.

In addition, complete atomic data are also important for understanding the nucleosynthesis in the NS merger. This is because, complete atomic data are necessary to calculate the opacity in the ejecta, which in turn determines the peak luminosity and the peak timescale of the kilonova. In particular, the presence of the lanthanide elements plays an important role in shaping the light curve of a kilonova. This is because these elements have many transitions at infrared wavelengths due to the small energy spacing, which redistribute the emission from optical to near-infrared wavelengths, making the color of the kilonova redder (the physics of the formation of emission spectra of kilonovae, in particular the role of lanthanide elements, is explained in Kasen et al. 2013; Barnes & Kasen 2013; Tanaka & Hotokezaka 2013). This property was indeed used to identify GW170817/AT2017gfo as a kilonova (e.g, Chornock et al. 2017; Drout et al. 2017; Kasliwal et al. 2017; Pian et al. 2017; Smartt et al. 2017; Tanvir et al. 2017; Utsumi et al. 2017).

To better understand the link between the emission properties of a kilonova and the nucleosynthesis therein, it is important to establish the atomic opacities of heavy elements, which requires complete atomic data (i.e., energy levels up to the excited levels and transition probabilities to the excited levels). As it is difficult to experimentally derive the complete data of energy levels and transitions, past works relied on atomic structure calculations (Kasen et al. 2013, 2017; Tanaka et al. 2018, 2020; Fontes et al. 2020; Banerjee et al. 2020). However, in such complete atomic calculations, including excited levels, the accuracy is typically low: ~ (10–20)% in the energy levels. As a result, there is a systematic uncertainty of a factor of about 2 in the opacity (Tanaka et al. 2018). To further reduce the uncertainty in the opacity, more complete and accurate atomic calculations are necessary.

In this paper, we perform ab initio atomic calculations for triply ionized Ce (Ce IV) and provide an extensive evaluation of the atomic data. We choose Ce³⁺ because it isone of the simplest lanthanide ions with a 4f electron in the ground configuration, and is therefore an ideal ion as a benchmark of atomic calculations for complicated lanthanide elements.

The atomic structure and transitions of Ce IV have been studied by several authors. Reader & Wyart (2009) measured the lines – most of which arise from inner-shell-excited configurations – and determined 112 energy levels. Zhang et al. (2001) measured the lifetimes for two excited levels (of the 5p⁶ 5d configuration) of Ce IV using the time-resolved laser-induced fluorescence technique. Savukov et al. (2003) calculated energies and transition rates for low-energy levels of Ce IV (and Pr V) using relativistic many-body perturbation theory. Safronova et al. (2015) used a configuration interaction approach with second-order perturbation theory and a linearized coupled-cluster all-order method to compute excitation energies of Ce IV along with other ions of Ce. Zilitis (2014) calculated oscillator strengths and lifetimes of the cesium isoelectronic sequence using the Dirac–Fock method. Recently, Carvajal Gallego et al. (2021) used the relativistic multiconfiguration Dirac-Hartree–Fock to compute energy spectra and radiative transition data for Ce IV, as well as Ce II and Ce III.

In the present work, we performed energy spectra computations using GRASP2018 (Froese Fischer et al. 2019) and HULLAC (Hebrew University Lawrence Livermore Atomic Code; Bar-Shalom et al. 2001) codes for the 148 energy levels of the even 5p⁶ 5d, 5p⁶ 6s, 5p⁵ 4f5d, 5p⁶ 6d, 5p⁶ 7s, 5p⁵ 4f6s, 5p⁶ 7d, 5p⁶ 8s, 5p⁶ 5g configurations and for the 77 levels of the odd 5p⁶4f, 5p⁶ 6p, 5p⁵ 4f², 5p⁶ 5f, 5p⁶ 7p configurations. Furthermore, we computed electric dipole (E1) transition data between these levels.

2 Computational procedure and scheme

2.1 GRASP2018 computations

The GRASP2018 package used for the computations is based on the multiconfiguration Dirac–Hartree–Fock (MCDHF) and relativistic configuration interaction (RCI) methods. In the MCDHF approximation, atomic state functions (ASFs) are given as linear combinations of symmetry adapted configuration state functions (CSFs), $$\Psi (\gamma PJM)={\sum }_{i=1}^{N_{\text{CSFs}}}{c}_i\Phi ({\gamma }_{i}PJM),$$(1)

where J and M are the angular quantum numbers, P is parity, N_CSFs is the number of CSFs, c_i are signed amplitudes of CSF contributions, and γ contains all additional quantum numbers required to uniquely characterize a CSF. The CSFs Φ(γ_i PJM) are built from products of one-electron Dirac orbitals. In the relativistic self-consistent field procedure, both the radial parts of the Dirac orbitals and the expansion coefficients were optimized to self-consistency.

In RCI computations, the atomic state function is expanded in CSFs, and only the expansion coefficients are determined by diagonalizing the Hamiltonian matrix. The RCI method was used to include the transverse-photon (Breit) interaction and the quantum electrodynamic (QED) corrections: vacuum polarization and self-energy. More details about MCHDF and RCI methods can be found in Froese Fischer et al. (2016) and Grant (2007).

In this work, ASFs were obtained as expansions over jj-coupled CSFs. To transform these ASFs into an LSJ-coupled CSF basis for labeling purposes, the method provided by Gaigalas et al. (2003, 2017) was used.

Generation of the initial wave functions begins from the MCDHF computation of the ground [Xe]4f configuration. These orbitals were kept frozen and used in further computations. In the next step, the 5d, 5f, 5g, 6s, 6p, 6d, 7s, 7p, 7d, 8s orbitals belonging to the configurations of the multireference (MR) set were optimized (these computations are referred to as MR below). The MR set consists of the 5p⁶5d, 5p⁶ 6s, 5p⁵ 4f5d, 5p⁶ 6d, 5p⁶ 7s, 5p⁵ 4f6s, 5p⁶ 7d, 5p⁶ 8s, 5p⁶ 5g even and 5p⁶4f, 5p⁶ 6p, 5p⁵ 4f², 5p⁶ 5f, 5p⁶ 7p odd configurations. A similar technique for generation of the initial wave functions was already applied to Nd ions (Gaigalas et al. 2019) and Er III (Gaigalas et al. 2020).

Unless stated otherwise, the inactive core used in the present calculations is [Kr]4d¹⁰5s². The even and odd states were optimized separately. In the following steps of the MCDHF computation, active spaces (AS) of CSFs were generated by allowing single (S) substitutions from the 5p shell to {8s, 7p, 7d, 5f, 5g} orbitals and single-double (SD) substitutions from the 4f, 5d, 5f, 5g, 6s, 6p, 6d, 7s, 7p, 7d, 8s shells to the orbital spaces (OS): OS₁ = {9s, 8p, 8d, 6f, 6g}, OS₂ = {10s, 9p, 9d, 7f, 7g, 6h}, and OS₃ = {11s,10p,10d,8f,8g,7h}. Thus, valence-valence (VV) electron correlations were included in the computations. When a new OS is being computed, the previous orbitals are frozen. The MCDHF calculations were performed in the extended optimal level (EOL) scheme (Dyall et al. 1989). Table 1 summarizes the calculations performed for even and odd configurations by showing their J and parity values and the ASFs that were included in the optimization process. Based on the orbitals from the MCDHF calculations, further RCI calculations were performed, including calculation of the Breit interaction and leading QED effects. By analyzing the mean radii of the orbitals, we find that the 4f orbitals are closer to the other n = 4 orbitals than either 5s or 5p. Therefore, the 5s orbital is considered to be valence as well. Further on, the 5s closed shell was opened for S substitutions to the {8s, 7p, 7d, 5f, 5g} orbitals in the RCI computations (VV2 strategy) to include additional correlation effects. In the last step, the MR set was extended to include additional important configurations based on the analysis of ASFs composition. In final computations, single-restricted-double (SrD) substitutions (allowing double (D) substitutions only from different shells) to the {8s, 7p, 7d, 5f, 5g} orbitals were added in the active space (VV3 strategy). The inclusion of electron correlation allowing substitutions from the 5s shell and the extension of the MR set increase the AS very rapidly, and so limitations on substitutions were applied.

2.2 HULLAC computations

HULLAC is a coherent set of programs that calculates the atomic structure and rates of radiative and collisional processes (Bar-Shalom et al. 2001). HULLAC uses fully relativistic orbitals to calculate the energy levels and radiative transition probabilities. The atomic relativistic states are obtained from the many-electron Hamiltonian $$H=H\text{{'}}+H_{\text{Breit}}+H_{\text{QED}},$$(2)

where H′ is the Dirac–Coulomb Hamiltonian, and H_Breit and H_QED are the Breit interaction and QED correction used as the second-order perturbation.

The single electron orbital wavefunction is obtained by solving the Dirac equation for an effective central field potential. The effective potential is obtained by solving the Poisson equation for a nuclear point charge Zδ(r) and the electron charge distribution (Bar-Shalom et al. 2001) $$\rho (r)=-4\pi r^{2}qN{[r^{l+1}\text{exp}(-\alpha r/2)]}^2,$$(3)

with the boundary condition that at r = ∞, the potential converges to (Z − q)∕r. Here, N is the normalization operator, q is the number of electrons, and α is the free parameter representing the average radius of the Slater type orbital. The free parameter α is obtained by minimization of the first-order configuration-average energies. For our calculation of Ce IV, the central potential as seen by the last electron is determined by removing one electron from the ground configuration 5p⁶ 4f (Kramida et al. 2020), which gives 5p⁶. The potential is also optimized for the ground configuration. We note that the energy levels for Ce IV are calculated for all the configurations of the MR set mentioned in Sect. 2.1, including the ground state.

The N-electron zero-order solution or the configuration state function (CSF) are derived as the anti-symmetrized product of the single-electron zero-order solution of the orbitals in the chosen coupling scheme (jj in this case). A Hamiltonian matrix, including the Breit interaction and the leading QED correction of multi-CSFs, is diagonalized as in the RCI method of GRASP2018. As a result, ASFs are obtained as linear combinations of the CSFs.

3 Results

3.1 Evaluation of energy spectra

In this work, excitation energies were computed for the 148 states of the even 5p⁶ 5d, 5p⁶ 6s, 5p⁵ 4f5d, 5p⁶ 6d, 5p⁶ 7s, 5p⁵ 4f6s, 5p⁶ 7d, 5p⁶ 8s, and 5p⁶ 5g configurations and for the 77 statesof the odd 5p⁶4f, 5p⁶ 6p, 5p⁵ 4f², 5p⁶ 5f, and 5p⁶ 7p configurations using the GRASP2018 and HULLAC codes. Table A.1 displays energy spectra and the atomic state function composition in LS-coupling from the GRASP2018 calculations for levels of these configurations. The levels are given with unique labels (Gaigalas et al. 2017). In the computed spectra there are levels, for which the assigned labels are not based on the largest contribution to the composition. These are marked in bold font in Table A.1. Energy levels were computed up to the 31 eV energy limit.

The Atomic Spectra Database (ASD) of the National Institute of Standards and Technology (NIST; Kramida et al. 2020) provides only 16 energy levels of the 5p⁶nl configurations for the Ce IV spectrum; these are based on Lang (1936). Reader & Wyart (2009) presented 112 experimental energy levels, most of them belonging to the inner-shell-excited configurations (5p⁵ 4f5d and 5p⁵ 4f6s). Energy levels given by Reader & Wyart (2009) fully agree with the values given in the NIST ASD (Kramida et al. 2020) except energies for the 5p⁶ 7d ²D_5∕2 and 5p⁶ 5g ²G_7∕2,9∕2 states. Reader & Wyart (2009) found the Lang (1936) 5g, 7d_5∕2, and 8s levels to be spurious and determined new values for the 5p⁶ 7d ²D_5∕2 state and for the states of the 5p⁶5g configuration.

To evaluate the accuracy of the results, the calculated energy levels were compared with data from the NIST ASD (Kramida et al. 2020) and experimental or other theoretical results. The comparison of computed (with GRASP2018 and HULLAC codes) energy levels and experimental (Reader & Wyart 2009) or theoretical ones (Savukov et al. 2003 and Safronova et al. 2015) with the NIST ASD, giving a difference (in cm⁻¹) ΔE = E_NIST − E_calc, is presentedin Table 2. The NIST ASD specifies one energy value (225 625 cm⁻¹) for the two unresolved levels of the 5p⁶5g configuration. When comparing this energy level with other works, the level with the lower energy was compared. In the case of Reader & Wyart (2009) and GRASP2018 results, the lower energy belongs to the level with J = 7∕2, while in the case of HULLAC, the lower level has J = 9∕2. As we can see from the table, computed energies of GRASP2018 agree well with recommended values from the NIST ASD. The root-mean-square (rms) deviations obtained for energy levels of the 5p⁶ nl configurations from the NIST data are 1270 cm⁻¹. The HULLAC energy levels in comparison with the NIST data are too high, and the rms deviations from the NIST data are 5780 cm⁻¹. When comparing the results of other theoretical investigations (Savukov et al. 2003 and Safronova et al. 2015) with the NIST ASD, the rms deviations are 2160 cm⁻¹ and 2240 cm⁻¹, respectively;but those authors only investigated low-lying levels.

It should be mentioned that the order of the levels of the 5p⁶ 6d ²D term in the GRASP2018 results disagrees with the NIST ASD. Considering the order of these states in HULLAC studies and in the calculations performed by Savukov et al. (2003), it is seen that the order agreeswith the NIST ASD, but the splitting is much smaller. In the HULLAC case, it is only 15 cm⁻¹, in Savukov et al. (2003) it is 720 cm⁻¹, while in the NIST ASD this splitting reaches 1715 cm⁻¹. Reader & Wyart (2009) show that this anomalous fine structure is the result of interaction between 5p⁶ 6d and the inner-shell excited configuration 5p⁵4f5d. Whereas there is good agreement between GRASP2018 and NIST data for most energy levels, and the order for other levels is correct, a deeper analysis was performed for the levels of the 5p⁶ 6d ²D term to clarify this disagreement. In Fig. 1, the convergence of these levels by including different correlations is presented. As seen from the figure, in the MR computations, the 5p ⁶ 6d ²D_3∕2 level is lower than 5p⁶6d ²D_5∕2. By including the correlations (VV2), the order of the levels is inverted and does not change by extending the orbital space. It was shown that the energy values converge to the NIST results by extending the orbital space from OS₁ to OS₃ (in the VV2 strategy), and the energies are very close to the NIST recommended values at OS₃. By including the additional correlations (VV3 strategy; see OS₃* in Fig. 1), the order of the states is inverted and the ²D_3∕2 level becomes the lower one. The figure also shows that the Breit and QED corrections have little influence on the energies. In Fig. 1, HULLAC results are also presented. HULLAC energy levels can be seen to be higher than energies given in the NIST ASD, and as the splitting between the fine-structure levels is very small, results for the two levels overlap.

We also analyzed the composition of these levels present our findings in Fig. 2 and Table 3. In Fig. 2 only the main contributors from different strategies are presented, while Table 3 shows the total contribution of each configuration to these two levels when different strategies are applied. After including correlations with the VV3 strategy, the weight of the 5p⁶ 6d configuration decreases in the 5p⁶6d ²D_3∕2 level, and the influence of the 5p⁵4f5d configuration increases. For the 5p⁶6d ²D_5∕2 state, the situation is vice versa. From Fig. 2 it is seen that the 5p⁶6d ²D_3∕2 state is quite pure in the HULLAC computations (about 88%), while in the final GRASP2018 results the 5p⁶6d contributionis smaller; it constitutes 58% in the composition and is close to that found by Reader & Wyart (2009). The 5p⁶6d ²D_5∕2 state is purer (about 80%) in the final GRASP2018 calculations, while in the HULLAC and Reader & Wyart (2009) investigations, it is 48% and 50%, respectively. From the analysis, it is seen that the change in the order of the two states in the final results is affected by the MR set extension and inclusion of additional correlations in the VV3 strategy. To study the accuracy of the GRASP2018 results, we also analyzed the convergence for the levels with the same orbital angular momenta (l = 3) but with higher principal quantum number (n = 7). Figure 3 presents the convergence of the 5p⁶ 7d²D term states when different correlations are included. Energy levels of the 5p⁶7d configuration in the final results agree very well with NIST data. The fine-structure levels of the 5p⁶7d ²D term have the same trends as those of the 5p⁶6d ²D term. The ²D_3∕2 and ²D_5∕2 levels are inverted when the VV2 strategy is used, and the order of the states changes when correlations are added with the VV3 strategy. The accuracy of the results and the quality of the wave functions of the 5p⁶6d ²D term are analyzed in Sect. 3.2.

Figures 4 and 5 present the comparison of computed (with the GRASP2018 and HULLAC codes) energy levels with energy levels given by Reader & Wyart (2009). The HULLAC and GRASP2018 energy levels were related by comparing the compositions in jj coupling with the data of Reader & Wyart (2009). Figures 4 and 5 show three distinct groups of energy levels. As seen from Fig. 4, the agreement with GRASP2018 for the levels of the 5p⁵4fnl configurations is poorer than it is for the 5p⁶nl configurations, the difference reaches 16 000 cm⁻¹ for some levels. The HULLAC results (see Fig. 5) are in good agreement with those of Reader & Wyart (2009) for the levels of the 5p⁵4f5d and 5p⁵4f6s configurations. The rms deviations for energy spectra of these configurations from the Reader & Wyart (2009) data are displayed in Table 4.

Fig. 1 Convergence of the 5p66d 2D3∕2 and 5p66d 2D5∕2 levels when different electron correlation effects are included. OS1–OS3 marks the results from the VV2 strategy, OS3* marks the results from the VV3 strategy.

Fig. 2 Fractional coefficients of the main contributors to the 5p66d 2D3∕2 and 5p66d 2D5∕2 levels including different electron correlation effects. OS1–OS3 marks the results from the VV2 strategy, OS3* marks the results from the VV3 strategy.

Fig. 3 Convergence of the 5p67d 2D3∕2 and 5p67d 2D5∕2 levels when different electron correlation effects are included. OS1–OS3 marks the results from the VV2 strategy, OS3* marks the results from the VV3 strategy.

Fig. 4 Differences between experimental energy levels of Reader & Wyart (2009) and those of the present GRASP2018 calculations (in cm−1). “HE” means levels with the highest calculated energies (>230 000 cm−1). All HE levelsbelong to the 5p54f5d configuration.

Fig. 5 Differences between experimental energy levels of Reader & Wyart (2009) and those of the present HULLAC calculations (in cm−1). “HE” means levels with the highest calculated energies (>230 000 cm−1). All HE levelsbelong to the 5p54f5d configuration.

3.2 Evaluation of transition data

The accuracy of the wave functions and the uncertainties of the transition parameters were evaluated by: (1) analyzing the dependence of the line strength S on the gauge parameter G (Rudzikas 2007); (2) analyzing the dependence of the cancelation factor (CF; Cowan 1981) on the gauge parameter G; (3) studying the dT parameter defined in Ekman et al. (2014), which shows the relative difference between the Babushkin and Coulomb gauges of the computed transition rates; (4) comparing the computed transition data with other available data. For this purpose transitions from different regions of energy spectra were chosen as illustrated in Fig. 6. Firstly, the accuracy of the wave functions was evaluated by studying the transition parameters as a function of G.

The dependence of the line strength S on G is parabolic, $$S(G)=a{G}^2+bG+c,$$(4)

where $$a=\frac{A}{2}{(M_B-M_C)}^2,$$(5) $$b=A\sqrt{2}{M}_C(M_B-M_C),$$(6) $$c=A{M}_C^2.$$(7)

A is a coefficient, and in the case where S(G) represents the line strength, is expressed as $$A=\frac{c^2(2J_{\text{lower}}+1)}{2\Delta E^2},$$(8)

where c is the speed of light and ΔE is transition energy. M_B and M_C are the reduced matrix elements $\langle \Phi ({\gamma }_{j}PJ)\Vert T^{(1)}\Vert \Phi ({{\gamma }^{\prime }}_k{P}^{\prime }{J}^{\prime })\rangle$ (see Eq. (3) in Gaigalas et al. 2020) in the Babushkin ($G=\sqrt{2}$) and Coulomb (G = 0) gauges, respectively. Generally, the quadratic equation has two roots. Equation (4) with the coefficients a, b, and c given by Eqs. (5), (6), and (7), has only one root given by Eq. (9). The more accurate the wave functions, the flatter (the smaller the a and the larger the G value for which S = 0 (G_{S = 0})) the gauge dependence of the line strengths (Rudzikas 2007; Gaigalas et al. 2010). The a value is smaller when the agreement between the reduced transition matrix elements of the two gauges is better. The minimum of the parabola is defined as $$G_{S=0}=\frac{\sqrt{2}{M}_C}{M_C-M_B},$$(9)

and is larger when the difference between the matrix elements of the two gauges is smaller. In the case where matrix elements fully agree in both gauges, the transition parameter is independent of G.

G_{S = 0} parameter can be expressed through the line strength: $$G_{S=0}=\frac{\text{sgn}(M_C)\sqrt{2S_C}}{\text{sgn}(M_C)\sqrt{S_C}-\text{sgn}(M_B)\sqrt{S_B}},$$(10)

where $\text{sgn}(M_B)=\frac{M_B}{\left|M_B\right|}$ and $\text{sgn}(M_C)=\frac{M_C}{\left|M_C\right|}$ shows the sign of reduced transition matrix elements in the Babushkin and Coulomb gauges, respectively. S_B and S_C are the line strengths in the Babushkin and Coulomb gauges, respectively. When the G_{S = 0} parameter isexpressed through the line strength, the sign of the reduced transition matrix element is very important. In the range − ∞ < G_{S = 0} < 0 and $\sqrt{2}$ < G_{S = 0} < ∞, matrix elements usually have the same sign in the two forms. When 0 < G_{S = 0} < $\sqrt{2}$, the matrix elements are of different sign in the two forms.

The minimum of the parabola is related to the dT parameter by the expression $$\left|G_{S=0}\right|=\left|\frac{\sqrt{2S_C}(\text{sgn}(M_C)\sqrt{S_C}+\text{sgn}(M_B)\sqrt{S_B})}{dT\times max(S_C,S_B)}\right|.$$(11)

The dependencies of the line strengths on the gauge parameter G for the studied transitions are displayed in Figs. 7–13. Here, the lines strengths from the present HULLAC computations and from other works are also plotted. As seen from the figures, the correlations included in the final computations (the VV3 strategy) are very important for all studied transitions and especially for transitions between highly excited levels or from highly excited levels to the ground configuration (Figs. 9–13). In Figs. 7–8, the gauge dependence changes monotonically upon increasing the orbital spaces or upon including additional correlations. For the studied transitions, the dependence is closer to linear when correlations included in the final computations (the VV3 strategy) are added. Also, the minimum position of the parabola (the G_{S = 0} values) – which shows the accuracy of the wave functions – for the VV3 strategy is larger than for other strategies. Thus, the wave functions of the final results should be more accurate than from the VV2 strategy. The dT parameter for the studied transitions was poorer for all choices of OS of the VV2 strategy than in the final results. The G_{S = 0} values together with the dT parameters of the VV3 strategy are presented in Table 5. The magnitude of dT is smaller for transitions that have more accurate wave functions. The minimum position of the parabola is also important for predicting how the results will change by improving the calculations. If the Babushkin and Coulomb gauges are close to the minimum position of the parabola, then improving the computations will not significantly change the results. If the Babushkin and Coulomb gauges are on the parabola, then improvement of the computations will lead to a relatively large change in the results.

In this section, we present the results of an evaluation of the quality of the wave functions for the 5p⁶6d ²D_3∕2 and 5p⁶6d ²D_5∕2 levels whichwas carried out because a disagreement was found in the order of these levels as compared with the NIST data. To this end, we analyzed the transitions in which the 5p⁶6d ²D term participates, and the dependence of the line strength S on the gauge parameter G is presented in Figs. 12–13. By analyzing the gauge dependence, we see that the wave functions of the 5p⁶6d ²D levels in the final computations (the VV3 strategy) are more accurate than they are in other strategies. From the analysis, it follows that the final results (the VV3 strategy, with the 5p⁶6d ²D_3∕2 level lower than 5p⁶6d ²D_5∕2) are more accurate than results from other strategies.

Further, to determine which gauge parameter G is more suitable, we studied the cancelation factor (CF; Cowan 1981 The dependencies of the CF on the gauge parameter G for the studied transitions are shown in Figs. 14–20. A small value of the CF, for example less than 0.1 or 0.05 (these values are given in Cowan (1981)), indicates that the calculated transition parameter is affected by strong cancelation effects. By analyzing these transitions, it was observed that transition data in the Babushkin gauge are less affected by cancelation effects than transition data in other gauges, except the gauge with G = 1, which gives even larger CFs than the Babushkin gauge in the $5{\text{s}}^{2}5{\text{p}}^{6}6\text{p}{ }^2{\text{P}}_{1/2}^{°}$ – 5s²5p⁶ 5d ²D_3∕2, 4$\text{f}5{\text{s}}^{2}5{\text{p}}^5{(}^3\text{D})5\text{d}{ }^2{\text{F}}_{5/2}$ – 4$\text{f}5{\text{s}}^{2}5{\text{p}}^6{ }^2{\text{F}}_{5/2}^{°}$, and $4\text{f}5{\text{s}}^{2}5{\text{p}}^5{(}^3\text{F})5\text{d}{ }^2{\text{F}}_{7/2}$–$4\text{f}5{\text{s}}^{2}5{\text{p}}^6{ }^2{\text{F}}_{7/2}^{°}$ transitions. By analyzing the dependence of the CF on the gauge parameter G, it is seen that the trend of the dependence for the studied transitions is similar, the peak (the largest value of the CF) in all strategies is for G = 1 or $G=\sqrt{2}$, except the $4\text{f}5{\text{s}}^{2}5{\text{p}}^5{(}^3\text{D})5\text{d}{ }^2{\text{F}}_{5/2}$–$4\text{f}5{\text{s}}^{2}5{\text{p}}^6{ }^2{\text{F}}_{5/2}^{°}$ transition (see Fig. 17). For this transition, CFs with G > 1 are larger in the OS₁ –OS₃ spaces using the VV2 strategy. The situation changes when the VV3 are included in the computations. Based on this analysis, the transition parameters were computed in the Babushkin gauge.

Table A.2 presents a comparison of the computed wavelengths and transition probabilities with results from Reader & Wyart (2009), Carvajal Gallego et al. (2021), and Savukov et al. (2003). The rms of the relative differences of the GRASP2018 results from the experimental wavelengths (Reader & Wyart 2009) is about 2%, 11%, and 19% for transitions between 5p⁶nl − 5p⁶nl, 5p⁶nl − 5p⁵4fnl, and 5p⁶nl−HE, respectively.Carvajal Gallego et al. (2021) presented transition rates for 37 lines of Ce IV. However, these authors do not give the energy spectra in their paper for this ion. Therefore, only transitions for which the energy levels could be related with the NIST energies have been compared here.

Figure 21 presents a comparison of the computed line strengths from the present work and other theoretical calculations (Carvajal Gallego et al. 2021; Savukov et al. 2003) with results from Reader & Wyart (2009). The agreement between all theoretical methods is similar, except for three transitions from 5p⁶5d to the ground configuration, for which all other calculations do not agree with the results from Savukov et al. (2003). The best agreement with results from Reader & Wyart (2009) was achieved for transitions (with S >1), where the rms deviations for all theoretical methods are about (11–15)%. It should be mentioned that all transitions presented in Fig. 21 are between levels of the 5p⁶nl configurations. Figure 22 presents a comparison of the computed line strengths from GRASP2018 with results of Reader & Wyart (2009), where transitions between different groups of configurations are distinguished. The agreement between both results and the deviations within the same line strength range differ for transitions from different groups. This is especially notable for transitions with greater line strength. Figure 23 presents a comparison of transitions of the 5p⁶nl group between the results from two GRASP2018 strategies and the results from Reader & Wyart (2009). Here, most of the results are in better agreement with the VV2 strategy, which by the analysis presented above was evaluated as less accurate than the VV3 strategy.

Usually, calculated transition rates are more accurate for greater line strengths than for weaker ones. By analyzing the agreement of the line strength of the final GRASP2018 results between both forms, it was observed that ratios of the two forms have different distributions within the same range of line strengths, depending on the configurations involved. Some of the computed GRASP2018 transitions are presented in Fig. 24, where transitions between different groups of configurations are distinguished. Deviations of line strengths for the transitions from different groups differ within the same line strength range. This is also seen in Fig. 22, where the results from the present study are compared with those of Reader & Wyart (2009). This means that not all levels are described with equal accuracy, even levels from the same configuration group. When a large number of levels is computed, and the highest levels are close to the ionization limit, it is quite difficult to achieve the same accuracy for all computed levels. The goal of this work is to compute as many levels as possible, taking into account the complexity of the task and limitation of computer resources. Therefore, 225 energy levels were computed up to the 31 eV energy limit, which is close to the ionization energy (36.906 eV). Also, all possible transitions were calculated between these levels. From the above analysis, it follows that the accuracy of the transitions with similar line strengths should likely be different. Here, transitions were plotted against the G_{S = 0} parameter, which helps to identify the most accurate transitions.

The distribution of line strength S as a function of G_{S = 0} parameter is presented in Fig. 25. From the gauge dependence analysis, which shows the quality of the wave functions (Rudzikas 2007), it follows that the most accurate results have greater G_{S = 0}. Figure 25 shows that the most accurate results are in different ranges of line strengths. For example, if the G_{S = 0} parameter is fixed at a value of 100, the values of the line strength S vary in the range of (0.0001–20) a.u. Also, the line strengths of similar magnitude have different accuracy. For example, if the range of S = (0.01–0.1) a.u. is fixed, the G_{S = 0} parameter varies in the range (0.01–6000). Therefore, if transitions with similar line strengths are evaluated together, the accuracy will be overestimated or underestimated for some of them, depending on the G_{S = 0} parameter.

Taking all of the above into account, the uncertainties of the final GRASP2018 transition data were estimated based on quantitative and qualitative evaluation of the transition parameters. Figure 26 shows the connection between quantitative (natural logarithm of the ratio of line strength S in the Babushkin form to the line strength S in the Coulomb form) and qualitative (gauge dependence analysis) evaluation of the transition data. Figure 26 can be divided into three parts: (1) the upper part of the graph (blue points), where − ∞ < G_{S = 0} < 0; (2) the lower left part of the graph (red points), where $0<G_{S=0}<\sqrt{2}$; and (3) the lower right part of the graph (green points), where $\sqrt{2}<G_{S=0}<\infty$. In the case where G_{S = 0} is in the range between the Coulomb and Babushkin gauges (the lower left part of Fig. 26), the values of the transition parameter in the two forms are on different branches of the parabola represented by Eq. (4). In this case, the agreement between both gauges is random, and cannot be used as the accuracy indicator. In two other cases, when the value of G_{S = 0} is smaller than zero or greater than $\sqrt{2}$, the values in the Coulomb and Babushkin gauges are on the same branch of the parabola of Eq. (4), the ratio between the two shows real agreement and can be used for accuracy estimation. Therefore, only transitions belonging to the ranges with − ∞ < G_{S = 0} < 0 or $\sqrt{2}<G_{S=0}<\infty$ can be evaluated based on the ratio between the two gauges. The computed transition strengths from these two groups were divided into groups according to the natural logarithm of the ratio (see Fig. 27). Transitions with similar agreement between the two gauges (or with similar G_{S = 0} parameter) were assigned to the same accuracy class. The limits of the accuracy classes were chosen according to the NIST ASD (Kramida et al. 2020) terminology (AA ≤ 1%, A+ ≤ 2%, A ≤ 3%, B+ ≤ 7%, B ≤ 10%, C+ ≤ 18%, C ≤ 25%, D+ ≤ 40%, D ≤ 50%, and E > 50%). All the remaining transitions with $0<G_{S=0}<\sqrt{2}$ were not evaluated, and were assigned to the accuracy class E.

To evaluate the total uncertainties of A and gf values, the uncertainties of the transition energy should be added. Transition energy uncertainties should be scaled by a factor of 3 for A or by factor of 1 for gf, because A is proportional to the cube of transitionenergy ΔE, while gf is proportional to the first power of ΔE.

The transition data from the GRASP2018 calculations, such as wavelengths, weighted oscillator strengths, line strengths, and transition rates of E1 transitions, are given in Table 6 along with the dT parameter, the G_{S = 0} parameter, and the estimated accuracy. The full table is available at the CDS.

The computed lifetimes for some levels are compared with experimental values and other theoretical calculations in Table 7. The GRASP2018 lifetimes of the 5p⁶6d ²D are lower than the experimental values by 18% and 23% for ²D_3∕2 and ²D_5∕2, respectively.These lifetimes were recalculated using the experimental wavelengths for the radiative transition probabilities. This way, the errors stemming from uncertainties of the transition energies were eliminated, while the uncertainties coming from calculations of the line strength were estimated using the relative difference (τ_B − τ_C)∕min(τ_B, τ_C) (Col. 1(adj.) in Table 7). Recalculated GRASP2018 lifetimes of the 5p⁶6d ²D decreased a little, and the estimated uncertainty of calculations for these levels is about 29%. The HULLAC lifetimes for these levels are too low by 33% and 37% for ²D_3∕2 and ²D_5∕2, respectively.The lifetimes calculated by Savukov et al. (2003) are almost three times lower than the experimental ones. The lifetimes of the 5p⁶6p ²P levels in GRASP2018 have an uncertainty of about 3%, and they are in good agreement with other theoretical results.

Figures 28 and 29 show the statistical evaluation of oscillator strengths and wavelengths of E1 transitions from the GRASP2018 and HULLAC computations. GRASP2018 gives more transitions in the range of log(gf) < −3, while HULLAC gives more transitions with log(gf) > −3. The largest difference between the GRASP2018 and HULLAC results is in the range of log(gf) = (−3.4 to −3.2). In the statistical distribution of the wavelengths, the difference between GRASP2018 and HULLAC results is not very large.

Fig. 6 Scheme of transitions that were investigated to evaluate the accuracy of the wave functions and transition parameters. The numbering and labeling of the levels match Table A.1.

Fig. 7 Gauge dependence of the $5{\text{s}}^{2}5{\text{p}}^{6}6\text{p}{ }^2{\text{P}}_{1/2}^{°}$–5s2 5p6 5d 2D3∕2 E1 line strength for the Ce3+ ion. The dashed lines mark the values of the line strength from other works. This transition is assigned to the accuracy class A.

Fig. 8 Gauge dependence of the 5s2 5p6 5d 2D3∕2–$4\text{f}5{\text{s}}^{2}5{\text{p}}^6{ }^2{\text{F}}_{5/2}^{°}$ E1 line strength for the Ce3+ ion. The dashed lines mark the values of the line strength from other works. This transition is assigned to the accuracy class D+.

Fig. 9 Gauge dependence of the 5s2 5p6 5g 2G9∕2–$5{\text{s}}^{2}5{\text{p}}^{6}5\text{f}{ }^2{\text{F}}_{7/2}^{°}$ E1 line strength for the Ce3+ ion. The dashed lines mark the values of the line strength from other works. This transition is assigned to the accuracy class B.

Fig. 10 Gauge dependence of the $4\text{f}5{\text{s}}^{2}5{\text{p}}^5{(}^3\text{D})5\text{d}{ }^2{\text{F}}_{5/2}$–$4\text{f}5{\text{s}}^{2}5{\text{p}}^6{ }^2{\text{F}}_{5/2}^{°}$ E1 line strength for the Ce3+ ion. The dashed lines mark the values of the line strength from other works. This transition is assigned to the accuracy class B.

Fig. 11 Gauge dependence of the $4\text{f}5{\text{s}}^{2}5{\text{p}}^5{(}^3\text{F})5\text{d}{ }^2{\text{F}}_{7/2}$–$4\text{f}5{\text{s}}^{2}5{\text{p}}^6{ }^2{\text{F}}_{7/2}^{°}$ E1 line strength for the Ce3+ ion. The dashed lines mark the values of the line strength from other works. This transition is assigned to the accuracy class D+.

Fig. 12 Gauge dependence of the 5s2 5p6 6d 2D3∕2–$5{\text{s}}^{2}5{\text{p}}^{6}6\text{p}{ }^2{\text{P}}_{3/2}^{°}$ E1 line strength for the Ce3+ ion. The dashed lines mark the values of the line strength from other works. This transition is assigned to the accuracy class A.

Fig. 13 Gauge dependence of the $5{\text{s}}^{2}5\text{p{^}6\text{}}5\text{f}{ }^2{\text{F}}_{5/2}^{°}$–5s2 5p6 6d 2D5∕2 E1 line strength for the Ce3+ ion. The dashed lines mark the values of the line strength from other works. This transition is assigned to the accuracy class B+.

Fig. 14 Dependence of the cancelation factor on gauge for the $5{\text{s}}^{2}5{\text{p}}^{6}6\text{p}{ }^2{\text{P}}_{1/2}^{°}$–5s2 5p6 5d 2D3∕2 transition for the Ce3+ ion.

Fig. 15 Dependence of the cancelation factor on gauge for the 5s2 5p6 5d 2D3∕2–$4\text{f}5{\text{s}}^{2}5{\text{p}}^6{ }^2{\text{F}}_{5/2}^{°}$ transition for the Ce3+ ion.

Fig. 16 Dependence of the cancelation factor on gauge for the 5s2 5p6 5g 2G9∕2–$5{\text{s}}^{2}5{\text{p}}^{6}5\text{f}{ }^2{\text{F}}_{7/2}^{°}$ transition for the Ce3+ ion.

Fig. 17 Dependence of the cancelation factor on gauge for the $4\text{f}5{\text{s}}^{2}5{\text{p}}^5{(}^3\text{D})5\text{d}{ }^2{\text{F}}_{5/2}$–$4\text{f}5{\text{s}}^{2}5{\text{p}}^6{ }^2{\text{F}}_{5/2}^{°}$ transition for the Ce3+ ion.

Fig. 18 Dependence of the cancelation factor on gauge for the $4\text{f}5{\text{s}}^{2}5{\text{p}}^5{(}^3\text{F})5\text{d}{ }^2{\text{F}}_{7/2}$–$4\text{f}5{\text{s}}^{2}5{\text{p}}^6{ }^2{\text{F}}_{7/2}^{°}$ transition for the Ce3+ ion.

Fig. 19 Dependence of the cancelation factor on gauge for the 5s2 5p6 6d 2D3∕2–$5{\text{s}}^{2}5{\text{p}}^{6}6\text{p}{ }^2{\text{P}}_{3/2}^{°}$ transition for the Ce3+ ion.

Fig. 20 Dependence of the cancelation factor on gauge for the $5{\text{s}}^{2}5{\text{p}}^{6}5\text{f}{ }^2{\text{F}}_{5/2}^{°}$–5s2 5p6 6d 2D5∕2 transition for the Ce3+ ion.

Fig. 21 Comparison of computed line strengths of E1 transitions in Ce IV with results from Reader & Wyart (2009), Carvajal Gallego et al. (2021), and Savukov et al. (2003).

Fig. 22 Comparison of the computed line strengths of E1 transitions in Ce IV with results from Reader & Wyart (2009). Transitions between different groups of the configurations are distinguished.

Fig. 23 Comparison of the computed line strengths of E1 transitions in Ce IV from different strategies of the GRASP2018 with results from Reader & Wyart (2009).

Fig. 24 Comparison of the computed line strengths of E1 transitions between the Babushkin and Coulomb forms in Ce IV (the VV3 strategy, the final GRASP2018 results). Transitions between different groups of the configurations are distinguished.

Fig. 25 Distribution of line strength S as a funcion of the GS = 0 parameter.

Fig. 26 Natural logarithm of the ratio of the line strength S in the Babushkin form to the line strength S in the Coulomb form, as a function of the GS = 0 parameter. The figures show the results of present transition calculations (each dot represents a certain transition). Panel (a) shows three distinct parts, each plotted on a separate graph: (1) the upper part of the graph (blue points), where − ∞ <GS = 0 < 0 (or part b); (2) the lower left part of the graph (red points), where $0<G_{S=0}<\sqrt{2}$ (or part c); and (3) the lower right part of the graph (green points), where $\sqrt{2}<G_{S=0}<\infty$ (or part d).

Fig. 27 Natural logarithm of the ratio of the line strength S in the Babushkin form to the line strength S in the Coulomb form, as a function of the GS = 0 parameter inthe range where ln(SB∕SC) reflects the real accuracy of the calculation (right part of Fig. 26 a). The limits of the accuracy classes are marked in different colors. The red lines correspond to the AA accuracy class limit, and the green ones correspond to the A+ accuracy class limit. The other limits of the accuracy classes are also marked.

Fig. 28 Statistical evaluation of oscillator strength of E1 transitions from the present GRASP2018 and HULLAC computations.

Fig. 29 Statistical evaluation of wavelengths of E1 transitions from the present GRASP2018 and HULLAC computations.

Fig. 30 Expansion opacity (left) and the Planck mean opacity (right) calculated using the atomic data from the GRASP2018 and HULLAC calculations.

4 Opacities

The atomic data calculated from GRASP2018 and HULLAC are used here to calculate the bound-bound opacity of the materialejected in a binary NS merger. Here we assume that the ejected material consists of pure Ce and calculate the opacity of Ce IV by solving the ionization balance equations (see below). Although such a case does not represent a realistic condition, it helps us to understand the impact of the accuracy of atomic calculations on the opacity.

In NS mergers, the matter is expanding at high velocity with a high velocity gradient. As shown by Karp et al. (1977), the opacity can be enhanced by the expansion of the system. We use the expansion opacity formalism derived by Eastman & Pinto (1993): $${\kappa }_{\text{exp}}(\lambda )=\frac{1}{ct\rho }{\sum }_l\frac{{\lambda }_l}{\Delta \lambda }(1-{\text{e}}^{-{\tau }_l}),$$(12)

where ρ is the density of the ejecta at time t, λ_l is the transition wavelength (l is the index of the line within the wavelength grid in the opacity calculations Δλ), and τ_l is the Sobolev optical depth at the transition wavelength. We use a linearly spaced wavelength grid of Δλ = 10 Å. Emission from the NS merger spans a wavelength range from ultraviolet to near-infrared (Metzger et al. 2010; Tanaka & Hotokezaka 2013; Kasen et al. 2013), and therefore we calculate the opacity in the wavelength range of 100−25 000 Å. Thus, our calculations have a wavelength resolution of Δ λ∕λ ~ 0.1–0.0004. We consider the density of the ejecta to be ρ = 10⁻¹³ g cm⁻³, which is the typical value at t = 1 day after the merger for a typical ejected mass of 0.01 M_⊙. The Sobolev optical depth τ_l is calculated as $${\tau }_l=\frac{\pi {\text{e}}^2}{m_{\text{e}}c}{n}_{i,j}{\lambda }_l{f}_{l}t\frac{g_l}{g_0}{\text{e}}^{-\frac{E_l}{kT}},$$(13)

where e and m_e are the charge and the mass of an electron, k is the Boltzmann constant, n_i,j is the number density of the jth ion of the ith element, E_l, g_l, g₀, and f_l are the energy, the statistical weight of the lower level and the ground state, and the oscillator strength of the transition. We note that the entire transition data set has been used for the calculation of opacity without any threshold imposed based on the transition strength. We assume local thermal equilibrium (LTE) for the calculation. The population in different energy levels is determined by the Boltzmann statistic, and the fraction of ionization for the given density is provided by solving the Saha ionization equation.

The expansion opacity is convolved with a blackbody function to calculate the Planck mean opacity as $${\kappa }_{\text{mean}}=\frac{{\displaystyle {\int }_0^{\infty }{\kappa }_{\text{exp}}}(\lambda )B_{\lambda }(T)\text{d}\lambda }{{\displaystyle {\int }_0^{\infty }{B}_{\lambda }}(T)\text{d}\lambda }.$$(14)

The expansion opacity and the Planck mean opacity calculated from the GRASP2018 and HULLAC data are shown in Fig. 30. We note that Ce reaches the peak of the 3+ ionization stage around the temperature range of T ~ (10 000–15 000) K for the given density (ρ = 10⁻¹³ g cm⁻³). Beyond T = 20 000 K, the number density of the Ce³⁺ is almost negligible. The expansion opacity at T = 15 000 K is higher at shorter wavelengths. This is because the number of transitions is high at shorter wavelengths. The opacity spectra show different features for the cases of GRASP2018 and HULLAC data, which reflect the differences in the energy level distributions and transition wavelengths in the two calculations.

However, the overall properties of the opacities are similar: the Planck mean opacities derived from the atomic data calculated with GRASP2018 and HULLAC agree within 20 %. This is because mainly the strongest transitions contribute to the opacity, and these are not significantly different in the two calculations (Fig. 28). The Planck mean opacity increases with temperature, and reaches its maximum at around T = 15 000 K, beyond which it decreases. The maximum value of the Planck mean opacity is κ_mean = 0.05 cm^2 g^-1 and 0.06 cm² g⁻¹ with the data of GRASP2018 and HULLAC, respectively.

The opacity presented in this work is much higher than found in previous studies on opacity of Ce IV (Tanaka et al. 2020; Carvajal Gallego et al. 2021). This is because the previous works used only a subset of the configurations used in this work. Therefore, the number of transitions in those previous studies was small (16 and 37 in Tanaka et al. 2020 and Carvajal Gallego et al. 2021, respectively) in comparison to this work, where the number of transitions is 5056. This incompleteness of the atomic data affects the result of opacity calculation significantly.

To explore the effect of completeness of the atomic data on opacity in more detail, we performed a convergence test on the opacity calculation using only a subset of configurations (right panel in Fig. 31). We note that only the HULLAC data are used for this test. Our results show that removal of high-energy levels, that is, the energy levels of the configurations such as 5p⁶ 8s, 5p⁶5g, and so on down to 5p⁵ 4f 6s (see the left panel of Fig. 31), does not affect the opacity (see the blue curves in the right panel of Fig. 31). However, further removal of the energy levels results in a drop in the opacity (the red curves in the right panel of Fig. 31).

Transitions from relatively low energy levels (e.g., those from 5p⁵4f 5d, 5p⁵4f²) to the high-energy levels corresponding to the configurations 5p⁶ 5f, 5p⁶7s, and 5p⁶6d are important, and hence removal of these energy levels affects the opacity. Transitions between the energy levels of the 5p⁵4f 5d and 5p⁵4f² contribute significantly to the opacity. Transitions from the low-lying levels (e.g., from 5p⁶6p) to the levels of the two configurations above them are also significant. Hence, without these energy levels, opacity is lower. Transitions between the low-lying energy levels have the most influence on the opacity. Therefore, removing all the energy levels of the configurations up to 5p⁶6p (the left panel of Fig. 31) causes almost an order of magnitude drop in opacity from the default value (calculated using the atomic data corresponding to the full set of configurations). Our convergence test highlights the importance of the complete atomic data for a realistic estimate of opacities.

Fig. 31 Effects of completeness of the atomic data on the Planck mean opacity. Left panel: Distribution of the energy levels from individual configurations as obtained from HULLAC. The number of energy levels inside 0.2 eV bins is represented as in the color scale. Right panel: Planck mean opacity as a function of temperature as calculated using the atomic data from HULLAC. The thick blue line is the opacity calculated using the full set of configurations, whereas the other lines correspond to the opacity calculated using only a subset of configurations. The opacity is unaffected by the removal of the energy levels corresponding to the configurations 5p68s, 5p65g, and so on down to 5p5 4f 6s (the blue curves). However, further removal of the energy levels (e.g., those corresponding to the configurations 5p65f, 5p66d and so on) results in a drop in the opacity (see red curves) from the default value calculated using the full atomic data (the thick blue curve).

5 Conclusions

We computed energy spectra for 225 levels and E1 transition data for the Ce³⁺ ion using the GRASP2018 and HULLAC codes. The HULLAC energy levels are higher than suggested by the NIST data: the rms deviations of the HULLAC results from the NIST data are 5780 cm⁻¹. The energies computed with GRASP2018 are in good agreement with recommended values from the NIST ASD: the rms deviations are 1270 cm⁻¹. There is discrepancy in the order of the 5p⁶6d ²D_3∕2,5∕2 levels between the final GRASP2018 results and the NIST data. By analyzing the convergence and the quality of the wave functions for these levels, we find that our final results (with the 5p⁶6d ²D_3∕2 level lower than 5p⁶6d ²D_5∕2) are more accurate than results from other strategies.

Here, we developed the gauge dependency method, which can be used to evaluate the quality of the wave functions. This qualitative evaluation of the wave functions along with the agreement between both gauges was used to estimate the uncertainties of the computed line strengths. We also show that this method can be used to identify the transitions for which the ratio between the Babushkin and Coulomb gauges is a meaningful indicator of accuracy and those transitions for which the agreement between the two gauges is a random event not related to the real accuracy of the calculation.

As shown by the analysis of dependencies of line strengths on the gauge parameter G, the correlations included in the final computations (the VV3 strategy) are very important for all studied transitions. As the gauge dependencies are closer to linear when correlations included in the final computations (the VV3 strategy) are added, the wave functions of the final results should be more accurate than from other strategies. The analysis also showed thattransition data in the Babushkin gauge are less affected by cancelation effects than data in other gauges.

Using the atomic data from GRASP2018 and HULLAC, the opacity in the ejecta from the NS merger is evaluated. The opacity spectra show different features, reflecting the difference in the energy level distributions and transition wavelengths in both data sets. However, the overall trend in the opacity is similar: the Planck mean opacities calculated with the GRASP2018 data agree with those with the HULLAC data within 20%. This means that the present atomic calculations provide a reasonable estimate of the overall opacity. The calculated opacity is significantly higher than those found in previous studies (Tanaka et al. 2020; Carvajal Gallego et al. 2021). We believe this is due to to the completeness of the atomic data corresponding to the higher number of configurations. However, we note that complete ab initio atomic calculations (up to the 31 eV energy limit for Ce IV) as presented in this paper are necessary for more elements and ions to fully understand the astrophysical opacity of lanthanide elements.

Appendix A Tables of atomic data for Ce IV.

The following Tables A.1 and A.2 list the results of the calculations discussed in the text.

References

All Tables

All Figures

	Fig. 1 Convergence of the 5p66d 2D3∕2 and 5p66d 2D5∕2 levels when different electron correlation effects are included. OS1–OS3 marks the results from the VV2 strategy, OS3* marks the results from the VV3 strategy.
In the text

	Fig. 2 Fractional coefficients of the main contributors to the 5p66d 2D3∕2 and 5p66d 2D5∕2 levels including different electron correlation effects. OS1–OS3 marks the results from the VV2 strategy, OS3* marks the results from the VV3 strategy.
In the text

	Fig. 3 Convergence of the 5p67d 2D3∕2 and 5p67d 2D5∕2 levels when different electron correlation effects are included. OS1–OS3 marks the results from the VV2 strategy, OS3* marks the results from the VV3 strategy.
In the text

	Fig. 4 Differences between experimental energy levels of Reader & Wyart (2009) and those of the present GRASP2018 calculations (in cm−1). “HE” means levels with the highest calculated energies (>230 000 cm−1). All HE levelsbelong to the 5p54f5d configuration.
In the text

	Fig. 5 Differences between experimental energy levels of Reader & Wyart (2009) and those of the present HULLAC calculations (in cm−1). “HE” means levels with the highest calculated energies (>230 000 cm−1). All HE levelsbelong to the 5p54f5d configuration.
In the text

	Fig. 6 Scheme of transitions that were investigated to evaluate the accuracy of the wave functions and transition parameters. The numbering and labeling of the levels match Table A.1.
In the text

	Fig. 7 Gauge dependence of the $5{\text{s}}^{2}5{\text{p}}^{6}6\text{p}{ }^2{\text{P}}_{1/2}^{°}$–5s2 5p6 5d 2D3∕2 E1 line strength for the Ce3+ ion. The dashed lines mark the values of the line strength from other works. This transition is assigned to the accuracy class A.
In the text

	Fig. 8 Gauge dependence of the 5s2 5p6 5d 2D3∕2–$4\text{f}5{\text{s}}^{2}5{\text{p}}^6{ }^2{\text{F}}_{5/2}^{°}$ E1 line strength for the Ce3+ ion. The dashed lines mark the values of the line strength from other works. This transition is assigned to the accuracy class D+.
In the text

	Fig. 9 Gauge dependence of the 5s2 5p6 5g 2G9∕2–$5{\text{s}}^{2}5{\text{p}}^{6}5\text{f}{ }^2{\text{F}}_{7/2}^{°}$ E1 line strength for the Ce3+ ion. The dashed lines mark the values of the line strength from other works. This transition is assigned to the accuracy class B.
In the text

	Fig. 10 Gauge dependence of the $4\text{f}5{\text{s}}^{2}5{\text{p}}^5{(}^3\text{D})5\text{d}{ }^2{\text{F}}_{5/2}$–$4\text{f}5{\text{s}}^{2}5{\text{p}}^6{ }^2{\text{F}}_{5/2}^{°}$ E1 line strength for the Ce3+ ion. The dashed lines mark the values of the line strength from other works. This transition is assigned to the accuracy class B.
In the text

	Fig. 11 Gauge dependence of the $4\text{f}5{\text{s}}^{2}5{\text{p}}^5{(}^3\text{F})5\text{d}{ }^2{\text{F}}_{7/2}$–$4\text{f}5{\text{s}}^{2}5{\text{p}}^6{ }^2{\text{F}}_{7/2}^{°}$ E1 line strength for the Ce3+ ion. The dashed lines mark the values of the line strength from other works. This transition is assigned to the accuracy class D+.
In the text

	Fig. 12 Gauge dependence of the 5s2 5p6 6d 2D3∕2–$5{\text{s}}^{2}5{\text{p}}^{6}6\text{p}{ }^2{\text{P}}_{3/2}^{°}$ E1 line strength for the Ce3+ ion. The dashed lines mark the values of the line strength from other works. This transition is assigned to the accuracy class A.
In the text

	Fig. 13 Gauge dependence of the $5{\text{s}}^{2}5\text{p{^}6\text{}}5\text{f}{ }^2{\text{F}}_{5/2}^{°}$–5s2 5p6 6d 2D5∕2 E1 line strength for the Ce3+ ion. The dashed lines mark the values of the line strength from other works. This transition is assigned to the accuracy class B+.
In the text

	Fig. 14 Dependence of the cancelation factor on gauge for the $5{\text{s}}^{2}5{\text{p}}^{6}6\text{p}{ }^2{\text{P}}_{1/2}^{°}$–5s2 5p6 5d 2D3∕2 transition for the Ce3+ ion.
In the text

	Fig. 15 Dependence of the cancelation factor on gauge for the 5s2 5p6 5d 2D3∕2–$4\text{f}5{\text{s}}^{2}5{\text{p}}^6{ }^2{\text{F}}_{5/2}^{°}$ transition for the Ce3+ ion.
In the text

	Fig. 16 Dependence of the cancelation factor on gauge for the 5s2 5p6 5g 2G9∕2–$5{\text{s}}^{2}5{\text{p}}^{6}5\text{f}{ }^2{\text{F}}_{7/2}^{°}$ transition for the Ce3+ ion.
In the text

	Fig. 17 Dependence of the cancelation factor on gauge for the $4\text{f}5{\text{s}}^{2}5{\text{p}}^5{(}^3\text{D})5\text{d}{ }^2{\text{F}}_{5/2}$–$4\text{f}5{\text{s}}^{2}5{\text{p}}^6{ }^2{\text{F}}_{5/2}^{°}$ transition for the Ce3+ ion.
In the text

	Fig. 18 Dependence of the cancelation factor on gauge for the $4\text{f}5{\text{s}}^{2}5{\text{p}}^5{(}^3\text{F})5\text{d}{ }^2{\text{F}}_{7/2}$–$4\text{f}5{\text{s}}^{2}5{\text{p}}^6{ }^2{\text{F}}_{7/2}^{°}$ transition for the Ce3+ ion.
In the text

	Fig. 19 Dependence of the cancelation factor on gauge for the 5s2 5p6 6d 2D3∕2–$5{\text{s}}^{2}5{\text{p}}^{6}6\text{p}{ }^2{\text{P}}_{3/2}^{°}$ transition for the Ce3+ ion.
In the text

	Fig. 20 Dependence of the cancelation factor on gauge for the $5{\text{s}}^{2}5{\text{p}}^{6}5\text{f}{ }^2{\text{F}}_{5/2}^{°}$–5s2 5p6 6d 2D5∕2 transition for the Ce3+ ion.
In the text

	Fig. 21 Comparison of computed line strengths of E1 transitions in Ce IV with results from Reader & Wyart (2009), Carvajal Gallego et al. (2021), and Savukov et al. (2003).
In the text

	Fig. 22 Comparison of the computed line strengths of E1 transitions in Ce IV with results from Reader & Wyart (2009). Transitions between different groups of the configurations are distinguished.
In the text

	Fig. 23 Comparison of the computed line strengths of E1 transitions in Ce IV from different strategies of the GRASP2018 with results from Reader & Wyart (2009).
In the text

	Fig. 24 Comparison of the computed line strengths of E1 transitions between the Babushkin and Coulomb forms in Ce IV (the VV3 strategy, the final GRASP2018 results). Transitions between different groups of the configurations are distinguished.
In the text

	Fig. 25 Distribution of line strength S as a funcion of the GS = 0 parameter.
In the text

Fig. 26 Natural logarithm of the ratio of the line strength S in the Babushkin form to the line strength S in the Coulomb form, as a function of the GS = 0 parameter. The figures show the results of present transition calculations (each dot represents a certain transition). Panel (a) shows three distinct parts, each plotted on a separate graph: (1) the upper part of the graph (blue points), where − ∞ <GS = 0 < 0 (or part b); (2) the lower left part of the graph (red points), where $0<G_{S=0}<\sqrt{2}$ (or part c); and (3) the lower right part of the graph (green points), where $\sqrt{2}<G_{S=0}<\infty$ (or part d).

In the text

Fig. 27 Natural logarithm of the ratio of the line strength S in the Babushkin form to the line strength S in the Coulomb form, as a function of the GS = 0 parameter inthe range where ln(SB∕SC) reflects the real accuracy of the calculation (right part of Fig. 26 a). The limits of the accuracy classes are marked in different colors. The red lines correspond to the AA accuracy class limit, and the green ones correspond to the A+ accuracy class limit. The other limits of the accuracy classes are also marked.

In the text

	Fig. 28 Statistical evaluation of oscillator strength of E1 transitions from the present GRASP2018 and HULLAC computations.
In the text

	Fig. 29 Statistical evaluation of wavelengths of E1 transitions from the present GRASP2018 and HULLAC computations.
In the text

	Fig. 30 Expansion opacity (left) and the Planck mean opacity (right) calculated using the atomic data from the GRASP2018 and HULLAC calculations.
In the text

Fig. 31 Effects of completeness of the atomic data on the Planck mean opacity. Left panel: Distribution of the energy levels from individual configurations as obtained from HULLAC. The number of energy levels inside 0.2 eV bins is represented as in the color scale. Right panel: Planck mean opacity as a function of temperature as calculated using the atomic data from HULLAC. The thick blue line is the opacity calculated using the full set of configurations, whereas the other lines correspond to the opacity calculated using only a subset of configurations. The opacity is unaffected by the removal of the energy levels corresponding to the configurations 5p68s, 5p65g, and so on down to 5p5 4f 6s (the blue curves). However, further removal of the energy levels (e.g., those corresponding to the configurations 5p65f, 5p66d and so on) results in a drop in the opacity (see red curves) from the default value calculated using the full atomic data (the thick blue curve).

In the text