Signatures of r-process Elements in Kilonova Spectra

To investigate which elements produce strong absorption lines in the kilonova spectra, we systematically calculate the strength of bound–bound transitions for given density, temperature, and element abundances. The line strength is approximated by the Sobolev optical depth (Sobolev 1957) for each bound–bound transition,

$$\begin{eqnarray}&&{\tau }_{l}=\displaystyle \frac{\pi {e}^{2}}{{m}_{e}c}{f}_{l}{n}_{i,j}t{\lambda }_{l}\end{eqnarray} \tag{ 1 }$$

for homologously expanding ejecta. The Sobolev approximation is valid for expanding matter with a large radial velocity gradient. Here n_i,j is the number density of the lower level of the transition (ith ionized element in the jth excited state), and f_l and λ_l are the oscillator strength and the transition wavelength of the bound–bound transition, respectively. As in previous work on kilonovae (e.g., Barnes & Kasen 2013; Tanaka & Hotokezaka 2013), we assume local thermodynamic equilibrium (LTE); we assume a Boltzmann distribution for the population of excited levels and solve the Saha equation to obtain ionization states.

Atomic data are essentially important to evaluate the line strength. For theoretical calculations of kilonova light curves, atomic data from theoretical calculations have often been used (e.g., Kasen et al. 2013; Tanaka et al. 2018; Banerjee et al. 2020; Fontes et al. 2020). This is useful in terms of completeness of data because the opacity of the ejecta should be correctly evaluated for the light-curve calculations (see Appendix A). However, while such theoretical data may give a reasonable estimate for the total opacity, they are not accurate in transition wavelengths, and thus not suitable for element identification. In this work, since we focus on the imprints of elemental abundances in kilonova spectra, we construct the latest line list based on the Vienna Atomic Line Database (VALD; Piskunov et al. 1995; Kupka et al. 1999; Ryabchikova et al. 2015). This database is suitable to identify lines because the transition wavelengths are calibrated with experiments and semiempirical calculations. It should be, however, noted that the line list is not necessarily complete, in particular, in the NIR region, although kilonovae are brighter in the NIR at late times. The impact of the incompleteness is discussed in Section 4.3.

For the abundance in the ejected matter from an NS merger, we use a multicomponent free-expansion (mFE) model in Wanajo (2018). This model is constructed as an ensemble of the parameterized outflows with constant velocity, initial entropy, and initial Y_e, which fit the r-process residuals of the solar abundances (Prantzos et al. 2020). The ranges of velocity (in units of c), entropy (in units of Boltzmann constant per nucleon, k_B/nucleon), and Y_e are taken to be 0.05–0.30, 10–35, and 0.01–0.50 with intervals of 0.05, 5, and 0.01, respectively, as in Wanajo (2018). Here we consider three models (left panel of Figure 1): (1) a model that fits the r-process residuals for A ≥ 88, where A is mass number, i.e., including those heavier than the first r-process peak isotopes (A = 80–84); (2) a model that fits those for A ≥ 69 (including the first r-process peak isotopes) and 3% of those for A ≥ 100; and (3) a model that fits those for A ≥ 88 and 1% of those for A < 110. Hereafter we refer to these as the Solar (S), Light (L), and Heavy (H) models, respectively. Note that the S model is the same as mFE-b in Wanajo (2018), but the r-process residuals have been updated to those in Prantzos et al. (2020). The L model exhibits a similar abundance pattern to that of a metal-poor star with weak r-process signature (e.g., HD 122563, Honda et al. 2006; right panel of Figure 1). The H model represents a putative case, e.g., with a contribution from only dynamical ejecta. Note that the minimum mass number A = 88 for the S and H models corresponds to the dominant isotope of Sr, the element that has been measured in all r-process-enhanced stars (e.g., Cowan et al. 2021). We also performed the same calculations with the minimum mass number replaced by A = 85 (excluding the first r-peak and lighter isotopes) and 90 and confirmed that our results are unaffected by this choice.

Zoom In Zoom Out Reset image size

