Electron-Ion Recombination in Lanthanides Part I. Radiative Recombination Rate Coefficients in Doubly Ionised Lanthanides

1. Introduction

The origin of the elements remains an open question in physics and (nuclear) astrophysics. Elements up until iron are synthetised mainly in stellar burning, being energetically unfavourable to go beyond. It is thought that a large amount of the heavier-than-iron nuclei are produced via the rapid neutron-capture process (r-process). The identification of the astrophysical environments capable of sustaining the continuous neutron fluxes required for the r-process has therefore been a major research objective for several decades [1,2,3,4].

Since the 1970s, neutron star mergers (NSMs) have been theorised as possible astrophysical sites for heavy element nucleosynthesis [1,2], mainly because the elevated temperatures and high fraction of neutrons would be indicative that heavy nuclei were being synthesised there via the r-process. The renowned 2017 detection, by the LIGO and VIRGO collaborations [5], of the gravitational wave event GW170817 [6,7] and respective electromagnetic counterpart, the astronomical transient AT2017gfo [8,9], also known as kilonova, now makes it possible to probe which elements are produced in such extreme neutron-rich events. On the other hand, the detected kilonova’s spectrum exhibits spectral features on the (near-)infrared band [10,11], compatible with the atomic structure of elements with open f-shells, such as lanthanides and actinides [3,4,12]. This makes the observation of (near-)infrared spectral lines a key diagnostic of r-process ejecta. Furthermore, elements with open f-shells have very dense sets of atomic energy levels, leading to an even denser line spectra, consequently yielding large opacities [13,14,15,16,17].

To confirm the production of heavy elements, such as lanthanides and actinides, in NSMs it is necessary to identify spectral lines (or spectral features) in the kilonova spectrum [18], being done through Collisional-Radiative Modelling (CRM). However, to perform such identification, it is necessary to know a wide array of atomic parameters for the elements to be identified; for the lanthanides and actinides, there is a lack of atomic data for this identification [19,20]. Obtaining these atomic parameters for lanthanides and actinides proves to be challenging, both experimentally and theoretically. Experimentally, lanthanides are rare (lanthanides are commonly known as rare earth elements) and expensive; Theoretically, their open f-shells configurations are computationally expensive and demanding regarding memory and computation time.

The ejecta in the AT2017gfo kilonova spectrum underwent rapid temporal evolution [11], with different atomic processes becoming dominant at different times after the collision. This behaviour resulted from the rapid decline in density and temperature, which progressively altered the balance between radiative and collisional processes. Therefore, different atomic parameters are needed to model kilonovae at different times. In the early times, in photospheric phases, local thermodynamic equilibrium (LTE) conditions are a valid assumption [19,21]: here, needed atomic data consist on atomic energy levels and allowed transition rates, with collisional processes (such as electron-impact excitation and electron-ion recombination [22,23,24]) being statistically approximated, since collisional processes dominate over radiative ones. For the later epochs, in nebular phases, LTE breaks down due to the decrease in temperature and density, entering non-LTE conditions [20,25,26,27]. In the non-LTE regime, collisional processes need to be explicitly considered (since collisional processes no longer dominate over radiative ones), as well as forbidden transitions rates: the previously needed atomic energy levels and allowed transitions rates are also needed. Efforts in modelling kilonovae lead to better estimates on the abundances of the elements produced in these astrophysical events [19,20,21], contributing to answer the question of where are heavy elements produced in the Universe.

Among the atomic processes relevant under non-LTE conditions, electron-ion recombination plays a key role in determining the ionisation balance and, consequently, the emergent spectrum of the ejecta. This work presents explicit calculations of radiative recombination, for doubly ionised lanthanides, to be used in CRM. Radiative recombination is the direct/background process in electron-ion recombination; the resonant process, dielectronic recombination, is also important to be considered in CRM, but will not be covered in this work, being covered in a future Part II. The calculations were carried out for doubly ionised lanthanides using atomic codes, considering radiative and dielectronic recombination as separated processes, under the independent treatment of radiative recombination and dielectronic recombination methodology [28,29,30].

