Opacity of the highly ionized lanthanides and the effect on the early kilonova

We investigate the effect of the presence of lanthanides (Z = 57- 71) on the kilonova at t~hours after the neutron star merger for the first time. For this purpose, we calculate the atomic structures and the opacities for selected lanthanides: Nd (Z = 60), Sm (Z = 62), and Eu (Z = 63). We consider the ionization degree up to tenth (XI), applicable for the ejecta at t ~ a few hours after the merger, when the temperature is T ~ 10^5 K. We find that the opacities for the highly ionized lanthanides are exceptionally high, reaching k_exp~1000 cm^2/g for Eu, due to the highly dense energy levels. Using the new opacity, we perform radiative transfer simulations to show that the early light curves become fainter by a (maximum) factor of four, in comparison to lanthanide-free ejecta at t~0.1 day. However, the period at which the light curves are affected is relatively brief due to the rapid time evolution of the opacity in the outermost layer of the ejecta. We predict that for a source at a distance of ~100 Mpc, UV brightness for lanthanide-rich ejecta shows a drop to ~21-22 mag at t~0.1 day and the UV peaks around t~0.2 day with a magnitude of ~19 mag. Future detection of such a kilonova by the existing UV satellite like Swift or the upcoming UV satellite ULTRASAT will provide useful constraints on the abundance in the outer ejecta and the corresponding nucleosynthesis conditions in the neutron star mergers.