Although these models include the abundances with Z = 1–110, we use only those with Z = 20–100 in our calculations at t = 1.5 days as shown in the right panel of Figure 1. The heaviest elements with Z ≥ 101 are excluded because their mass fractions are very small (∼10⁻⁶ to 10⁻⁴), and there are no atomic data for such heavy elements. The light elements with Z ≤ 19 are also excluded because their mass fractions are also small, on order of 10⁻⁴, and they do not affect our results (see also Perego et al. 2020, for the effects of lightest elements). We summarize the mass fractions of selected elements for our models in Table 1. The distributions of Y_e for these models are shown in Figure 2 (see also Appendix B for the distributions of velocity and entropy). As can be anticipated from the abundance patterns (Figure 1), the distributions for the L and H models are dominated by higher (>0.3) and lower (<0.3) values of Y_e than those for the S model.

Zoom In Zoom Out Reset image size

The strength of bound–bound transitions at t = 1.5 days is displayed in Figure 3. As typical values in the ejecta (see Section 3), we evaluate the Sobolev optical depth for the density of ρ = 10⁻¹⁴ g cm⁻³ and temperature of T = 5000 K (top panels) and 7000 K (bottom panels).

Zoom In Zoom Out Reset image size

We find that the Sr ii triplet lines (red lines at 10000–11000 Å) are strong for the L model in the case of T = 5000 K. Other notable features are the Ca ii triplet lines (blue lines around 8500 Å), which are stronger than the Sr ii lines. Note that Ca and Sr abundances in this model are dominated by ⁴⁸Ca and ⁸⁸Sr (see the left panel of Figure 1). For higher temperature (T = 7000 K; bottom panel), these lines become significantly weaker.

In the S model, which has more heavy elements, either of the Sr ii (red) or Ca ii (blue) lines are not necessarily as strong as other lines of heavy elements, such as Tb iii and Th iii, when assuming T = 5000 K. For higher temperature (T = 7000 K), the Sr ii and Ca ii lines become significantly weaker as in the L model, while the strength of the Tb iii and Th iii lines is not largely affected. The results for the H model are similar to those for the S model.

The behavior of the line strength can be understood as the dependence of the Sobolev optical depths on temperature and density, as well as abundance distribution. Figure 4 shows the Sobolev optical depth for ranges of the density ρ = 10⁻²⁰ to 10⁻¹¹ g cm⁻³ and temperature T = 1000–10,000 K. For the Ca ii and Sr ii lines, the strength of the middle line of each triplet (8452.088 Å for Ca ii and 10327.31 Å for Sr ii) is shown. Regarding Ce iii lines, we choose one of the lines (9872.344 Å) as representative.

Zoom In Zoom Out Reset image size

When we focus on the case with the density of ρ = 10⁻¹⁴ g cm⁻³ as in Figure 3, the strength of the Ca ii and Sr ii lines is maximal around T ∼ 4000 K. Then, the lines become weaker for higher temperature as these elements are ionized to the doubly ionized state. This is the reason why the models with T = 5000 K show the stronger Sr ii and Ca ii lines compared to those with T = 7000 K. On the other hand, the strength of the Ce iii, Tb iii, and Th iii lines is maximal around T ∼ 6000 K, and thus they do not show a large difference between T = 5000 and 7000 K. As a result, the line strength of these heavy elements becomes much more dominant than that of the Sr and Ca lines for T = 7000 K for the S and H models, in which these heavy elements are abundant.

It is natural that not only Sr ii but also Ca ii triplet lines become strong, because these elements have similar atomic structures. Figure 5 shows the energy levels and the transition wavelengths for these ions (NIST ASD; Kramida et al. 2019). Both elements belong to group 2 in the periodic table and have only one electron in the outermost shell when they are singly ionized. The electron occupies the s shell in the ground state and has a relatively small number of excited levels. Therefore, the transition probability of each bound–bound transition tends to be high.

Zoom In Zoom Out Reset image size