2. Radiative Recombination

This work focuses on electron-ion recombination processes, particularly on direct radiative recombination: in short, radiative recombination (RR). This is the process in which a free electron, from the continuum, is captured into a bound state of an atom/ion (hereafter, consider “ion”, as it is a more general term). As the free electron is captured, the ion’s ionisation, q, decreases, and the electronic configuration rearranges itself, which leads to the emission of radiation, further ensuring energy conservation (see equation (1)).

$${\mathrm{e}}^{-}+{\left[{\mathrm{X}}^{q+}\right]}_i⟶{\left[{\mathrm{X}}^{(q-1)+}\right]}_f+\gamma .$$

(1)

The ion that captures the electron is commonly referred to as the recombining ion, and the ion with the newly captured electron is the recombined ion; Figure 1 shows a schematic of the process described.

To compute this process, atomic codes are used. In this work, we carried out RR calculations with the Flexible Atomic Code (FAC) [31], benchmarking simple cases with Autostructure (AS) [32]. Both codes calculate atomic processes through a mean-field approach by means of local central potentials, which are computationally advantageous (in memory usage and computing time) when compared with ab initio atomic codes. Furthermore, following up on recent studies [33,34,35], we compared the resulting RR rate coefficients with the approximation presented by T.S. Axelrod (1980) [36]. Last but not least, we fitted the RR rate coefficient curves to the formula of Verner and Ferland (1996) [37].

It is important to note that, despite the explicit calculation of RR rate coefficients being important for CRM, the resonant contribution of these electron-ion recombination processes, dielectronic recombination (DR), often dominates at specific temperature ranges. For CRM, both RR and DR are needed: to account for this, the background contribution of RR rate coefficients is summed to the resonant contributions of DR rate coefficients in what is called the independent treatment of RR and DR [28,29,30], which assumes these are separated processes.

2.1. Calculation of RR Rate Coefficients with Atomic Codes

The quantity of interest for CRM, with respect to RR, is the RR rate coefficient, ${\alpha }_{\mathrm{RR}}$. It is the velocity-averaged RR cross-section, ${\sigma }_{\mathrm{RR}}$, from a given recombining ion’s level into a level of the recombined ion, such that [30]

$${\alpha }_{\mathrm{RR},\mathrm{lev}}\left(T\right)\equiv 〈v_{\mathrm{e}}{\sigma }_{\mathrm{RR},\mathrm{lev}}\left(v_{\mathrm{e}}\right)〉={\int }_0^{\infty }{v}_{\mathrm{e}}{f}_{\mathrm{MB}}(v_{\mathrm{e}},T){\sigma }_{\mathrm{RR},\mathrm{lev}}\left(v_{\mathrm{e}}\right)d{v}_{\mathrm{e}}.$$

(2)

The RR cross-section is typically studied as a function of the incoming electron’s velocity or energy. In atomic codes like FAC and AS, it is obtained by the principle of detailed balance. In this context, detailed balance relates RR and its time-inverse process, photoionisation (PI) [30]. This can be achieved with

$${\sigma }_{\mathrm{PI}}{\displaystyle \frac{g_f}{{p_{\mathrm{e}}}^2}}={\sigma }_{\mathrm{RR}}{\displaystyle \frac{g_i}{{p_{\gamma }}^2}},$$

(3)

in which $g\equiv 2J+1$ denotes the multiplicity of the atomic state, $p_{\mathrm{e}}$ is the electron momentum and $p_{\gamma }$ is the photon momentum. Lastly, it is necessary to define the PI cross-section, ${\sigma }_{\mathrm{PI}}$, which can be done (in atomic units) by

$${\sigma }_{\mathrm{PI}}=2\pi {\displaystyle \frac{{\left(\alpha \omega \right)}^{2k-1}}{(2k+1)g_f}}|〈f\left|\right|T^{\left(k\right)}\left|\right|i〉{|}^2,$$

(4)

