Theoretical investigation of energy levels and transitions for Ce III with applications to kilonova spectra

Doubly ionized cerium (Ce$^{2+}$) is one of the most important ions to understand the kilonova spectra. In particular, near-infrared (NIR) transitions of Ce III between the ground (5p$^6$ 4f$^2$) and first excited (5p$^6$ 4f 5d) configurations are responsible for the absorption features around 14,500 A. However, there is no dedicated theoretical studies to provide accurate transition probabilities for these transitions. We present energy levels of the ground and first excited configurations and transition data between them for Ce III. Calculations are performed using the GRASP2018 package, which is based on the multiconfiguration Dirac-Hartree-Fock and relativistic configuration interaction methods. Compared with the energy levels in the NIST database, our calculations reach the accuracy with the root-mean-square (rms) of 2732 cm$^{-1}$ or 1404 cm$^{-1}$ (excluding one highest level) for ground configuration, and rms of 618 cm$^{-1}$ for the first excited configuration. We extensively study the line strengths and find that the Babushkin gauge provide the more accurate values. By using the calculated gf values, we show that the NIR spectral features of kilonova can be explained by the Ce III lines.


INTRODUCTION
Binary neutron star (NS) mergers have been thought to be the origin of rapid neutron capture (r-process) elements.In 2017, the gravitational wave (GW) from the NS merger (GW170817) was successfully detected, and the associated electromagnetic counterpart (AT2017gfo) was observed (Abbott et al. 2017).The observed properties of AT2017gfo in UV/optical/near-infrared (NIR) wavelengths were found to be consistent with what expected as a kilonova (e.g., Kasen et al. 2017;Tanaka et al. 2017;Kawaguchi et al. 2018;Rosswog et al. 2018), thermal emission powered by radioactive decay of r-process nuclei (e.g., Li & Paczyński 1998;Metzger et al. 2010;Tanaka 2016;Metzger 2019).This fact provided us with evidence that r-process nucleosyntehsis occurs in NS mergers.
To interpret the observational properties of kilonovae, atomic data for heavy elements synthesized in NS merger ejecta are necessary.In NS merger ejecta, photons interact with matter mainly via bound-bound transitions before escaping from the system.Boundbound opacity plays a major role to determine the behavior of kilonova light curves (e.g., Kasen et al. 2013;Tanaka & Hotokezaka 2013).This requires the complete knowledge of the bound-bound transitions for all the heavy elements.Since such data are not avail-⋆ E-mail: gediminas.gaigalas@tfai.vu.lt able experimentally, theoretical atomic calculations for heavy elements have been performed to construct the atomic data (e.g., Kasen et al. 2013Kasen et al. , 2017;;Gaigalas et al. 2019Gaigalas et al. , 2020;;Gaigalas et al. 2022;Tanaka et al. 2018Tanaka et al. , 2020;;Fontes et al. 2020;Banerjee et al. 2020Banerjee et al. , 2022)).Such theoretical data have been used to evaluate the opacities in kilonova ejecta, which provides the foundation of light curve modeling of kilonovae.
On the other hand, to interpret the detailed spectral features of kilonovae, more accurate atomic data are necessary.Domoto et al. (2022) have pointed out the importance of doubly ionized Ce (Z = 58) in the spectra of kilonovae.The Ce III transitions produce distinct absorption features in the spectra at NIR wavelengths: three transitions at ∼ 16000 Å between energy levels of 5p 6 4f 2 − 5p 6 4f 5d cause the features.In fact, the features caused by Ce III nicely match with the observed features in the spectra of AT2017gfo.However, since the transition probabilities (g f -values) of the transitions are experimentally unknown, they adopted the theoretical values (Tanaka et al. 2020) whose accuracy is not certain.Domoto et al. (2023) further tested the identification of Ce III in the spectra of AT2017gfo, by estimating the g f -values of the three Ce III lines using absorption lines in stellar spectra.Such, socalled "astrophysical g f -values" broadly agree with those used in Domoto et al. (2022), which currently support the identification of Ce III in the spectra of AT2017gfo.
In fact, spectra of Ce 2+ ion have been studied by several experiments.Sugar (1965) performed observations of the spectrum of Ce III from 757 to 11 091 Å and discovered one hundred twenty-six newly energy levels, including revised values of previously known levels.Johansson & Litzén (1972) also measured the wavelengths of the 5p 6 4f 2 − 5p 6 4f 5d transition in a region between 11000 -26 000 Å.These data of the Ce 2+ ion were critically evaluated by Martin et al. (1978) and are given in the Atomic Spectra Database (ASD) of the National Institute of Standards and Technology (NIST, Kramida et al. 2024).Andersen & Sørensen (1974) measured the lifetimes for six excited levels (of the 5p 6 4f 6p configuration) using the beam-foil method.Radiative lifetimes of nine levels of the 5p 6 4f 6p configurations were measured by Li et al. (2000) using the time-resolved laser-induced fluorescence technique.They also presented transition probabilities, these were obtained from branching fractions calculated by the Cowan code (Cowan 1981) and the experimental lifetimes.
Atomic data of this ion have also been studied with semiempirical and ab-initio theoretical calculations.Mainly the ground configuration or low excited configurations were investigated.Bord et al. (1997) used the Cowan code to calculate oscillator strengths.Wyart & Palmeri (1998) investigated the energy levels and transition parameters using parametric fit, and reported ten new levels and classified more than 70 new lines.Biémont et al. (2002) studied the importance of core-polarization effects on oscillator strength in Ce III using the relativistic Hartree-Fock (HFR) method.Quinet & Biémont (2004) used the HFR method to calculate the Landé g-factors for doubly ionized lanthanides (Z=57-71).Stanek & Migdałek (2004) used the multiconfiguration Dirac-Fock method to study transition parameters for 6s 2 1 S 0 − 6s 6p 1 P 1 , 3 P 1 transitions in rare earth ionized systems (from La +1 through Nd +4 ).Li et al. (2014) used semi-empirical methods to compute M1 and E2 transition data within the ground configurations of some Ba-like and Dylike ions.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 the levels of ground and first excited configurations.Froese Fischer & Godefroid (2019) investigated the effect of electron correlation on the energy levels of the ground configuration [Xe]4f 2 using the Grasp code.Carvajal Gallego et al. (2021) used the relativistic multiconfiguration Dirac-Hartree-Fock to compute energy spectra and radiative transition data.
In this paper, we perform ab-initio calculations for Ce 2+ with the Grasp2018 (Froese Fischer et al. 2019) code by focusing on the ground (5p 6 4f 2 ) and first excited (5p 6 4f 5d) configurations.Then, toward the application to NIR spectral feature of kilonovae, we compute electric dipole (E1) transitions.This paper is organized as follows.In Section 2, we describe our computational procedures.We show our results in Section 3, by providing intensive evaluation of energy level data and transition data.In Section 4, we apply our atomic data to the kilonova spectra and discuss the impact to the NIR spectral features.Finally, we give conclusions in Section 5.