We confirm that the observed absorption feature around λ ∼ 8000 Å in AT2017gfo is consistent with the Sr ii triplet. Figure 6 shows the comparison between the spectrum of AT2017gfo at t = 1.5 days (Pian et al. 2017) and the result of the L model (ρ = 10⁻¹⁴ g cm⁻³, T = 5000 K). In this model, the wavelengths of the transitions are blueshifted owing to the adopted expansion velocity of v = 0.2 c, which is suggested by the observations of AT2017gfo (e.g., Pian et al. 2017; Smartt et al. 2017). With this velocity, the wavelengths of the Sr ii triplet (red) lines match the observed broad absorption feature, as also demonstrated by Watson et al. (2019). Note that Watson et al. (2019) included elements with Z ≥ 33 in their calculations, i.e., Ca (Z = 20) was not included. Although our L model suggests that Ca ii can also exhibit strong absorption lines, the spectrum of AT2017gfo does not show such a feature at λ ∼ 6800 Å, an expected wavelength of the Ca ii triplet with v = 0.2 c. Implications of the absence of the Ca ii lines in the observed spectrum are discussed in Section 4.2.

Zoom In Zoom Out Reset image size

In this section, we calculate realistic synthetic spectra of kilonovae. The analysis in the previous section has evaluated the line strength for one-zone models. However, the ejecta from NS mergers have spatial distribution in density and temperature. Also, the temperature structure is controlled by radioactive heating and radiative transfer. Therefore, we perform radiative transfer simulations by taking the ejecta structure and radioactive heating into account. We use a wavelength-dependent radiative transfer simulation code (Tanaka & Hotokezaka 2013; Tanaka et al. 2014, 2017, 2018; Kawaguchi et al. 2018). The photon transfer is calculated by the Monte Carlo method. To compute the opacity for bound–bound transitions, we adopt the expansion opacity (Karp et al. 1977) and use the formula from Eastman & Pinto (1993):

$$\begin{eqnarray}&&{\kappa }_{\exp }(\lambda )=\displaystyle \frac{1}{{ct}\rho }\displaystyle \sum _{l}\displaystyle \frac{{\lambda }_{l}}{{\rm{\Delta }}\lambda }(1-{{\rm{e}}}^{-{\tau }_{l}}),\end{eqnarray} \tag{ 2 }$$

where τ_l is the Sobolev optical depth (Equation (1)). In the equation, the summation is taken over all the transitions within the wavelength bin Δλ. The Sobolev optical depths are evaluated in the same manner as in Section 2: we use the line list from the VALD database and assume LTE for ionization and excitation. Hereafter, we focus on the spectra at 1.5 days after the merger because insufficient atomic data in the NIR region start to affect the emission at later phases (see Section 4.3 and Appendix A).

For the ejecta density structure, we assume a single power law (ρ ∝ r⁻³) for the velocity range of the ejecta v = 0.05–0.3 c. The total ejecta mass is set to be ${M}_{\mathrm{ej}}=0.03{M}_{\odot }$. We use the three models for the abundance distributions as described in Section 2.1. The heating rate of radioactive nuclei is consistently taken from each model. The thermalization efficiency of γ-rays and radioactive particles follows the analytic formula by estimating characteristic timescales (Barnes et al. 2016).

Figure 7 shows our results of synthetic spectra at t = 1.5 days after the merger. To show the contribution of different elements, we also plot the Sobolev optical depths in the ejecta at v = 0.2 c. The wavelengths of lines are blueshifted according to v = 0.2 c. This choice broadly captures the absorption features of the synthetic spectra, which indicates that the line-forming regions in the ejecta are formed around this velocity for our models.

Zoom In Zoom Out Reset image size

The spectrum of AT2017gfo at t = 1.5 days (gray curve) is also plotted, which is vertically shifted for comparison. Our synthetic spectra do not necessarily reproduce the overall spectral shape: the L model underproduces the NIR flux, while the S and H models underproduce the optical flux. This is probably due to the assumption of the homogeneous abundance distribution in our models. It appears that a multicomponent model is needed to reasonably reproduce the observed spectrum of AT2017gfo (e.g., Kawaguchi et al. 2018). In this paper, therefore, we mainly focus on spectral features made by some elements.