where $|\langle f||T^{\left(k\right)}{\left|\right|i\rangle |}^2$ is the line strength, with $T^{\left(k\right)}$ being the multipole operator and k the rank of the multipole inducing the transition [30,31,38]; $\omega$ is the angular frequency of the ejected photon and $\alpha$ is the fine-structure constant. In this work, only electric dipole (E1) allowed transitions were considered, as these are the largest contributors in kilonovae [39].

Up until this point, FAC and AS use the same formalism. The main differences are in how these two atomic codes handle the relativistic interactions in the atomic system. A note should be made in which both FAC and AS use the Distorted Wave (DW) approximation to address this process [31,40,41]. DW is widely used to treat RR and DR as separated, independent process, in contrast with R-matrix codes that treat both process together/simultaneously [28,29], under the unified treatment of electron-ion recombination processes.

2.1.1. The Flexible Atomic Code

FAC [31] is a fully relativistic atomic code, meaning that it employs Dirac’s formalism through the Dirac-Coulomb Hamiltonian and a local central potential: a deeper description of FAC’s atomic structure formalism can be found in the literature [17,31,41,42,43,44]. FAC uses the relativistic $jj$ coupling, and describes the atomic state functions via a relativistic configuration interaction (CI) as a linear combination of configuration state functions, the mixing coefficients arising from diagonalising the Hamiltonian.

Coming back to RR and the detailed balance equation (3), to explicitly relate the RR and PI cross-section, FAC uses the relativistic Milne relation (in atomic units) [31]

$${\sigma }_{\mathrm{RR}}={\sigma }_{\mathrm{PI}}{\displaystyle \frac{g_f}{g_i}}{\displaystyle \frac{{\left(\alpha \omega \right)}^2}{\epsilon \left(\epsilon {\alpha }^2+2\right)}},$$

(5)

where $\epsilon$ is the energy of the incoming electron, and the remaining quantities have already been defined.

FAC outputs the RR and PI cross-sections. The RR rate coefficients were calculated in a post-processing code developed during this work, for a temperature range adequate to study this process in different temperature conditions ($T\in [{10}^0,{10}^9]$ K); our post-processing code was benchmarked in a previous study by Garcia et al. [45] with W^$13+$, showing good agreement with other calculations.

2.1.2. Autostructure

AS [32] is a semi relativistic atomic code. It uses a non-relativistic Hamiltonian, ${\mathcal{H}}_{\mathrm{NR}}$, with relativistic interactions as corrections, ${\mathcal{H}}_{\mathrm{RC}}$. AS’s Hamiltonian with the relativistic corrections is commonly referred to as the Breit-Pauli Hamiltonian, being thoroughly described in the literature [13,41,46,47,48,49,50]. Moreover, the Thomas-Fermi-Dirac-Amaldi (TFDA) model potential from AS [51] was used in this work; the details of the TFDA potential can also be consulted in the literature [13,30,38,46,47,50,52]. AS also uses a CI to describe the atomic state functions, although, in contrast with FAC, the angular part was computed in intermediate coupling.

Regarding radiative recombination, AS links the PI and RR cross sections also through the Milne relation, but in its non-relativistic limit [30] taking $\epsilon {\alpha }^2\ll 2$ in equation (5).

AS’s output prints the RR and PI cross-sections. The RR rate coefficients were computed using the adasrr post-processing script developed by N.R. Badnell, available online [53]. This post-processing script computes RR rate coefficients for different levels of the recombining ion over the ADAS (Atomic Data and Analysis Structure [54]) temperature range, $T\in q^2\times [{10}^1,{10}^7]$ K, with q being the ion’s ionisation (or its residual charge) [55,56,57].

2.2. The Axelrod Approximation

In 1980, T.S. Axelrod proposed in his PhD thesis [36] a semi-empirical expression to account for RR rate coefficients. Historically speaking, the astrophysics community has been exploiting this approximation due to the lack of explicit atomic calculations of RR rate coefficients, and ease of implementation [33,34,35]. Nonetheless, Mulholland et al. [58] and Ferreira da Silva et al. [59] recently showed that, for the case of electron-impact excitation processes, Axelrod’s (in that case, Van Regemorter [60] and Axelrod’s) approximation underestimates this contribution.