COMPUTATIONAL PROCEDURE AND SCHEME
The calculations are performed using the Grasp2018 package, which is based on the multiconfiguration Dirac-Hartree-Fock (MCDHF) and relativistic configuration interaction (RCI) methods.More details about these methods can be found in Fischer et al. (2016) and Grant (2007).
Table 1.Summary of computed eigenvalues for each J of the even and odd configurations in extended optimal level scheme.

J
Parity Eigenvalues J Parity Eigenvalues In the paper, only the main steps of the computational procedure are described.The initial wave functions were generated in the same way as in the previous computations of lanthanide ions (Gaigalas et al. 2019(Gaigalas et al. , 2020;;Radžiūtė et al. 2020Radžiūtė et al. , 2021;;Rynkun et al. 2022;Gaigalas et al. 2022).At first, the MCDHF computations of the ground [Xe]4f 2 configuration were performed.These orbitals were kept frozen and used for the next steps of the computations.Further, the 5d orbitals belonging to the first excited configuration set were optimized.
The calculations for even and odd states were performed simultaneously.In the next steps of the MCDHF computation, active spaces (AS) of CSFs were generated by allowing single-double (SD) electron substitutions from the 5s, 5p, 4f, 5d shells to the orbital spaces (OS): OS 1 = {6s, 6p, 6d, 5f, 5g}, OS 2 = {7s, 7p, 7d, 6f, 6g, 6h}, OS 3 = {8s, 8p, 8d, 7f, 7g, 7h}.When a new OS is computed, the previous orbitals are frozen.The [Kr]4d 10 defines an inactive closed core and no substitutions were allowed from it.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, including the Breit interaction and leading QED effects, were performed.
RCI calculations were performed in the extended CSFs basis.The electron correlations were extended by opening closed 3d, 4s, 4p, 4d shells step by step to find the most appropriate computational scheme for energies and transition parameters calculations.The intermediate results are described and discussed in Section 3. The scheme chosen for the final calculations include electron correlations when SD substitutions were allowed from the 4d, 5s, 5p and 4f, 5d shells to the OS 3 , single-restricted double (SrD) substitutions were allowed from the 4s, 4p shells to the OS 1 (+4d SrD4s4p scheme).The restrictions on substitutions were applied, while allowing substitutions from the deeper core shells increases the AS very rapidly.The number of CSFs in the final even and odd state expansions distributed over the different J symmetries is 8975504 and 14265384 for even and odd parity, respectively.

RESULTS
The accuracy of the computed results was evaluated comparing with data from NIST ASD (Kramida et al. 2024) and other theoretical computations.In the following sections the influence of the correlations to the energy levels and transition data was studied.The most appropriate computational scheme was chosen for both (energies and transition parameters) calculations.

Evaluation of energy spectra
The influence of the correlations was studied by opening the closed shells step by step for substitutions.To reduce the computational resources, the importance of the correlations of the closed shells was studied for the levels of the ground configurations with J = 4 (since it is the ground level) and the levels of the first excited configuration with J = 3.The contributions of the correlations are presented in Fig. 1.The figure shows that correlations from the 4d shell are important.By opening the 4p shell the energies of the ground configuration almost do not change; the energies of the first excited configuration are too high.Restricting the substitutions from the 4p and 4s shells (allowing only S substitutions; +4d SrD4s4p scheme) improves the agreement with the NIST data.When the 3d shell is opened, the agreement remains similar.The importance of correlations for transition data was also studied and described in Section 3.2.Therefore, the +4d SrD4s4p scheme was chosen to calculate the energy levels for both configurations.Fig. 2 presents the comparison of the final results (using +4d SrD4s4p scheme) with the NIST ASD.As seen in the figure, the differences for the energy levels of two configurations up to 12000 cm −1 energy reaches 600 cm −1 .The disagreement of other energies reaches 2500 cm −1 , and the largest difference (8200 cm −1 ) is for the level of the ground configuration (4f 2 1 S 0 ).The level 4f 2 1 S 0 is remote from the other levels of the ground configuration, its energy is 41076 cm −1 .Therefore, other additional correlations are relevant for this level.The root-mean-square (rms) deviations obtained for the energy levels of the ground configuration from the NIST data are 2732 cm −1 , but excluding the level with the worst disagreement (4f 2 1 S 0 ), the rms is 1404 cm −1 .The rms for the first excited configuration is 618 cm −1 .
The final results were also compared with other theoretical calculations.The comparison is presented in Fig. 3.It should be mentioned that Safronova et al. (2015) studied only some levels of the ground and first excited configurations, while Froese Fischer & Godefroid (2019)  the NIST data with excluded level (4f 2 1 S 0 ) is 1777 cm −1 .The rms for energy levels of the ground configuration by Carvajal Gallego et al. (2021) from the NIST data with excluded level (4f 2 1 S 0 ) is 1392 cm −1 , and 565 cm −1 for the first excited configuration.The levels of the first excited configurations computed by Safronova et al. (2015) disagree about 4000 cm −1 .The final results of the energy spectra along with the atomic state function composition in LS -coupling are given in Table 2.

Evaluation of transition data
The importance of correlations for transition data was also studied.
As it was mentioned above, the importance of correlation effects was studied for 4f 2 (J = 4) and 4f 5d (J = 3) levels.In Figs. 4 and 5 are presented the contributions of the electron correlation effects to the line strengths for a few transitions.From the figures it seen that there are large disagreements between the Babushkin and Coulomb forms.By analyzing the impact of electron correlation to line strength it seen that the Babushkin form is more stable than the Coulomb form for all studied transitions.It is important to include correlations from 4d orbitals, and further opening the 4p and 4s core have a smaller effect and almost do not change by opening the 3d orbitals.The same trend is observed for the remaining studied transitions, too.Therefore, the +4d SrD4s4p scheme was chosen to calculate all E1 transitions between two configurations.Recently, a new method based on the stationary second-order Rayleigh-Schrödinger many-body perturbation theory in an irreducible tensorial form (RSMBPT) has been developed to estimate various correlations.This method was already tested on light and moderate complexity ions, taking into account the core-valence (CV) and core (C) correlations (Gaigalas et al. 2024a,b) to compute the energy levels.The expressions for estimating the core-core (CC) and valence-valence (VV) correlations have been derived and the papers are in preparation (Gaigalas et al. 2024c,d).Since using the regular Grasp2018 computational scheme (+4d SrD4s4p) the agreement between obtained line strengths in the Babushkin and Coulomb gauges is quite poor, the RSMBPT method was applied for Ce III calculations.These studies were carried out only for 4f 2 (J = 4) and 4f 5d (J = 3) levels to reduce the computational resources.In these computations we use the same wavefunctions as described above.We open the core till the 3s subshell.Thus using RSMBPT method we select the most important configurations by the CV, C, CC, VV correlations impact with the specified fraction (expressed in the percentage) of the total contribu-tion (in this case 97%).Other correlations, which can not be estimated using the RSMBPT method were included in the RCI calculations in a regular way.Further the transitions were computed between these levels.The results of these computations are marked as RCI (RSMBPT) 97%.The line strengths in two gauges together with the cancellation factors (CFs), G S =0 parameter and accuracy classes according to the QQE method (Rynkun et al. 2022;Gaigalas et al. 2022) are presented in Table 3 for few transitions.As seen from Table 3, by opening deeper core (up to 3s orbital), the line strength in Babushkin gauge has changed very little.Meanwhile the results and core correlations analysis show that the Coulomb gauge is unstable and very sensitive to correlations.
As is well known, the Coulomb gauge describes better the part of the radial wavefunctions that is closer to the nucleus.We have therefore tried to describe better this part by including the core correlations in the MCDHF calculations (Gustafsson et al. 2017).The wavefunctions were obtained in following way: the first steps getting radial wavefunctions for 1s,...,5d orbitals are the same as described above.The orbitals belonging to the OS 1−3 were optimized using MCDHF in the extended optimal level EOL scheme.AS was generated by allowing SD electron substitutions from the 4f, 5d shells and S substitutions from 3s, 3p, 3d, 4s, 4p, 4d, 4s, 5p to the OS: OS 1 = {6s, 6p, 6d, 5f, 5g}, OS 2 = {7s, 7p, 7d, 6f, 6g, 6h}, OS 3 = {8s, 8p, 8d, 7f, 7g, 7h}.Further using RSMBPT method the most important configurations by the CV, C, CC, VV correlations impact with the specified fraction (expressed in the percentage) of the total contribution (in this case with 97% and 98%) was selected and performed RCI computations including correlations, which can not be estimated using the RSMBPT method.The transitions were computed between these levels and these results are marked as RCI (RSMBPT) 97% new wavefunctions and RCI (RSMBPT) 98% new wavefunctions.The results are presented in Table 3.As seen from the table, using the RSMBPT method with recalculated wavefunctions the agreement between two forms is better.It should also be noted that the values of the line strengths are close to those of the +4d SrD4s4p strategy in the Babushkin gauge.Comparing the CF in both forms using different strategies, it is seen that the CF is larger in the Babushkin gauge in all cases.A small value of the CF (less than 0.1 or 0.05 (Cowan 1981)) indicates that the calculated transition parameter is affected by strong cancellation effects.The values of CF for these given transitions are CF B > 0.205, and CF C > 0.00139.We also see that CF C changes much more than CF B using different computational schemes, indicating that the Coulomb gauge is much more sensitive to the correlations.So, after all the investigations and analysis, the Babushkin gauge should be more accurate and be closer to the exact value.Using the RSMBPT method for Ce III computations we do not get good energy differences between levels of the ground and excited configurations.Therefore, further investigations using the RSMBPT method for more complex ions are needed, as such computations are extremely demanding even using standard schemes.
Since no experimental transition data for the studied transitions of Ce 2+ are available, the obtained data are compared with other theoretical results.A comparison of the computed line strengths from the present work and the results of Carvajal Gallego et al. (2021), Biémont et al. (2002), Wyart &Palmeri (1998), andTanaka et al. (2020) (see also Domoto et al. 2022) is given in Table 4.There are disagreements between the theoretical results.Comparing the line strengths computed in this work with the results by Carvajal Gallego et al. (2021), Biémont et al. (2002), Wyart & Palmeri (1998), it is seen that they are in better agreement with the Babushkin form.Analyzing the ratios between the Babushkin and  the Coulomb gauges given by Carvajal Gallego et al. (2021), it is seen the large disagreement between two gauges in their results, too.
A comparison of the computed g f -values and other theoretical results is also shown in Fig. 6.Here, we show only the g f -values of the three transitions (whose wavelengths are shown in each panel) that are the strongest and important to interpret the spectral features of kilonovae (Domoto et al. 2022, see Section 4).The "astrophysical g f -values" (Domoto et al. 2023) are also shown in the figure.The final results in the Babushkin gauge are closer to the values from other works.
The final results (using the +4d SrD4s4p scheme) of the transition data, as wavelengths, weighted oscillator strengths, line strengths, and transition rates of the E1 transitions, are given in Table 5 along with the G S =0 parameter, and the estimated accuracy for line strengths in the Babushkin gauge according to the QQE method (Rynkun et al. 2022;Gaigalas et al. 2022).The uncertainties of the line strengths were estimated using the QQE method, since experimental data and other theoretical results with estimated error bars
are not available.The full table is available as online supplementary material.

APPLICATIONS TO KILONOVA SPECTRA
Here, we apply the final results of g f -values of the Ce III lines to calculate kilonova spectra.To calculate synthetic spectra of kilonovae, we use a wavelength-dependent radiative transfer simulation code (Tanaka & Hotokezaka 2013;Tanaka et al. 2014Tanaka et al. , 2017Tanaka et al. , 2018;;Kawaguchi et al. 2018Kawaguchi et al. , 2020)).The photon transfer is calculated by the Monte Carlo method.The setup of simulation is identical to that in Domoto et al. (2022) and Domoto et al. (2023), but we adopt the g f -values of the Ce III lines computed in this work (Table 5).For more details of the simulation, we refer the readers to Domoto et al. (2022).Among the transitions between the states of the ground configuration and the first excited configuration, it has been shown that the three transitions shown in Fig. 6 gives the strongest contribution to the NIR spectral features of kilonovae (Domoto et al. 2022).Thus, although we have all the transition probabilities of the computed transitions, we take the g f -values of these three lines to see their effects on the spectral features.As the line strengths suggest that the Babushkin gauge is more reliable, we adopt the g f -values in the Babushkin gauge in the simulations (i.e., the 5th column of Table 5).
Fig. 7 shows the comparison between the synthetic spectra (blue) and the observed spectra of AT2017gfo taken with the Very Large Telescope (VLT) at t = 1.5, 2.5, and 3.5 days after the merger (gray, Pian et al. 2017;Smartt et al. 2017).The observed spectra at t = 4.5 days after the merger taken with the Hubble Space Telescope (HST), which are not affected by telluric absorption, is also shown (black, Tanvir et al. 2017).In the synthetic spectra, the absorption features appear around 14500 Å, which are caused by the Ce III lines.These lines are blueshifted according to the velocity of the line-forming region at the NIR wavelengths (e.g., v ∼ 0.1 c at t = 2.5 days).Since our calculated g f -values are smaller than those adopted in Domoto et al. (2022), the absorption features in the new synthetic spectra (orange) become slightly weaker (see the inset of Fig. 7 for the enlarged view at t = 2.5 days).Nevertheless, the absorption features are still clearly present, and the strength of the feature is broadly consistent with the observed spectra of AT2017gfo.This gives the further support to the identification of the Ce III lines in the NIR spectra of kilonova.