We find that Sr ii and Ca ii produce absorption lines at λ ∼ 8000 Å and λ ∼ 6500 Å, respectively, in the spectrum of the L model (arrows in the top panel of Figure 7). This is reasonable since these triplet lines are stronger than other lines as shown in Section 2.2. The Sr ii absorption feature at λ ∼ 8000 Å is consistent with that seen in AT2017gfo, which is blueshifted according to v = 0.2 c. This also supports the identification of this feature attributed to the Sr ii lines (Watson et al. 2019).

For the S and H models, however, the Ca ii and Sr ii line features are overwhelmed by another broad absorption at around λ ∼ 7500 Å (arrows in the middle and bottom panels of Figure 7), which is mainly caused by Ce iii (Z = 58), Tb iii (Z = 65), and Th iii (Z = 90). The highest contribution comes from the Tb iii and Th iii lines (the rest-frame wavelength is 9115.425 Å for Tb iii and 9910.075 Å for Th iii), while the Ce iii lines give a minor contribution. This is mainly due to the large fraction of heavy elements in the S and H models. Additional effects by the temperature are discussed in Section 4.1.

Our results confirm that the absorption feature at λ ∼ 8000 Å in the spectrum of AT2017gfo can be caused by the Sr ii triplet, as reported by Watson et al. (2019). However, it is emphasized that the Sr ii lines are not always strongest: they give the strong absorption feature only in our lanthanide-poor (L) model but not in our lanthanide-rich (S or H) models. This suggests that the line-forming region of AT2017gfo at t = 1.5 days was lanthanide-poor (as discussed below).

Note that, in Watson et al. (2019), the absorption due to the Sr ii triplet is strongest at λ ∼ 8000 Å regardless of the existence of heavy elements. As the strength of the Sr ii lines is sensitive to the temperature (see Figure 4), this difference seems to be caused by the treatment of temperature in the ejecta. Our simulations take the structure of density and temperature in the ejecta into account and calculate the temperature from the radioactive heating rate in a consistent way. The temperature in the ejecta at v = 0.2 c is T ∼ 5200 K for the L model and T > 7000 K for the S and H models. This difference is mainly caused by the higher heating rate for the latter models. Since these models synthesize heavier elements than those in the L model, the average radioactive heating rate is higher (Wanajo 2018). Also, when the ejecta have more heavy elements, the overall opacity becomes higher, which suppresses the radiative cooling (Kasen et al. 2013; Tanaka & Hotokezaka 2013; Tanaka et al. 2018). Therefore, in the S or H models, the temperature in the ejecta is higher than that in the L model for a fixed velocity.

To see the effect of abundance distribution and temperature on the line strength, we show the density and temperature structures in the ejecta in Figure 8 and the Sobolev optical depths for each model as a function of velocity in Figure 9. When we focus on the Sr ii triplet lines, it is clear that they are strong in the L model. On the other hand, these lines are not necessarily strongest in the S and H models at the same velocity. Since the mass fraction of Sr is almost the same between the L and S models (see Table 1), the strength of Sr ii lines is mainly affected by the temperature difference (see also Figure 4). Also, provided that the temperature is lower (T ∼ 5000 K) in the S model, the Sr ii line strength becomes comparable to those of heavy elements, as shown in Figure 3. Thus, the feature by the Sr ii lines does not appear in the spectra of lanthanide-rich ejecta.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the line strength of heavy elements, the abundance distribution plays a major role. As shown in Figure 3, the strength of Ce iii, Tb iii, and Th iii lines is not largely changed between T = 5000 and 7000 K, and they are comparable to or even stronger than other lines. Therefore, if the ejecta have a large amount of heavy elements as in the S model, the spectra show the feature of these heavy elements.

The optical depths in the S and H models are qualitatively similar, and thus the synthetic spectra of these models also become similar. The strength of Ca ii and Sr ii lines for the H model is weaker than that for the S model because of their smaller abundances (bottom panel of Figure 9). Note that a significant difference between the spectra for the S and H models can appear at the later phase when the density and temperature in the ejecta decrease (see also Section 4.3).