For comparison purposes, this work compares the RR rate coefficients explicitly calculated with FAC and AS with Axelrod’s approximation:

$${\alpha }_{\mathrm{RR}}\left(T\right)={\alpha }_0{\left({\displaystyle \frac{T}{1\times {10}^4\mathrm{K}}}\right)}^{-3/4}{q}^2,$$

(6)

with q being the ionisation stage and ${\alpha }_0$ a predetermined, tabulated value. Brethauer et al. [35] set ${\alpha }_0=3\times {10}^{-12}$ cm³/s for all lanthanides. It is worth noting that equation (6) does not depend upon the atomic number, Z, but rather only on the ionisation stage, q: by considering the same ${\alpha }_0$ for all lanthanides, we are treating similarly lanthanides with atomic structure that are known to be different [17], such as those of La–Nd, Pm–Dy, and Ho–Yb.

2.3. Fitting RR Rate Coefficients

For the sake of data replication, it is very common [55,56,61,62,63] to fit RR rate coefficients with the analytical formula of Verner and Ferland (1996) [37]

$${\alpha }_{\mathrm{RR}}\left(T\right)=A{\left[\sqrt{T/T_0}{\left(1+\sqrt{T/T_0}\right)}^{1-B}{\left(1+\sqrt{T/T_1}\right)}^{1+B}\right]}^{-1}.$$

(7)

For low-charge ions (this work), B is taken as [61]

$$B\to B+Cexp(-T_2/T),$$

(8)

which means that for the case of low-charge ions there are six free fitting parameters: A, B, C, $T_0$, $T_1$ and $T_2$. A has units of cm³/s; B and C are dimensionless; and $T_0$, $T_1$ and $T_2$ have units of temperature. This fit was performed with the Markov chain Monte Carlo (MCMC) technique, and the quality of the fit was quantified via a relative error.

3. Calculation’s Details

Doubly ionised lanthanides were studied in this work with two atomic codes: FAC and AS. RR rate coefficients are typically calculated for a transition from an initial configuration to that same configuration with an additional electron in a bound orbital that can foster it, for example