SUMMARY AND CONCLUSIONS
We calculated the energy levels of the ground and first excited configurations for the Ce III using the Grasp2018 code.The energy differences between the final Grasp2018 results and the NIST ASD for two configurations up to 12000 cm −1 reach 600 cm −1 .The disagreement for other energies reaches 2500 cm −1 , and the largest difference (8200 cm −1 ) is for the level of the ground configuration (4f 2 1 S 0 ).The rms deviations obtained for the energy levels of the ground configuration from the NIST data are 2732 cm −1 , but excluding the level with the largest discrepancy (4f 2 1 S 0 ), the rms is 1404 cm −1 .The rms for the first excited configuration is 618 cm −1 .We also computed E1 transition data between the levels of the ground and first excited configurations.The uncertainties of the E1 line strengths in the Babushkin gauge are estimated based on the QQE method described in Rynkun et al. (2022); Gaigalas et al. (2022) giving the accuracy classes according to the NIST ASD (Kramida et al. 2024).The line strengths were also compared with the results of other calculations.The analysis of the line strengths shows that the Babushkin gauge should be the more accurate, so we suggest use the Babushkin gauge with the assigned accuracy for the obtained results although the line strengths are assigned with E accuracy class.
Finally, we performed radiative transfer simulations for kilonova spectra by using the calculated g f values.The synthetic spec-tra clearly show the absorption features around 14,500 Å, which is caused by the blueshifted Ce III lines.Therefore, our ab-initio atomic calculations support the identification of the Ce III lines in the NIR spectra of kilonova.