Our results suggest that the absorption features depend on both temperature and abundance distribution in the ejecta. This means that different spectral features can appear in future kilonovae depending on, e.g., the ejecta masses and abundance distribution. In fact, our S and H models show spectra with absorption features by heavy elements such as Ce, Tb, and Th. The lines of these doubly ionized elements have been identified in chemically peculiar stars by using theoretical gf values (e.g., Ce iii, Ryabchikova et al. 2006; Th iii, Ryabchikova et al. 2007). Identification of such elements provides the direct evidence that heavy r-process elements are indeed synthesized in NS merger ejecta.

It should be noted that the transition probabilities of the doubly ionized ions of heavy elements are rather uncertain. The transition probabilities of these lines have been obtained only by theoretical calculations (Ce iii, Biémont et al. 2002b; Th iii, Biémont et al. 2002a; Tb iii, Database on Rare Earths At Mons University, DREAM; Biémont et al. 1999), while those for singly ionized Ca and Sr are experimentally evaluated (e.g., Ca ii, Theodosiou 1989; Sr ii, Kurucz 2017). Since kilonova ejecta are expected to have a wide range of temperature and ionization states, more experimental and observational calibrations for transition probabilities are necessary (see also Section 4.3).

It is also cautioned that our models are calculated by assuming one-dimensional ejecta with homogeneous abundance distribution. It is not clear how multidimensional modeling will affect our results presented here, depending on the spatial distributions of, e.g., ejecta masses, velocities, and abundances. We leave such exploration to future work.

Our results suggest that the Sr ii and Ca ii triplet lines can be used as a tracer of high-Y_e ejecta. The L model shows a strong absorption feature of the Ca ii lines, while the spectra of AT2017gfo do not exhibit such features. This implies that the abundance of Ca in the line-forming region of AT2017gfo is smaller than that in our L model. We perform the same radiative transfer simulations for the L model by varying the mass fraction of Ca to reconcile our model with the case for AT2017gfo. The result indicates that the Ca ii lines disappear when the ratio X(Ca)/X(Sr) is smaller than ∼0.002.

To infer the relevant physical properties from this constraint, we show the mass fractions of Ca and Sr as a function of Y_e for the individual outflows (Section 2.1) with different velocities and entropies in Figure 10. We first focus on the range of ${Y}_{{\rm{e}}}=0.30\mbox{--}0.50$ as a relevant condition in the L model (Figure 2) and the velocity of only v = 0.2 c, which is consistent with the blueshift of the Sr ii triplet in AT2017gfo. For these conditions, the Sr abundance (dominated by ⁸⁸Sr) is always relatively high (color lines in the middle panel), which is produced both in nuclear equilibrium and by neutron capture over a wide range of Y_e. By contrast, Ca, which is dominated by ⁴⁸Ca for ${Y}_{{\rm{e}}}\lt 0.45$, is preferentially produced in nuclear equilibrium with Y_e near its proton-to-nucleon ratio of 20/48 = 0.417 (Meyer et al. 1996; Wanajo et al. 2013). For this reason, the condition of X(Ca)/X(Sr) < 0.002 can be achieved only with the relatively high entropies of s ≥ 25 k_B/nucleon (bottom panel), which inhibit nuclear equilibrium from being established. This implies that the blue component of AT2017gfo, which produces the Sr absorption line, comes from relatively high entropy ejecta. Such a condition may be realized in the viscously heated, post-merger ejecta (e.g., Just et al. 2015; Lippuner et al. 2017; Fujibayashi et al. 2020a) or in the shocked dynamical ejecta along the polar direction (e.g., Wanajo et al. 2014; Sekiguchi et al. 2015; Radice et al. 2018).

Zoom In Zoom Out Reset image size