• Recombining configuration [Xe] $4f^4$
• Recombining configuration [Xe] $4f^4+\left(nlj\right)\to$ [Xe] $\left\{\begin{array}{c}4f^5\hfill \\ 4f^{4}5d^1,4f^{4}5f^1,4f^{4}5g^1\hfill \\ 4f^{4}6s^1,4f^{4}6p^1,4f^{4}6d^1,...\hfill \\ 4f^{4}7s^1,...\hfill \\ ...\hfill \end{array}\right.$    .

In this work, the additional electron $\left(nlj\right)$ in the shell with quantum numbers n, l, and j is treated explicitly until $n=10$ and $l=5$, following Trzhaskovskaya et al. (2021) [64]. We verified (with FAC) that recombination into shells with increasing l contributed residually in comparison with those of smaller l (within the same n): this happens because the RR cross-sections decrease steeply as the orbital angular momentum increases [65]. For high-n states, we assume that the captured electron’s wavefunction are highly hydrogenoid. In FAC, these high-n states were treated with Kramers’ formula [61,66,67] from $n=11$ to $n=999$; in AS, these states were treated with a hydrogenic top-up approximation [68,69,70,71,72] up to $n=999$, behaving similarly to Kramers’ approximation. These hydrogenic approximations are valid because, for sufficiently large n, the captured electron is far from the core and the corresponding wavefunctions become essentially hydrogenic.

Initially, RR rate coefficients from the ground state level were calculated with FAC and AS for some simpler doubly ionised lanthanides. As there was no data found to compare with, we performed analogous calculations in AS. The FAC radial wavefunctions and central potential, in the RR calculation, were optimised to the ground state configuration of the recombined ion. In AS, no optimisation was used, taking the Thomas-Fermi scaling parameters [13,30,46,47,50,52] equal to the unity. With this set of Thomas-Fermi scaling parameters we were able to determine the correct ground state levels, as was done by N.R. Badnell [57,73] (see Table 1). Moreover, since AS is only being used as a benchmark before proceeding to larger calculations with FAC, no further attention was directed to AS’s TFDA potential optimisation. After verifying that the ground state levels of the recombining ions were being correctly identified, a study of the atomic energy levels of the recombined side was made. The latter are then compared with those of NIST Atomic Spectra Database (ASD) [74] in Figure 2. There is a general agreement with the NIST experimental values, although a formal study on these was not performed.

Moreover and coming back to RR rate coefficients, ${\alpha }_0$ from the Axelrod approximation (equation (6)) was fitted between $T\in [500,20000]$ K, which is considered typical kilonovae temperatures, from 1 day to 1 week after the collision: this region is shaded in grey in the Figure 3, Figure 4, Figure 5and Figure 9. The fitting procedure implemented was a non-linear least squares method, taking $q=2$ in the Axelrod approximation, as doubly ionised lanthanides are the main goal of study in this work.

Next, RR rate coefficients from all levels belonging to the ground state configuration were calculated with FAC. Following the previous case, FAC’s central potential and radial wavefunctions were optimised to the ground state configuration of the recombined ion. Additionally, ${\alpha }_0$ from the Axelrod approximation (equation (6)) and the free parameters from the Verner and Ferland expression (equation (7)) were fitted to the obtained curves. The former were fitted for the same temperature range as in the previous case, using a non-linear least squares method and $q=2$; the latter was fitted for the whole temperature range ($T\in [{10}^0,{10}^9]$ K) using a MCMC fitting technique.

Finally, as for the temperatures of interest there are more levels populated besides those of the ground state configuration, RR rate coefficients from excited configurations were calculated with FAC. Here, no fit was made, but rather the contribution of the different configuration to the final curve was studied.

4. Results and Discussion

In this section we present the results from the calculations described in the previous sections, using FAC and AS, discussing them in parallel.

4.1. RR Rate Coefficient from the Ground State Level

Starting with RR rate coefficient from the ground state level, the agreement between FAC and AS is shown in Figure 3 for doubly ionised lanthanum, cerium, praseodymium, holmium, erbium, thulium and ytterbium. This comparison shows that both calculations yield RR rate coefficients in the same order of magnitude. However, there are still some deviations. We attribute these differences for high temperatures, associated with recombination into low-$nl$ states, to the differences in the levels’ energy. As was reported in previous studies [50], AS tends to poorly calculate the energies of the lowest-lying states. Furthermore, we ensured that we were taking the same levels as ground state levels (those of Table 1), as AS also tends to mix-up the ground state level [50]. Lanthanides are relativistic atomic systems and we believed that FAC’s fully relativistic approach better suits these ions than AS’s semi-relativistic approach.

Figure 3. RR rate coefficients from the recombining ion’s ground state level calculated with FAC and AS.

Figure 3. RR rate coefficients from the recombining ion’s ground state level calculated with FAC and AS.

Extending the RR rate coefficient from the ground state level calculation to all doubly ionised lanthanides (except gadolinium for computational issues, not enough RAM memory), we got the results shown in Figure 4. Figure 4 has a fitted Axelrod curve, where ${\alpha }_0$ was optimised to fit the FAC calculation in the $T\in [500,20000]$ K range, shown in Table 2: the uncertainty to ${\alpha }_0$ is given by covariance matrix of the fitting procedure. The fitting of the Axelrod approximation, equation (6), only makes vertical translations to the curve (y-intercept in log-log scale), since its slope (in log-log scale) is fixed at $-3/4$ by the approximation. Despite this, the fits are in general agreement, yielding a high coefficient of determination (Table 2).

Figure 4. RR rate coefficients from the recombining ion’s ground state level, calculated with FAC. The fitted parameters of the Axelrod approximation are shown in Table 2.

Figure 4. RR rate coefficients from the recombining ion’s ground state level, calculated with FAC. The fitted parameters of the Axelrod approximation are shown in Table 2.

Next, we moved from RR from the ground state level to RR from the ground state configuration. This means that the free electron can be captured into any level belonging to the ground state configuration, not only the ground state level.

4.2. RR Rate Coefficient from the Ground State Configuration

For RR rate coefficient from the ground state configuration, i.e. all levels from the ground state configuration, a procedure analogous to the previous one was followed. The RR rate coefficients curves for doubly ionised lanthanides, are shown in Figure 5, together with the newly fitted Axelrod approximation curves (fitted ${\alpha }_0$ presented in Table 3), and Verner and Ferland fit (fitting parameters displayed in Table 4). Similarly to before, the fitted ${\alpha }_0$ from the Axelrod approximation, equation (6), yield a high coefficient of determination, suggesting the fit is representative for data in the fitting interval ($T\in [500,20000]$ K). As for the Verner and Ferland fit, equation (7), the fits systematically yield a relative error below $4\%$, being around $2\%$ for the majority of the cases. The uncertainty of the parameters is given by the MCMC technique. We would like to mention that making this fit with the standard non-linear least squares method proved to be challenging due to the presence of local minima, and due to the method’s sensitivity regarding the initial parameters.

Figure 5. RR rate coefficients from the recombining ion’s ground state configuration levels, calculated with FAC. The fitted parameters of the Axelrod approximation are shown in Table 3. The parameters of the Verner and Ferland fit are shown in Table 4.

Figure 5. RR rate coefficients from the recombining ion’s ground state configuration levels, calculated with FAC. The fitted parameters of the Axelrod approximation are shown in Table 3. The parameters of the Verner and Ferland fit are shown in Table 4.

Despite the fact that the Verner and Ferland fit and the FAC calculations overlap (Figure 5), the relative errors below $4\%$ (shown in Table 4) can influence the light-curve luminosities. Brethauer et al. [35] showed that changes of 1 order of magnitude in Axelrod’s ${\alpha }_0$ influence the luminosity light-curves; our errors of $\pm 4\%$ correspond to $\pm 0.017$ orders of magnitude changes in the ${\alpha }_{\mathrm{RR}}\left(T\right)$ curves.

Following this, we show scaling in the different fitted parameters, both for RR rate coefficients from the ground state level and configuration of the recombining ion.

4.3. Scaling

The fitted ${\alpha }_0$ parameters from the Axelrod approximation of RR rate coefficients from the levels of the ground state configuration are larger than those of the ground state level, particularly for the lanthanides with half-open f-shell ground state configurations: see Table 2 and Table 3. This occurs because there are more levels being considered in the former than in the latter. Moreover, since most of the studied ions have complementary electronic configurations [74,75,76], there is a symmetry between the number of levels originating from the ground state configuration (see Figure 6(a)), for the lower-Z and the higher-Z doubly ionised lanthanides (Ce III and Er III, Pr III and Ho III, Nd III and Dy III, and Pm III and Tb III). By inspection of Figure 6, it can be understood how the number of levels (Figure 6(a)) influences the fitted ${\alpha }_0$ parameters (Figure 6(b)).

However, by normalising the fitted ${\alpha }_0$ of RR rate coefficient from the levels of the ground state configuration to the total number of levels in that configuration, the fitted ${\alpha }_0$ become comparable to those of the ground state level, as displayed in Figure 7. Furthermore, there is an increase of these fitted ${\alpha }_0$ with the atomic number of the lanthanides: a feature not captured by the Axelrod approximation. We believe that La III does not follow the trend in Figure 7 because its ground state configuration is not a $4f^1$, but rather a $5d^1$.

The rise of the fitted ${\alpha }_0$ values (Figure 7) can be related to the increase in the atomic number, or with the filling of the open f-shell, as can be understood with help of Table 1. If the first hypothesis holds, it would mean that, in the temperature range studied ($T\in [500,20000]$ K), ions with larger Z are more prone to capture free electrons from the continuum. This would be similar to the trends of the ionisation energy throughout the periodic table (it increases from left to right within the same period of the periodic table), but the ionisation energy in doubly ionised lanthanides does not follow the general trend in the periodic table [74]. If, instead, the filling of the open f-shell is responsible for the increase of the fitted ${\alpha }_0$ values, then it would mean that the electronic configuration to which the free electron is captured affects the probability (the cross-section and, consequently, the rate coefficients) of the process to occur (again, in the temperature range studied in this work).

It is not obvious that the fitted ${\alpha }_0$ of the Axelrod approximation for RR from the ground state level should overlap with those of the ground state configuration levels. The RR cross-section (and subsequent RR rate coefficient) of the non-zero excitation energy levels are not necessarily similar in shape, behaviour or order of magnitude between themselves and comparing with the ground state level. Additionally, it can be understood by inspection of Figure 7 that the non-zero excitation energy levels yield, on average, a larger RR rate coefficient.

Before moving on, it is worth noting that the fitted ${\alpha }_0$ parameters for the ground state level are approximately one order of magnitude smaller than those proposed by Brethauer et al. [35]. With a careful analysis of Figure 4 and Figure 5, as well as the $R^2$ coefficient shown in Table 2 and Table 3, it is easy to understand that the Axelrod approximation is a viable estimative for RR rate coefficients in the regime of lower temperatures, given an adequate choice of ${\alpha }_0$. To get an even better fit, the dependence on the temperature can also be fitted, instead of taking $T^{-3/4}$.

Figure 8 shows how the fitted six free parameters from the Verner and Ferland expression vary with atomic number. In contrast with Axelrod’s ${\alpha }_0$, these six free parameters do not exhibit a clear trend as a function of the atomic number, and the main trend is in the A parameters, as it is analogous to the ${\alpha }_{\mathrm{RR}}\left(T\right)$ curve height and, therefore, to the fitted ${\alpha }_0$ shown in Figure 6(b). As mentioned before, the uncertainties plotted are those of the MCMC technique, which yield large uncertainties for three of the cases. Several fits were performed for each ion (given the random nature of MCMC) and those with smaller relative error were chosen to present in this work.

Still in Figure 8, a visual inspection reveals possible correlations among the B, C, $T_1$, and even $T_0$ parameters. Tentatively, B and C appear to follow similar trends, while $T_1$ exhibits an inverse behaviour relative to both of them. The trend for $T_0$ is more difficult to ascertain owing to the large error bars, although it seems comparable to that of $T_1$, except for the last two data points. A closer examination of how the six free parameters influence the Verner and Ferland fit [37] shows that, although not redundant, some parameters affect the fit in analogous ways: a feature that may be reflected in the parameter correlations noted above.

4.4. Population-Weighted RR Rate Coefficients

For the kilonovae temperature range, one can make an estimate of which atomic energy levels, of the recombining ion, will be populated for a given temperature. Furthermore, it is possible to understand which configurations (as a set of levels) of the recombining ion will also be populated. To make such an estimative it is usually necessary to solve the rate equations; for simplicity, this work assumed LTE conditions, even for times in which these might not be a good approximation. In this fashion, it is possible to determine the population ($\mathcal{P}$) of the levels and configuration of the recombining ion by an adapted form of the excitation Boltzmann equation

$$\begin{array}{ccc}\hfill {\mathcal{P}}_{\mathrm{lev}}\left(T\right)={\displaystyle \frac{g_{\mathrm{lev}}{e}^{-E_{\mathrm{lev}}/k_{B}T}}{\mathcal{Z}\left(T\right)}}& \hfill & \hfill {\mathcal{P}}_{\mathrm{cfg}}\left(T\right)=\sum _{\mathrm{lev}\in \mathrm{cfg}}{\mathcal{P}}_{\mathrm{lev}}\left(T\right),\end{array}$$

(9)

where the subscripts “lev” and “cfg” denote “level” and “configuration”. $\mathcal{Z}\left(T\right)$ is the partition function, such that $\mathcal{Z}\left(T\right)={\sum }_{\mathrm{i}}^{\mathrm{all}\mathrm{bound}\mathrm{levels}}{g}_{\mathrm{i}}exp\left(-E_{\mathrm{i}}/k_{B}T\right)$. The population, $\mathcal{P}$, depends upon the temperature, T, because the higher the temperature, the more populated can high-energy atomic levels be populated.

Figure 9 shows the results of the methodology implemented to include excited configurations in the RR rate coefficients curves. This was done for doubly ionised lanthanum, cerium, praseodymium, erbium, thulium and ytterbium. By analysis of Figure 9, we can understand that the ground state configuration dominates for the lower-temperature regime, as expected. While proper non-LTE modelling is important to verify these results, we predict that for $T\lesssim 1000$ K (typical non-LTE temperatures), only the ground state configuration should be relevant.

Figure 9. Population-weighted RR rate coefficients for the most populated configurations of the recombining ion until 20000 K.

Figure 9. Population-weighted RR rate coefficients for the most populated configurations of the recombining ion until 20000 K.

As temperature increases, the other configurations become progressively more relevant. We would like to note that the RR rate coefficients of some excited configuration are larger than that of the ground state configuration simply because the RR cross sections are also larger: we did not find any particular reason for this. Moreover, we find it necessary to emphasise that the configurations shown in Figure 9 are only representative until $2\times {10}^4$ K, meaning that the populations of the configurations chosen only sum to $\approx 1$ until this temperature limit.

As the temperature reaches approximately ${10}^{3.5}$ K, Figure 9 shows that the first excited configuration tends to dominate the RR rate coefficient curve. Additionally, not only does it dominate the curve, but how much it dominates tends to increase with the atomic number. We also note that the temperature at which the bumps starts also tends to increase with the atomic number. Nevertheless, a formal analysis was not performed because the RR rate coefficients for excited configurations were only studied for six ions: these are too few to perform a systematic analysis, specially when there are no ions representative of half-open f-shells.

A second bump at higher temperatures, $\sim {10}^6$ K, can be seen in Figure 4 and Figure 5. It has been described in the literature [55,56] as recombination for low-$nl$ states, characteristic from low-charge many-electron ions. It occurs due to the highly non-hydrogenic screening of the low-$nl$ (low-energy) states’ wavefunctions, meaning these are not as effectively screened from the nucleus as higher $nl$ levels [55,56].

5. Conclusions

FAC and AS are in general agreement for simple cases of RR rate coefficients from the ground state level of doubly ionised lanthanides. We attribute the observable small deviations to the different atomic structure calculations, due to the differences in how these atomic codes include relativistic interactions in the Hamiltonian.

We fitted the RR rate coefficients with the Axelrod approximation, yielding fitted ${\alpha }_0$ parameters different from those used in previous works [25,35] by one order of magnitude. Furthermore, we verified scaling with the atomic number in the fitted ${\alpha }_0$ parameter of the Axelrod approximation, both for RR rate coefficients from the ground state level and configuration. For RR rate coefficients from the levels of the ground state configuration, the simple concept of complementary configuration played a key role in understanding the trend.

Although the parameters from the Verner and Ferland fit also tend to follow specific trends, we do not consider that such statement can be made with the data presented. These parameters are key for data replication/reproducibility.

It was not possible to perform the RR calculation for Gd III due to lack of computational memory (RAM). Its $4f^{7}5d^1$ ground state configuration generates lots more levels than a $4f^8$ configuration. Future work might lay in calculating only $4f^1$ (in lanthanum) and $4f^8$ (in gadolinium), and verify if these follow the trends shown in this work. Such analysis was not performed, as the present study was limited to configurations populated in the kilonova temperature range.

Lastly, we propose a procedure to account for excited configurations. However, the main message is that the RR cross-sections for some excited configurations are larger than those of the ground state configuration, also yielding larger RR rate coefficients. This means that as higher configuration become populated, these might dominate in specific temperature ranges because the RR cross-section is higher than those of the remaining configurations. It was also understood that to model the late epochs of kilonovae, in temperatures below 1000 K, the ground state configuration yields the most representative contribution to the RR rate coefficient curve.

Following the independent treatment of RR and DR, $\alpha \left(T\right)={\alpha }_{\mathrm{RR}}\left(T\right)+{\alpha }_{\mathrm{DR}}\left(T\right)$, the inclusion of the resonant dielectronic recombination in the total recombination rate coefficient curves is left our for a future work currently in preparation, as a Part II of this work.