Figure 1 .
Figure1.Differences between NIST ASD energy levels and those of the present Grasp2018 calculations using different computational schemes (in cm −1 ).

Figure 2 .Figure 3 .
Figure 2. Differences between NIST ASD energy levels and those of the present Grasp2018 calculations (in cm −1 ).

Figure 6 .
Figure 6.Comparison of the computed g f -values (black) with results from Carvajal Gallego et al. (2021), Biémont et al. (2002), Wyart & Palmeri (1998), Tanaka et al. (2020), and Domoto et al. (2023) (gray).For the computed values in this work, the results using +4d SrD4s4p scheme as the final values are shown.B and C indicate the results in the Babushkin and Coulomb gauges, respectively.Note that the filled symbols show the calibrated values based on the line strengths, considering the differences in theoretical energy levels to experimental energy levels, while the open symbols show the values as the results are.

Figure 7 .
Figure7.Comparison between the synthetic spectra using the final g fvalues computed in the Babushkin gauge (blue) and the observed spectra of AT2017gfo taken with VLT (gray,Pian et al. 2017;Smartt et al. 2017) at t = 2.5 and 3.5 days after the merger as well as with HST (black,Tanvir et al. 2017) at t = 4.5 days after the merger.The orange lines are the synthetic spectra ofDomoto et al. (2022).Gray shaded areas show the regions of strong atmospheric absorption.Spectra are vertically shifted for visualization.The inset shows the enlarged view of the synthetic spectra in the NIR region at t = 2.5 days.The curves in light colors show the original results, while those in dark colors show smoothed spectra for visualization.

Table 2 .
Atomic state function composition (up to three LS components with a contribution > 0.02 of the total atomic state function) in LS -coupling and energy levels (in cm −1 ) for Ce III.Energy levels are given relative to the ground state.

Table 4 .
Comparison of computed line strengths (S in a.u.) with results from Carvajal Gallego et al. (2021); Biémont et al. (2002); Wyart & Palmeri (1998); Tanaka et al. (2020).S B is the line strength in the Babushkin gauge, S C is the line strength in the Coulomb gauge.B/C is the ratio Babushkin over Coulomb gauges.