It should be noted that the conditions discussed above (v ∼ 0.2 c and s ≥25 k_B/nucleon) differ from the velocity and entropy distributions for our L model (the top panels of Figure 14 in Appendix B). However, this does not substantially affect our results shown in Sections 3 and 4.1 except for the Ca features. As cautioned in Wanajo (2018), our model obtained by fitting a given reference abundance distribution is not necessarily the unique solution. We show the final abundances for the individual outflows with the fixed velocity of ∼0.2 c and the entropy of ∼10–25 k_B/nucleon in Figure 11. It is clear that a similar fit to that for the L model can be obtained with those parameters, except for the range around Z = 20–40. The distributions of Y_e for those entropies are similar to those for the L model, which tend to, however, be slightly shifted to the higher-Y_e side for a higher entropy.

Zoom In Zoom Out Reset image size

Note also that we assume a solar r-process-like abundance distribution over a given range of atomic mass numbers to obtain our models (Section 2.1). If we relaxed this constraint, the ratio of X(Ca)/X(Sr) <0.002 would also be obtained with s ∼10 k_B/nucleon and ${Y}_{{\rm{e}}}\lesssim 0.35$ (Figure 10, bottom panel).

Our results imply that the Ca triplet line can be observed in future NS merger events. In particular,the post-merger ejecta with relatively high Y_e and relatively low entropy may have high ⁴⁸Ca abundance (Fujibayashi et al. 2020a). For such a case, we expect the Ca absorption feature with a moderate blueshift since the post-merger ejecta have relatively low expansion velocity (v ∼ 0.05 c). The detection of the Ca line in a future event will provide unique information on the ejecta condition as examined above for AT2017gfo.

Finally, it is interesting to note that ⁴⁸Ca is one of the isotopes whose origin remains a mystery. This isotope is synthesized most efficiently in the conditions of ${Y}_{{\rm{e}}}\sim 0.40\mbox{--}0.42$ and ≲15 k_B/nucleon (Meyer et al. 1996; Wanajo et al. 2013). Such a condition may be achieved in high-density Type Ia supernovae (Woosley 1997), electron-capture supernovae (Wanajo et al. 2013; Jones et al. 2019), low-mass core-collapse supernovae (Wanajo et al. 2018), collapsars (Fujibayashi et al. 2020b), or NS mergers (Wanajo 2018; Fujibayashi et al. 2020a). If the absorption line caused by the Ca ii triplet appears in future kilonova spectra, we will be able to confirm the production of ⁴⁸Ca by NS mergers.

The spectra of AT2017gfo show absorption features in the NIR region (e.g., Pian et al. 2017). We have performed systematic calculations and radiative transfer simulations including the NIR lines, but it is difficult to unambiguously identify the absorption features. This is because the transition data are not sufficient in the NIR wavelengths compared with those in the UV to optical wavelengths. Atomic data for the NIR transitions have been limited mainly owing to a lack of high-resolution stellar spectroscopic observations and laboratory experiments in the NIR wavelengths.

Recently, new measurements have become available owing to the development of NIR spectroscopy such as the APOGEE survey (Allende Prieto et al. 2008). When the abundances in certain Galactic stars have been obtained from their optical spectra, one can estimate the transition probabilities for NIR lines by measuring these stars in the NIR wavelengths. In fact, two such analyses about Ce ii (Z = 58) and Nd ii (Z = 60) have been performed for the APOGEE survey, which give the transition probabilities for the 9 lines of Ce ii (Cunha et al. 2017) and 10 lines of Nd ii (Hasselquist et al. 2016) in H band.

We perform here the same simulations as in Section 3, but including these new lines. The synthetic spectrum for the S model with the new Ce ii and Nd ii lines is shown in Figure 12. The result shows an additional absorption feature at λ ∼ 13000 Å, which is caused by the new Ce ii lines. These Ce ii lines are strong at around T ∼ 4000–5000 K, similar to the Sr ii lines. The H model also shows the same absorption feature, while the L model does not show it because of a small fraction of lanthanides (Table 1). Due to incompleteness of the line list in the NIR wavelengths, there is always a possibility that other unknown lines show stronger absorption lines. Therefore, it is not evident whether the Ce ii lines can be observed in future (S- or H-model-like) events. Our analysis, however, demonstrates that atomic data in the NIR wavelengths are important to fully decode kilonova spectra, in particular at the later phase when kilonovae are brighter in NIR.

Zoom In Zoom Out Reset image size