INTRODUCTION
It has long been hypothesized that the heavy elements are synthesized in the neutron star mergers (e.g., Lattimer & Schramm 1974;Eichler et al. 1989;Freiburghaus et al. 1999;Korobkin et al. 2012;Wanajo et al. 2014). The radioactive decays of the freshly synthesized heavy elements give rise to a transient in the ultraviolet, optical, and near-infrared wavelengths, called kilonova (e.g., Li & Paczyński 1998;Kulkarni 2005;Metzger et al. 2010). In fact, such a kilonova AT2017gfo (e.g., Coulter et al. 2017;Yang et al. 2017;Valenti et al. 2017) from the neutron star merger has already been observed by the follow up observations of the gravitational wave event GW170817 (Abbott et al. 2017), confirming the neutron star mergers to be a site of the r-process nucleosynthesis.
The light curve of AT2017gfo was bright in UV and optical band at the epoch of detection (e.g., Coulter et al. 2017;Yang et al. 2017;Valenti et al. 2017). The light curve evolved to be fainter at optical and brighter at NIR in the timescale of t ∼ a week (e.g., Cowperthwaite et al. 2017;Smartt et al. 2017;Drout et al. 2017;Utsumi et al. 2017). The late time (t > 1 day) light curve is well explained by kilonova (e.g., (Kasen et al. 2017;Tanaka et al. 2017;Shibata et al. 2017;Perego et al. 2017;Rosswog et al. 2018;Kawaguchi et al. 2018)). However, the origin of the early emission (t < 1 day) has not reached a consensus (Arcavi 2018). Several models exist for the early kilonova. For example, the early emission might be powered by the radioactive decays of the heavy elements, similar to the later phase (Waxman et al. 2018;Villar et al. 2017;Banerjee et al. 2020). Alternatively, the early emission might result from the interaction of the relativistic jet with the surrounding ejecta (Kasliwal et al. 2017;Piro & Kollmeier 2018;Nativi et al. 2021;Klion et al. 2021) or by β-decay of the free-neutrons (Metzger et al. 2015). It is important to understand the early kilonova because the early emission can reveal the abundance in the outer ejecta because photons only from the outer layer can escape at an early time. Hence, understanding the kilonova starting from an early time is crucial to understanding the abundance pattern from the outer to the inner ejecta.
The major uncertainty in modelling the early kilonova comes from the lack of the detailed opacity of the r-process elements. The shape of the kilonova light curve is mainly determined by the opacity in the ejecta (e.g., Metzger et al. 2010;Kasen et al. 2013;Tanaka & Hotokezaka 2013). Hence, modelling the kilonova requires detailed opacity. Past studies have shown that the bound-bound transitions contribute the most to the opacity in the kilonova (Kasen et al. 2013;Tanaka & Hotokezaka 2013;Fontes et al. 2015Fontes et al. , 2020Wollaeger et al. 2017;Tanaka et al. 2018Tanaka et al. , 2020. Calculations of the bound-bound opacity require the atomic data, which are largely unavailable for the early kilonova. This is due to the fact that, at an early time, the elements are highly ionized, maximum up to tenth ionization (or XI in spectroscopic notation; hereafter, the spectroscopic notation is used to describe the ionization) at T ∼ 10 5 K, the typical temperature at t ∼ 0.1 day. Most of the earlier studies have performed the atomic calculations only up to the ionization IV, which are suitable for the opacity at t ≥ 1 day. Hence, to derive opacity at the early time, the atomic structure calculations for the highly ionized heavy elements are necessary. The expansion opacity as a function of wavelength at T = 70000 K (left panel) and the Planck mean opacity as a function of temperature (right panel) at ρ = 10 −10 g cm −3 , and t = 0.1 day for lanthanides Nd, Sm, and Eu. Opacity of the light r-process element Cd is also shown for comparison (Banerjee et al. 2020). Banerjee et al. (2020) performed the atomic calculations for the light r-process elements at the ejecta condition suitable at the early time. However, their study does not include lanthanides (elements with Z = 57 − 71). Lanthanides are the elements with open 4fshell. Previous studies on the opacity at t ≥ 1 day have shown that the presence of lanthanides can significantly affect the opacity and the light curve.
In this paper, we perform the first atomic opacity calculation for the selected lanthanides: Nd (Z = 60), Sm (Z = 62), and Eu (Z = 63) up to the ionization XI and study the impact on the early kilonova emission. We show our new atomic and opacity calculations in Section 2. Using the new opacity, we perform radiative transfer simulations for lanthanide-rich ejecta in neutron star mergers in Section 3. The validity of standard method of opacity calculation in the expanding media (expansion opacity formalism, Sobolev 1960) for the highly ionized lanthanides is discussed in Section 4.1. We also investigate the future prospects to detect lanthanide-rich kilonova in Section 4.2. Finally, we summarize our conclusions in Section 5. AB magnitude system is adopted throughout the article.

OPACITY IN NEUTRON STAR MERGER
In the neutron star merger ejecta, different processes such as electron scattering, free-free, bound-free, and the bound-bound transitions contribute to the opacity. Here we calculate the different opacity components for the selected lanthanides, Nd (Z = 60), Sm (Z = 62), and Eu (Z = 63) following Banerjee et al. (2020). The ejecta conditions are assumed to be suitable for the early time, i.e., the density of the ejecta is taken to be ρ = 10 −10 g cm −3 , which is the typical density at t ∼ 0.1 day for an ejecta mass of M ej ∼ 0.01M and the elements are considered to be ionized up to ∼ XI corresponding to the typical temperature of T ∼ 10 5 K.
First, we estimate the electron scattering, free-free, and the bound-free opacities for lanthanides. We find that the electron scattering and the free-free opacity for lanthanides (ionized up to ∼ XI) are ∼ 3 ×10 −2 cm 2 g −1 and ∼ 2 ×10 −3 cm 2 g −1 , respectively, which are not very different from those for light r-process elements. Note that we assume the electron temperature is the same as the ejecta temperature. Also, we confirm that the bound-free opacity is not important in our chosen wavelength range (λ = 100 − 35000Å). This is because the fraction of the photons with energy beyond the ionization potential is not significant. A similar conclusion was made for the early bound-free opacity for the light r-process elements (Banerjee et al. 2020).
Next, we explore the bound-bound opacity for lanthanides, for which the atomic structure calculations are necessary. Since the complete atomic data calibrated with experiments are unavailable for the highly ionized lanthanides, we perform the theoretical atomic structure calculations as described in the following section. We use HULLAC (Hebrew University Lawrence Livermore Atomic Code, Bar-Shalom et al. 2001) for the atomic calculations. HULLAC uses fully relativistic orbitals to calculate the energy levels and radiative transition probabilities. A set of orbital functions is obtained by solving the single electron Dirac equation with a parametric, central field potential, which includes both nuclear field and the spherically averaged electron-electron interaction. The central field potential for a given electron charge distribution can be obtained by solving the Poisson equation with the boundary condition that the potential converges to (Z − q)/r at r = ∞ for an element with atomic number Z and the number of the electrons q. The nuclear charge is assumed to be a point one (Zδ(r)). The charge density distribution of the electrons is expressed as (Bar-Shalom et al. 2001): where N is the normalization factor and α is a free parameter related to the average radii (< r >) of the Slater type orbital as (2l + 3)/ < r >. The free parameter α is obtained by minimization of the first order configuration averaged energies for selected configurations. The all-electron zero-order solution or the configuration state function (CSF) is the anti-symmetrized products of the orbitals in a chosen coupling scheme (jj coupling in this case). After constructing the zero-order wave function, the magnetic and the retardation effect of the interaction from the other electrons (Breit term) and the quantum electrodynamic energy correction are taken into account. Finally, the atomic wave function is constructed using the linear combination of the CSFs.
We perform the atomic calculations including the excited configurations for the opacity. Such calculations require the ground configurations. This is because HUL-LAC solves the Dirac equation with a certain central potential derived by the electron distribution in the ground state. However, the ground configurations are not well established for the highly ionized lanthanides (ionization ≥ V). Therefore, for the three lanthanides (Nd, Sm, and Eu) ionized to the degree V -XI, we use the ground configurations estimated with HULLAC. The detailed strategy of the calculations to estimate the ground configurations is provided in Appendix A. For the all the other ions, we use the ground configurations as provided in NIST atomic spectra database (Kramida et al. 2020). Using the ground configurations together with the excited configurations, the atomic calculation for the opacity is performed. We provide all the configurations used in Table A4. The ground configurations and all the configurations used for minimization are shown in bold.
It is noted that assuming the central potential has a parametric form (Equation (1)) in HULLAC can be a potential source of uncertainty in atomic calculations. Different atomic codes adopt different approaches: for example, Flexible Atomic Code (FAC, Gu 2008) determines the potential in a self-consistent manner. Nevertheless, comparison of atomic data calculated by HUL-LAC and more ab-initio GRASP2018 code (Froese Fischer et al. 2019) shows good match for various ions with ionization < V (e.g., Tanaka et al. 2018;Gaigalas et al. 2019;Rynkun et al. 2022). Similar comparison for highly ionized lanthanides remains within the scope for the future work.

Energy level distribution
The energy level distribution obtained from our atomic calculations is shown as a function of ionization in Figure 1. The color scale represents the level density in the 0.2 eV energy bin. Since our main purpose is to calculate the bound-bound opacity, we show the energy levels only below the ionization threshold. Also, the energy level distribution of the light r-process element Cd (Z = 48, Banerjee et al. 2020) is shown for comparison.
The energy level distribution of lanthanides is denser in comparison with the light r-process element Cd (Figure 1). This behavior is common over all the ionizations (for the low-ionized cases, see Kasen et al. 2013;Tanaka & Hotokezaka 2013;Fontes et al. 2015Fontes et al. , 2020Wollaeger et al. 2017;Tanaka et al. 2018Tanaka et al. , 2020. For highly ionized lanthanides, this trend is more prominent in the cases of the Sm and Eu than Nd. This is due to the large number of available states for the electrons (high complexity measures) in the middle of lanthanide series.
For a given element, the density of the energy levels increases with the ionization. The energy distribution is densest around the ionization of VII -IX for lanthanides. For instance, for Eu, the number of energy levels increases significantly from the ionization VI -VIII. This is caused by the presence of the two open shells (4f -and 5p-shell, see Table A4) in highly ionized lanthanides. We can understand this by comparing the ground configurations of Eu VI and VIII. Eu VI has only one open shell in the ground configurations (4f 4 5p 6 , Table A4), whereas Eu VIII has two open shells (4f 4 5p 4 , Table A4). The same is true for the excited configurations in both ions. Hence, Eu VIII has a higher energy level density than Eu VI.
For a given ionization, the energy level distribution varies depending on the elements. This is due to the difference in the number of 4f -electrons (Table A4). For example, Eu IX has four 4f -electrons in the ground configuration (4f 4 5p 3 , Table A4), whereas Nd IX has only one (4f 5p 3 , Table A4). As a result, Eu IX shows a higher energy level density than Nd IX ( Figure 1).

bound-bound opacity
Using the new atomic data, we calculate the boundbound opacity for the three highly ionized lanthanides (Nd, Sm, and Eu). In an expanding media such as that in the cases of supernova and kilonova, the opacity is calculated by adopting the expansion opacity formalism ( Karp et al. 1977). To calculate the expansion opacity, the contribution from the multiple lines is averaged within a chosen wavelength bin. The strength of the individual transition is calculated by assuming Sobolev approximation (Sobolev 1960).
In the expanding media, the photons are continuously redshifted, causing the photons to progressively coming into resonance with different lines. Note that the resonance occurs over a certain wavelength range rather than at a particular wavelength because of the intrinsic profiles of the atomic lines. In neutron star merger ejecta, the line profile is predominantly determined by the thermal motion. If the thermal widths of the lines are negligible in comparison with the line spacing, we can evaluate the strength of each transition by the Sobolev optical depth, which is not dependent on the line profile function. More on the validity of the Sobolev approximation is discussed in Section 4.1.
When the Sobolev approximation is valid, the expansion opacity is calculated as follows. If there are N strong lines inside an arbitrarily chosen wavelength bin of ∆λ, the velocity gradient required to redshift the photons from one line to another is given as: Such a velocity gradient corresponds to a mean free path of ∆vt at a time t. The corresponding absorption coefficient within the wavelength bin of λ to λ+∆λ is written as (Kasen et al. 2013): In this expression, only strong lines are considered. To include the contribution from the weak lines, a modified version derived by Eastman & Pinto (1993) is used: where λ l is the transition wavelength in a chosen wavelength interval of ∆λ. The Sobolev optical depth at the transition wavelength (τ l ) is calculated as where n l is the number density of the lower level of the transition, and f l is the oscillator strength of the transition. Then, we can calculate the expansion opacity as the absorption coefficient per unit mass density: Using this formalism, we calculate the opacity for a single element ejecta with the density ρ = 10 −10 g cm −3 at t ∼ 0.1 day. We assume local thermodynamic equilibrium (LTE) to calculate the ionization fraction of the elements by solving the Saha ionization equation and to determine the population of the excited levels via Boltzmann statistics.
The low density of the neutron star merger ejecta (ρ ∼ 10 −10 g cm −3 even at t ∼ 0.1 day) is not enough to establish LTE via collisional processes alone. Nevertheless, LTE can be established via radiative transitions in the optically thick regions inside the photosphere, especially in the early time when most of the ejecta are optically thick.
LTE might not be valid if the non-thermal processes from the radioactive decay significantly affect ionization and excitation. However, Kasen et al. (2013) find that the ratio of the non-thermal to thermal excitation rate at t ∼ 1 day, when the radioactive power released per particle is ∼ 1 eV s −1 and the typical transition energy is ∼ 1 eV at the temperature T ∼ 5000 K is negligible (the ratio is ∼ 10 −8 ). Extending the calculation to early time at t ∼ 0.1 day, when the radioactive power released per particle is ∼ 100 eV s −1 (Metzger et al. 2010) and the typical transition energy can be as high as ∼ 10 eV at the temperature T ∼ 10 5 K, we find the ratio is not significant (a rough estimate shows the ratio is about ≤ 10 −8 ). A similar argument can be made for the non-thermal ionization at the early time. Hence, it is expected that the non-thermal processes do not make the system largely deviated from LTE at the timescale of interest. At a later time, when the ejecta become less dense and more transparent, larger deviation from LTE is expected (for more discussion on non-LTE opacity, see Pognan et al. 2022).
The expansion opacity for lanthanides are exceptionally high (left panel of Figure 2). For instance, the expansion opacity at its peak reaches κ exp ∼ 1000 cm 2 g −1 for Eu at T ∼ 70000 K. On the other hand, the opacity of the light r-process element Cd can reach only up to κ exp ∼ 1 cm 2 g −1 under the same condition. This is due to the significantly higher number of the energy levels in the highly ionized lanthanides (Figure 1).
The expansion opacities show a strong wavelength dependence, with a higher value at shorter wavelengths (left panel of Figure 2). This is caused by the larger number of transitions at the shorter wavelengths. Moreover, the opacity for lanthanides show distinct peaks at short wavelengths (e.g., see λ ∼ 500Å and λ ∼ 1200Å in Figure 2), which is due to the fact that there are numerous strong transitions at these wavelengths.
The temperature dependence of the expansion opacity is estimated by convolving it with the blackbody function to calculate the Planck mean opacity (see right panel of Figure 2). The Planck mean opacities for different elements show distinct peaks at temperatures T ∼ 5000 K and T ∼ 70000 K. At the high temperature, Eu has the maximum opacity among the other two lanthanides. On the other hand, the opacity for Nd is highest at low temperature.
The opacity at high temperature reflects the density of energy levels. At the temperature T ∼ 70000 K, where the opacity peaks appear, lanthanides are ionized to ∼ VII -IX. At this ionization range, the level density of lanthanides is the highest (Section 2.1). We argue that, at high temperature, relatively high energy levels contribute to the opacity. This is in contrast with the lower ionized (i.e., low temperature) case, where mostly the lower lying energy levels are important (see Tanaka et al. 2020). This is because, at low temperature, only the low lying levels are populated (by Boltzmann distribution). However, at high temperature, even the relatively higher lying energy levels can be populated, making the transitions between the high energy levels possible. Hence, the density of the levels in a wider energy range is important at high temperature.
It is worth mentioning that the opacity is affected by the completeness of the atomic data as our results show that even higher lying levels are important for the opacity at high temperatures. Hence, we investigate whether our atomic data include essential transitions (i.e., whether our atomic data are sufficiently complete for the opacity). We find that the atomic data are mostly complete for opacity. More details can be found in Appendix B.

Model
In neutron star merger, the heaviest elements are produced mainly in the tidal ejecta component, i.e., in the ejecta that is mostly distributed towards the equatorial plane (e.g., Bauswein et al. 2013;Just et al. 2015;Sekiguchi et al. 2015;Kullmann et al. 2022;Just et al. 2022). Hence, the kilonova observed in the equatorial direction is likely to show the effect of the presence of lanthanides. We calculate the light curve for such a kilonova from a neutron star merger using a timeand wavelength-dependent Monte Carlo radiative transfer code (Tanaka & Hotokezaka 2013;Tanaka et al. 2017;Kawaguchi et al. 2018). The code calculates the multicolor light curves and spectra for a given a density structure and electron fraction (Y e ) distribution assuming the homologously expanding motion of the ejecta. The radioactive heating rate of r-process nuclei is calculated according to Y e , by using the results from Wanajo et al. (2014). The code adopts a time-dependent thermalization factor from Barnes et al. (2016). Our simulation considers the wavelength range λ ∼ 100 − 35000Å.
We adapt a spherical ejecta model (Metzger et al. 2010) with a power-law density structure ρ ∝ r −3 with a velocity range of v = 0.05c to 0.2c and a total ejecta mass of M ej = 0.05M (same as the fiducial model of Banerjee et al. 2020). The abundance in the ejecta is assumed to consist of the single lanthanide element (Nd, Sm, or Eu, considered as different models) with a mass fraction X La = 0.1. Such a lanthanide fraction is obtained in an ejecta with Y e = 0.20, the typical value for Y e for the equatorial ejecta. The remainder of the ejecta are considered to have the light r-process abundance. The abundance for the light r-process elements is determined by using the results from Wanajo et al. (2014) for Y e = 0.30 − 0.40. A flat mass distribution is considered for each value in the Y e range. We renormalise the abundance to match the total mass fraction of the light r-process elements to be 0.9. Note that we assume the heating rate is only from the light r-process elements because including single lanthanide does not change the heating rate significantly.
Here we mention that performing the radiative transfer simulation using the complete linelist is not feasible since the number of transitions is extremely high. For instance, the linelist for Eu can consist up to ∼ 0.3 billion lines (Table A4). In contrast, the total number of lines for the light r-process abundance is ∼ 10 million. Hence, we make a reduced linelist for lanthanide elements Sm and Eu for the ionization > V. For this purpose, we randomly select the transitions from the original linelist by keeping the statistical properties the same. The detailed scheme and the validity are discussed in Appendix C.

Results
Figure 3 (right panel) shows the bolometric luminosities for different models. We find that at t ∼ 0.1 day, the bolometric luminosities for lanthanide-rich ejecta (L bol ∼ 0.5 − 1 × 10 42 erg s −1 , different for different models) is fainter than lanthanide-free ejecta (L bol ∼ 2 × 10 42 erg s −1 ) by a factor of four to two, depending on the models. The light curves for lanthanide-rich ejecta rise afterwards and show no difference with that of lanthanide-free case at t ∼ 0.4 day. Finally, the luminosities for lanthanide-rich models drop and show deviation again at around t ∼ 1 day.
The shape of the early bolometric light curve is determined by the opacity in the outermost layer in the ejecta (v ≥ 0.19c). This is because, in the early time, the diffusion sphere lies at the outermost layer of the ejecta. At t = 0.1 day, the temperature of the outermost layer provides the suitable condition (T < 70000 K) to reach the ionization range of VII -IX, where the opacity peaks appear for lanthanides (left panel of Figure 3). Such a rise in the opacity in the outermost layer causes the luminosity to drop at t ∼ 0.1 day.
As the ejecta expands further, the temperature in the outer layer decreases, crossing the temperature (and ionization) range where the opacities peak. This results in the luminosity rising at t ∼ 0.4 day. Finally, at around t ∼ 1 day, the outermost layer of the ejecta cools down enough (T < 10000 K) so that the opacity peaks of the low-ionized lanthanides appear (see left panel of Figure  3). Hence, the luminosities drop again at t ∼ 1 day. Hence, we see that the extremely high opacities for the highly ionized lanthanides affect the light curves only for a brief period of time (t ≤ 0.4 day), reflecting the rapid temperature evolution.
Different lanthanides show the different extent of drops in the luminosity at different times. For instance, at t = 0.1 day, the luminosity of Eu-rich ejecta is the most affected, whereas at t = 1 day, the drop is the most significant for the ejecta containing Nd. The extents of the drop are determined by the peak opacities for lanthanides at different temperatures (right panel of Figure 2). For instance, at high temperature, i.e., the condition at t = 0.1 day, the opacity is maximum for Eu. In contrast, at a low temperature, i.e., the condition at t = 1 day, Nd has the maximum opacity. This explains the light curve shapes in the presence of the different lanthanides. Figure 4 shows the typical spectra in the presence of lanthanides at t = 0.1 and t = 1 days. The spectra are almost featureless. Since the accuracy of our atomic data is not enough and we use the reduced linelist, which affects the detailed spectral feature, we do not attempt to discuss the individual elemental signature. Instead, we note that the spectra rapidly evolves from UV to the optical wavelength range from t = 0.1 to t = 1 days, consistent with the expectation. The same trend is observed in all models.

DISCUSSIONS
Our work shows that the opacities for the highly ionized lanthanides are exceptionally high owing to the extremely dense energy levels. Moreover, we show the luminosity is suppressed in the early time in the presence of lanthanides in the ejecta. In this section, we discuss the validity of the expansion opacity formalism for the highly ionized lanthanides. We also discuss the future detection prospects for the early kilonova from lanthanide-rich ejecta.

Validity of the Sobolev approximation
The expansion opacity at a chosen wavelength interval of ∆λ is derived by taking the cumulative contribution from all the lines inside ∆λ with the assumption that there is sufficient space between the strong lines. If the intrinsic line profiles overlap, i.e., if the line spacing (∆λ line = ∆λ/N ) is comparable to the thermal width of the line (∆λ th ), such treatment can not represent the true opacity. Since our calculations show the high number of transitions caused by the high density of the lines (Section 2.1), we check whether the line profiles do not overlap and if the Sobolev approximation is valid.
By following Kasen et al. (2013), we define a critical opacity when the thermal width of the line is equal to the line spacing, i.e., ∆λ th = ∆λ/N = ∆λ line . Under such condition, the velocity required to redshift the photon between two consecutive lines (Equation (2)) can be given by the thermal velocity v th : The value of the v th is ∼ 4km s −1 at the typical temperature T ∼ 10 5 K at the time t ∼ 0.1 day using an average atomic mass number of A = 150. Using Equation (7), Equation (4), and Equation (6), we can derive the critical opacity as follows: The critical opacity is independent of the wavelength, while it depends on the temperature and density of the ejecta at a particular time. If the expansion opacity exceeds the critical opacity, the intrinsic line spacing is smaller than the line width, i.e., the lines overlap with each other. Then the expansion opacity using the Sobolev approximation for the radiative transfer can not represent the true opacity in the expanding media.
At a time t ∼ 0.1 day after the neutron star merger, for a density ρ ∼ 10 −10 g cm −3 , and a temperature of T ∼ 70000 K, κ crit = 3 cm 2 g −1 as shown in Equation (8). The expansion opacity, under the same condition, can reach up to κ exp ∼ 1000 cm 2 g −1 at far-UV (λ ≤ 2000Å, left panel of Figure 2), exceeding the value of the critical opacity. Hence, using the expansion opacity at t ∼ 0.1 day for lanthanides can not represent the true opacity at far-UV. Consequently, the light curves are possibly affected in the far-UV wavelengths. Nevertheless, our calculation is most likely to remain unaffected at λ ≥ 2000Å, which is the detection range of the existing UV instruments like Swift (Roming et al. 2005). The alternative treatment of opacity calculation and its implication to the early kilonova light curve will be discussed in future work.

Future prospects
In this section, we discuss the prospects of observing an early kilonova from lanthanide-rich ejecta. Figure 5 shows the magnitudes in the three different Swift UVOT filters (Roming et al. 2005) for a source at 100 Mpc for different models. Our results show that the UV brightness for lanthanide-rich ejecta drops at t ∼ 0.1 day, reaching ∼ 21 − 22 mag depending on models. The brightness increases afterwards, reaching ∼ 19 mag at t ∼ 0.2 day. Finally, the brightness decreases to > 22 mag after t ∼ 1 day.
The extents of drop and the slope of the light curves are different for the different models. For instance, at t ∼ 0.1 day, the magnitudes for Eu-rich ejecta are the faintest, whereas Nd-rich ejecta show faintest magnitudes at t ∼ 1 day. Moreover, the presence of Eu makes the light curve rise faster at t ∼ 0.1 day, whereas the presence of the Nd makes the light curve fall faster at t ∼ 1 day. This is because of the differences in the opacity in the outermost layer in the presence of the different lanthanides, as discussed in Section 3.2. The early UV signals for lanthanide-rich ejecta are bright enough to be detected by Swift (with a limiting magnitude of 22 mag for an exposure time of 1000 s, Roming et al. 2005), provided the kilonova is discovered early enough so that the prompt observation can be started. Such a kilonova is also a good target for the upcoming wide-field UV satellite ULTRASAT (limiting magnitude of 22.4 mag for 900 s of integration time, Sagiv et al. 2014). Future detection of such a kilonova will provide the clues to the abundance pattern in the outer layer of the ejecta, which can give useful constraints on the nucleosynthesis condition in the neutron star merger. Moreover, detecting a rapid rise at t ∼ 0.1 day in UV is likely to be an indicator of the presence of lanthanide element with a similar property to Eu. On the other hand, detection of a rapid decline at t ∼ 1 day in UV will likely to point that Nd-like lanthanide element is present in the ejecta.
We note that our calculation is limited to provide the trustable result only for the observation at the λ ≥ 2000Å at t ∼ 0.1 day. This is because we calculate the opacity using the expansion opacity formalism which not valid at λ ≤ 2000Å at t ∼ 0.1 day (Section 4.1).

CONCLUSIONS
To investigate the early kilonova emission from lanthanide-rich neutron star merger ejecta, we perform the atomic opacity calculation for the three lanthanides Nd (Z = 60), Sm (Z = 62), and Eu (Z = 63). For the atomic calculation, we consider the ionization up to XI, which is the maximum ionization at a typical condition of T ∼ 10 5 K at t = 0.1 day. Our opacity calculations with the new atomic data show that lanthanide opacity can be exceptionally high, reaching κ exp ∼ 1000 cm 2 g −1 for Eu (left panel of Figure 2), due to the dense energy levels in the highly ionized lanthanides (Figure 1).
Using the new opacity, we perform the radiative transfer simulations to calculate the early kilonova from lanthanide-rich neutron star merger ejecta. Our models assume that the abundance of the ejecta is the mixture of the single lanthanide (Nd, Sm, or Eu with a fraction of X La = 0.1) and the light r-process elements. Such a lanthanide-rich kilonova may replicate the ejecta condition for a kilonova observed at equatorial direction.
We find that in the presence of lanthanides, the bolometric light curves show a brief period of luminosity drop at t ∼ 0.1 day by a (maximum) factor of four in comparison to lanthanide-free case. The luminosities rise to the same value as lanthanide-free case at t ∼ 0.4 day and finally drops again at t ∼ 1 days (right panel of Figure 3). The shape of the light curve is determined by the opacity in the outermost layer in the ejecta. The opacity there changes as the the temperature (ionization) changes with the expansion of the ejecta (left panel of Figure 3). The extents of the luminosity drop are different depending on lanthanide element present since the maximum opacity are different for different lanthanides (Section 2.2).
The UV light curves show the same trends as the bolometric light curves. For a source at 100 Mpc the UV brightness drops to ∼ 21 − 22 mag at t ∼ 0.1 day (Figure 5). The brightness increases afterwards to reach ∼ 19 mag at t ∼ 0.2 day, beyond which, the brightness decreases to > 22 mag around t ∼ 1 day ( Figure 5). The extents of drop and the slopes of the light curves are different for the different models. We show that it is possible to detect the early kilonova even for lanthaniderich ejecta by Swift (Roming et al. 2005), if the kilonova is discovered early enough. Also, such a kilonova can be detected using the upcoming wide-field UV satellite ULTRASAT (Sagiv et al. 2014). Detection of such a kilonova in the early time will provide the abundance pattern in the outer ejecta and can put constraints on the nucleosynthesis condition in the neutron star mergers.
We note that our spherical ejecta model with a homogeneous abundance pattern is designed to represent the kilonova observed in the equatorial direction is simple. In reality, the ejecta structure is likely to be more complicated. For example, many simulations predict the presence of a faster moving outer layer (e.g., Kyutoku et al. 2014). In such cases, although the unique features in the light curves will still be present, the significance and the timescale might be different from those observed in our models. Exploring such possibilities in the future is of interest. Moreover, it is emphasized that our models can provide the trustable light curves only at λ ≥ 2000Å. This is because we use the expansion opacity which is not justified in the far-UV wavelengths (λ ≤ 2000Å) if the highly ionized lanthanides are present. The alternative opacity treatment and the effect on the early kilonova will be explored in future work.

ACKNOWLEDGMENTS
Numerical simulations presented in this paper were carried out with Cray XC-B at the Center for Computational Astrophysics, National Astronomical Observatory of Japan; and at the computer facility in the Yukawa Institute for Theoretical Physics (YITP), Kyoto University, Japan. The authors would like to thank the anonymous reviewer for his/her comments, which helped to improve the manuscript. SB wants to thank S. Saito for useful discussion. This research was supported by the Grant-in-Aid for Scientific Research from JSPS (19H00694,20H00158,21H04997) and MEXT (17H06363).

A. GROUND CONFIGURATIONS OF LANTHANIDES
In this section, we describe our strategy to estimate the ground configurations for the highly ionized lanthanides. In the highly ionized lanthanides (≥ V), the outermost shells are either 4f -or 5p-shell. Hence the electrons are removed either from 4f -or 5p-shell for further ionization. NIST atomic spectra database (NIST ASD, Kramida et al. 2020) provides the ground configurations for highly ionized lanthanides assuming 5p-electron is removed for further ionization in most of the cases for ≥ V . However, such results are based on the relatively simplified theoretical calculations with different approximations (Carlson et al. 1970;Rodrigues et al. 2004;Sugar & Kaufman 1975;Martin et al. 1978). Thus, we calculate the ground configurations of highly ionized lanthanides using HULLAC.
We prepare different test cases for the atomic calculations in the following way. HULLAC calculates the atomic orbitals by solving the Dirac equation with a central potential, which is determined and optimized based on the given electron distribution. We perform the atomic calculations for various central potentials: (1) calculated by changing the electron distribution in 4f and 5p orbitals and (2) optimizing for energy levels belonging to different set of configurations. For each case, we identify the CSF with the largest mixing coefficient for the lowest energy level. Then, the configuration generating the CSF of the ground level is taken to be the ground configuration. If a certain configuration appears to be the ground one in all the calculations for a particular ion, we regard it as the ground configuration.
We explain our strategy by taking Eu as an example. The ground configuration of Eu IV as suggested by NIST ASD is 4f 6 5p 6 (experimentally verified). We test whether the ground configuration of Eu V are 4f 5 5p 6 (as provided in the NIST ASD) or 4f 6 5p 5 , corresponding to the 4f and 5p electron removal from Eu IV ion. As the total number of electrons in 4f and 5p orbitals is 11, the effective potential on a single electron is constructed based on the distribution of 10 electrons. If the ground configuration is given as 4f 5 5p 6 , the single electron potentials are constructed based on the electron distributed either as 4f 5 5p 5 or 4f 4 5p 6 . Then we optimize the potential for 4f 5 5p 6 . We denote these two cases as Case A and Case B in Table A2. Similarly, if the ground configuration is 4f 6 5p 5 , the potential can be constructed for the electrons distributed either as 4f 6 5p 4 or 4f 5 5p 5 , and the potential can be optimized for 4f 6 5p 5 . These are denoted as Case C and Case D. The energy levels are calculated for both candidate configurations in all the cases. The atomic calculations with these four different cases show that the ground state always belongs to the configuration 4f 5 5p 6 . Therefore, we regard 4f 5 5p 6 as the ground configuration for Eu V.
Similar calculations for the other ions of Eu show the convergence of the ground configurations for most of the ions (shown in bold in Table A1). For some ions, however, the results of the four cases do not converge. An example is Eu VII as shown in Table A3. In such a case, we perform another set of calculation by employing both candidate configurations for the energy minimization (Cases A , B , C in Table A3). If the ground configurations converge in these cases (A , B , C ), we choose that configuration as the ground configuration. In the case of Eu VII, 4f 4 5p 5 is the ground configuration in all the Cases A , B , C , and hence, it is regarded as the ground configuration. The ions which require these additional calculations are indicated with asterisk in Table A1. The same strategy is applied to evaluate ground configurations of other highly ionized lanthanides (Nd (Z = 60) and Sm (Z = 62)). All the ions reach convergence either after Cases A-D or Cases A , B , C .
The ground configurations obtained from HULLAC are different from those provided in the NIST ASD, which shows that, 5p-electron removal starts from ion V for further ionization. In contrast, our results show that 5p-electron removal starts from ionization VI or VII, depending on the elements. For ionization VI and higher, NIST ASD adopts the results from the theoretical calculations by Carlson et al. (1970), which provides the ground configuration by removing the consecutive least bound electron. The least bound orbitals are determined from the solution of the relativistic Hartree-Fock wavefunction for the neutral atoms. On the contrary, we calculate atomic energy levels for individual ions with an effective central field potential by taking electron-electron interaction into account. Hence, we choose to use the ground configurations obtained from the HULLAC instead of the NIST ASD for the final calculation for the opacity including excited confiurations (Table A4). Table A1. Summary of ground configuration calculation for Eu for ionization V -XI Ionization 5p 6 5p 5 5p 4 5p 3 5p 2 5p V 4f 5 5p 6 4f 6 5p 5 VI 4f 4 5p 6 4f 5 5p 5 VII 4f 3 5p 6 4f 4 5p 5 * VIII 4f 3 5p 5 4f 4 5p 4 IX 4f 3 5p 4 4f 4 5p 3 X 4f 3 5p 3 4f 4 5p 2 XI 4f 3 5p 2 4f 4 5p The result for ground configurations are shown in bold for each ionization The configurations with asterix requires two configurations for optimization The energy levels are calculated for the configurations 4f 5 5p 6 and 4f 6 5p 5 The energy levels are calculated for the configurations 4f 3 5p 6 and 4f 4 5p 5  In this section, we explore the completeness of our new atomic data for calculating the opacity. For this purpose, we perform the convergence test on the opacity by using the atomic data corresponding to only a subset of configurations for IX ions of lanthanides ( Figure B1). We choose IX ion because this is one of the major contributors to the opacity at high temperature (∼ 70000 K) due to the highly dense energy levels. Figure B1 shows the distribution of the energy levels for individual configurations (left) and the Planck mean opacities calculated using a subset of atomic data (right). The opacity using the default configuration set, i.e., all the configurations mentioned in Table A4, is represented by the thick black curve.
The opacities for Eu IX remain unaffected as long as the energy levels belonging to the configurations up to 4f 4 5p 2 6p are included. This implies that the transitions to or from the energy levels belonging to 4f 4 5p 2 7s do not have a significant impact on the Planck mean opacity, mostly due to the negligible population in these relatively high energy levels (e.g., the energy levels belonging to 4f 4 5p 2 7s are > 80 eV, left panel of Figure B1). However, further removal of energy levels introduces about a factor of ∼ 2 difference in the opacity for Eu IX. Similar trends are found for Sm IX. On the other hand, for the case of Nd IX, the opacities are affected up to a factor of ∼ 2 if the energy levels belonging to the configurations up to 4f 5p 2 6p are not included. This is due to relatively low energy levels of the excited configurations in Nd. However, as the number density of the levels is not extremely high for Nd IX, the opacity itself is small, as compared with Eu IX and Sm IX. Therefore, we conclude that our opacities are mostly converged, and further addition of excited configurations will have a negligible effect on opacity.
C. SCHEME FOR REDUCED LINELIST For the highly ionized lanthanides, the number of lines is exceptionally high as described in Section 2.1. The numbers of transition for one ion reaches ∼ 0.1 billion in some cases (Table A4). With such a large linelist, performing radiative transfer simulation becomes infeasible. For the purpose of making the radiative transfer simulation possible, we create a reduced linelist by randomly choosing a single line out of n sample number of lines. If the reduced linelist preserves the statistical properties of the original linelist, i.e., if the statistical distributions of the transition wavelengths, radiative Up to 4f 4 5p 2 6p Up to 4f 4 5p 2 6s Up to 4f 4 5p 2 5d Figure B1. The distribution of the energy levels within the energy bin of 0.2 eV for the individual configurations as obtained from HULLAC (left panel) and the corresponding Planck mean opacities (right panel) for Nd IX, Sm IX, and Eu IX ions (top to the bottom panels). Different colors of the curves in right panel correspond to the opacities calculated by including different subsets of the configurations as shown in the legend.
transition probabilities, and the statistical weights of the energy levels are preserved, the opacity with the reduced linelist can approximately reproduce the original result, provided the contribution from the selected lines is enhanced by a factor of n sample (see Equation (C1)).