Ejected material from neutron star mergers gives rise to electromagnetic emission powered by radioactive decays of r-process nuclei, the so-called kilonova or macronova. While properties of the emission are largely affected by opacities in the ejected material, available atomic data for r-process elements are still limited. We perform atomic structure calculations for r-process elements: Se (Z = 34), Ru (Z = 44), Te (Z = 52), Ba (Z = 56), Nd (Z = 60), and Er (Z = 68). We confirm that the opacities from bound–bound transitions of open f-shell, lanthanide elements (Nd and Er) are higher than those of the other elements over a wide wavelength range. The opacities of open s-shell (Ba), p-shell (Se and Te), and d-shell (Ru) elements are lower than those of open f-shell elements, and their transitions are concentrated in the ultraviolet and optical wavelengths. We show that the optical brightness can be different by mag depending on the element abundances in the ejecta such that post-merger, lanthanide-free ejecta produce brighter and bluer optical emission. Such blue emission from post-merger ejecta can be observed from the polar directions if the mass of the preceding dynamical ejecta in these regions is small. For the ejecta mass of 0.01 , observed magnitudes of the blue emission will reach 21.0 mag (100 Mpc) and 22.5 mag (200 Mpc) in the g and r bands within a few days after the merger, which are detectable with 1 m or 2 m class telescopes.
The direct detection of gravitational waves (GWs) opened the era of GW astronomy (Abbott et al. 2016a, 2016c, 2017). A next important step will be the identification of their electromagnetic (EM) counterparts to further study the astrophysical nature of the GW sources (Nissanke et al. 2013), as sky localization by GW detectors is not accurate enough to pin down their positions (Abbott et al. 2016d). In fact, extensive EM follow-up observations have been performed for the detected GW events so far (Abbott et al. 2016b).
From compact binary mergers including at least one neutron star (NS), i.e., NS–NS mergers and black hole (BH)–NS mergers, various EM signals are expected over a wide wavelength range (e.g., Metzger & Berger 2012; Rosswog 2015). One of the most promising EM transients is the so-called "kilonova" or "macronova," which is the emission powered by the radioactive decay of newly synthesized r-process nuclei (Li & Paczyński 1998; Kulkarni 2005; Metzger et al. 2010). For recent reviews of kilonova emission, see Fernández & Metzger (2016), Tanaka (2016), and Metzger (2017). Kilonova emission is a good candidate for optical and near-infrared follow-up observations after the detection of GWs (Cowperthwaite et al. 2016; Kasliwal et al. 2016; Morokuma et al. 2016; Smartt et al. 2016; Soares-Santos et al. 2016; Yoshida et al. 2017).
Kilonova emission from the dynamical ejecta () of NS mergers is likely to have a luminosity of with a timescale of about 1 week, which is expected to peak at red optical or near-infrared wavelengths (Barnes & Kasen 2013; Kasen et al. 2013; Tanaka & Hotokezaka 2013; but see also Fontes et al. 2017; Wollaeger et al. 2017). This is due to the high opacities of r-process elements in the ejecta, especially those of lanthanide elements (Kasen et al. 2013). In fact, the short GRB 130603B showed a near-infrared excess in the afterglow (Berger et al. 2013; Tanvir et al. 2013), which was interpreted as a kilonova signal (see also Hotokezaka et al. 2013b; Piran et al. 2014). In addition, possible kilonova candidates have been reported for GRB 060614 (Jin et al. 2015; Yang et al. 2015) and GRB 050709 (Jin et al. 2016).
After the dynamical mass ejection, NS–NS mergers and BH–NS mergers are expected to have further mass ejection by viscous heating (Dessart et al. 2009; Fernández & Metzger 2013; Fernández et al. 2015a, 2015b; Shibata et al. 2017) that originates from magnetohydrodynamic turbulence (Price & Rosswog 2006; Kiuchi et al. 2014, 2015; Giacomazzo et al. 2015; Ciolfi et al. 2017; Siegel & Metzger 2017), and subdominantly by neutrino heating (Wanajo & Janka 2012; Perego et al. 2014; Fujibayashi et al. 2017) and nuclear recombination (Fernández & Metzger 2013). These components are as a whole denoted as "post-merger" ejecta in this paper. The post-merger ejecta can consist of less neutron-rich material than the dynamical ejecta (Just et al. 2015; Martin et al. 2015; Wu et al. 2016; Lippuner et al. 2017) since neutrino absorption can make the ejected material less neutron rich or the electron fraction (number of protons per nucleon) higher. If the ejecta are free from lanthanide elements, the emission from post-merger ejecta can be brighter and bluer, which can be called "blue kilonova" (Metzger & Fernández 2014; Kasen et al. 2015). However, due to the lack of atomic data for r-process elements, previous studies assume opacities of Fe for lanthanide-free ejecta. To predict emission properties of kilonova, systematic atomic data for r-process elements are important (see Kasen et al. 2013; Fontes et al. 2017; Wollaeger et al. 2017).
In this paper, we newly perform atomic structure calculations for selected r-process elements. Using these data, we perform radiative transfer simulations and study the impact of element abundances to kilonova emission. In Section 2, we show methods and results of our atomic structure calculations. In Section 3, we calculate opacities with these atomic data and discuss the dependence on the elements. We then apply our data to radiative transfer simulations in Section 4 and show light curves of kilonova from dynamical and post-merger ejecta of NS mergers. Finally, we give a summary in Section 5.
2. Atomic Structure Calculations
We perform atomic structure calculations for Se (Z = 34), Ru (Z = 44), Te (Z = 52), Ba (Z = 56), Nd (Z = 60), and Er (Z = 68). These elements are selected to systematically study the opacities of elements with different open shells: Ba is an open s-shell element, Se and Te are open p-shell elements, Ru is an open d-shell element, and Nd and Er are open f-shell elements. We focus on neutral atom and singly and doubly ionized ions because these ionization states are most common in kilonovae at day after the merger (Kasen et al. 2013; Tanaka & Hotokezaka 2013).
In Figure 1, these elements are shown with three different abundance patterns in the ejecta of NS mergers. While relativistic simulations of NS mergers predict wide ranges of between 0.05 and 0.45, the detailed distributions depend on the NS masses and their ratios as well as the adopted nuclear equations of state (Sekiguchi et al. 2015, 2016). In this paper, we assume a flat mass distribution between and 0.40 as representative of dynamical ejecta. As shown in Figure 1 (orange line), the dynamical ejecta consist of a wide range of r-process elements from the first (Z = 34) to third (Z = 78) abundance peaks. Note that the mass fractions of lanthanide and actinide elements are smaller than those in the dynamical ejecta model by Wollaeger et al. (2017), where abundances are largely suppressed at and enhanced at due to the dominance of a low- component.
For the post-merger ejecta, the exact abundance patterns are subject to uncertainties because of the difficulties in taking into account neutrino irradiation (e.g., Perego et al. 2014; Fujibayashi et al. 2017), viscous-driven turbulent motion (e.g., Giacomazzo et al. 2015; Kiuchi et al. 2015; Ciolfi et al. 2017; Siegel & Metzger 2017), and also neutrino oscillation (e.g., Frensel et al. 2017). Therefore, we consider single models of 0.25 (green) and 0.30 (blue) for simplicity. The former represents a case that contains the second (Z = 52) abundance peak and a small amount of lanthanides (~ 5 × 10−3 in mass fraction). The latter is a lanthanide-free model without elements of . For all of the models, the nucleosynthesis abundances of each are taken from Wanajo et al. (2014).
For the atomic structure calculations, we use two different codes, HULLAC (Bar-Shalom et al. 2001) and GRASP2K (Jönsson et al. 2013). The HULLAC code, which employs a parametric potential method, is used to provide atomic data for many elements while the GRASP2K code, which enables more ab initio calculations based on the multiconfiguration Dirac–Hartree–Fock (MCDHF) method, is used to provide benchmark calculations for a few elements. Such benchmark calculations are important because systematic improvement of the accuracies is not always obtained with the HULLAC code especially when little data are available in the NIST Atomic Spectra Database (ASD; Kramida et al. 2015). By using these two codes, we also study the influence of the accuracies of atomic calculations on the opacities. Tables 1 and 2 summarize the list of ions for atomic structure calculations. In the following sections, we describe our methods to calculate the atomic structures and transition probabilities.
Table 1. Summary of HULLAC Calculations
|Ion||Configurations||Number of Levels||Number of Lines||Subset 1a||Subset 2b|
|Se i||, c, , ,||3076||973,168||2395||654|
|Se ii||, c, , ,||2181||511,911||1978||584|
|Se iii||, c, , ,||922||92,132||2286||882|
|Ru i||, , , ,||1545||250,476||49,181||20,350|
|Te i||, ,||329||14,482||410||348|
|Te ii||, ,||253||9167||705||569|
|Te iii||, ,||57||419||249||227|
|Nd i||, , , ,||31,358||70,366,259||12,365,070||2,804,079|
|Nd ii||, , , ,||6888||3,951,882||3,682,300||1,287,145|
|Nd iii||, , , ,||2252||458,161||303,021||136,248|
|Er i||, ,||10,535||9,247,777||443,566||129,713|
|Er ii||, , , ,||5333||2,432,665||1,713,258||489,383|
Notes. Configurations taken into account for the optimization of central-field potentials are indicated by bold letters (see the text).aNumber of lines whose lower level energy is , 10, 15 eV for neutral atom and singly and doubly ionized ions, respectively. bNumber of lines with in Subset 1. c stands for single orbitals in all nl shells with and l = 0 − .
HULLAC (Hebrew University Lawrence Livermore Atomic Code; Bar-Shalom et al. 2001) is an integrated code for calculating atomic structures and cross sections for the modeling of atomic processes in plasmas and emission spectra. The latest version (9-601k) of HULLAC is used in the present work to provide atomic data for Se i–iii, Ru i–iii, Te i–iii, Nd i–iii, and Er i–iii. In HULLAC, fully relativistic orbitals are used for calculations of atomic energy levels and radiative transition probabilities. The orbital functions are the solutions of the single electron Dirac equation with a local central-field potential U(r) which represents a nuclear field and a spherically averaged interaction with other electrons in atoms,
where and β are the 4 × 4 Dirac matrices, c the speed of light in atomic units, and an orbital energy associated with each principal quantum number n, azimuthal quantum number l, and angular quantum numbers j and m. N-electron configuration state functions (CSFs) are constructed by coupled anti-symmetrized products of the orbital functions using Racah algebra. The jj-coupling scheme is used in HULLAC.
The total Hamiltonian, which consists of the N-electron Dirac–Coulomb Hamiltonian,
the magnetic and retardation effects in the interelectronic interaction (the Breit interaction), and quantum electrodynamic (QED) energy corrections, is diagonalized with multi-CSFs (relativistic configuration interaction: RCI). In Equation (2), VN is the monopole part of the electron–nucleus Coulomb interaction. Atomic energy levels are obtained from eigenvalues of the total Hamiltonian. Atomic wave functions of the energy levels are expressed as linear combinations of CSFs. Electric-dipole transition probabilities between two energy levels are obtained from the transition moments of the wave functions in length (Babushkin) gauge.
In Table 1, the configuration set used in the diagonalization as well as the number of energy levels and transitions of each ion are summarized. The ground state configuration is indicated at the top. It is noted that the configuration set in the present calculations should be read as minimal. Extended sets of configurations are used for Se to get improved energy levels. For the calculations of Nd, we use configurations similar to those included in the calculations by Kasen et al. (2013) and Fontes et al. (2017). We additionally include the configurations of and for Nd i, and and for Nd iii. Note that we do not include , which has higher energy levels than and .
The optimal central-field potential is obtained such that energy levels of the ground state and a few excited states agree with those in the ASD. The electron charge distribution function of q electrons in an nl shell is expressed by the density of the Slater-type orbital as
where A is a normalization factor and α values represent the average radii of the Slater-type orbital. The central-field potential for this electron charge distribution and the nuclear charge distribution seen by an external electron are obtained from the Poisson equation with the boundary condition . Occupancy of each Slater-type orbital is naively chosen as the ground state configuration of the next higher charge state. The ground state configuration for each ion is as given in the ASD. Alternative occupancies will give different electron charge density distributions, which result in different central-field potentials. In some cases, such alternative occupancies are used to improve results. For Ru i, an alternative occupancy [Kr] gives deeper and quasi-degenerate and orbital energies, resulting in a better agreement with the energy levels of the ASD. Similarly, the alternative occupancies [Xe] and [Cd] are used for Nd ii and Er iii, respectively, in the present calculations.
The α values that minimize the first-order configuration average energies of the ground state and low-lying excited states are chosen. Such α values depend on excited state configurations added in the first-order energies to be minimized. We choose the excited state configurations by single and double substitutions of valence and sub-valence orbitals from the ground state configuration. The ground state for each ion as well as the excited state configurations taken into account for the energy minimization are indicated by bold letters in Table 1. Getting correct energy levels by this semi-empirical optimization takes less computational time with limited computational resources, although systematic improvement of the results without a benchmark is not always possible. The results of a few lowest excited energy levels deviate from those of the ASD by about 10% at most for Se and Te. However, we cannot obtain such close agreements for Ru, reflecting the complexity of the atomic structures with open d shells. Results of the energies for Nd ii–iii and Er ii–iii are shown in Figure 2 and discussed in Section 2.3.
GRASP2K (Jönsson et al. 2013) is used to provide atomic data for Ba ii–iii, Nd ii–iii, and Er ii–iii. GRASP2K is based on the MCDHF and RCI methods taking into account Breit and QED corrections (Grant 2007; Froese Fischer et al. 2016). Based on the Dirac–Coulomb Hamiltonian (Equation (2)), the atomic state functions (ASFs) are obtained as linear combinations of symmetry adapted CSFs. The CSFs are built from products of one-electron Dirac orbitals. Based on a weighted energy average of several states, the so-called extended optimal level (EOL) scheme (Dyall et al. 1989), both the radial parts of the Dirac orbitals and the expansion coefficients are optimized self-consistently in the relativistic self-consistent field procedure. In the present calculations, ASFs are obtained as expansions over jj-coupled CSFs. To provide the LSJ labeling system, the ASFs are transformed from a jj-coupled CSF basis into an LSJ-coupled CSF basis using the method provided by Gaigalas et al. (2017).
The MCDHF calculations are followed by RCI calculations, including the Breit interaction and leading QED effects. Note that, for Nd ii and Er ii, only MCDHF calculations are performed. Radiative transition data (transition probabilities, oscillator strengths) between two states built on different and independently optimized orbital sets are calculated by means of the biorthonormal transformation method (Olsen et al. 1995). For electric-dipole and quadrupole (E1 and E2) transitions, we use the Babushkin gauge as in the HULLAC calculations.
In Table 2, we give a summary of the MCDHF and RCI calculations for each ion. As a starting point, MCDHF calculations are performed in the EOL scheme for the states of the ground configuration. The wave functions from these calculations are taken as the initial one to calculate the even and odd states of multireference configurations. The set of orbitals belonging to these multireference configurations are referred to as Layer 0. After that, the even and odd states are calculated separately. Unless stated otherwise, in the present calculations the inactive core for each of ions is mentioned in Table 2. The CSF expansions for states of each parity are obtained by allowing single and double substitutions from the multireference configurations up to active orbital sets (see Table 2). The configuration space was increased step by step with increasing layer number. The orbitals of previous layers are held fixed, and only the orbitals of the new layer are allowed to vary.
Table 2. Summary of GRASP2K Calculations
|Ion||Inact.||Ground||Multireference Set||Active Set||Number of Levels|
Note. Summary of the MCDHF and RCI calculations indicating inactive core, multireference set, active set, number of calculated levels, and number of configuration state functions () in the final list for each parity.a The ground state and states of the configuration as the nh orbital were excluded from the computations for .
The included configurations for the Nd calculations are essentially the same as those in Kasen et al. (2013) and Fontes et al. (2017). Note that, as in the HULLAC calculations, we do not include for Nd iii, which has high energy levels. Also, the GRASP2K calculations do not include and , which are included in the HULLAC calculations. These differences do not have a big impact on the opacities as the energy levels from the configurations are rather high.
Figure 2 shows the derived lowest energy for each electron configuration of the Nd ii, Nd iii, Er ii, and Er iii ions. For Nd ii, both HULLAC and GRASP2K calculations show reasonable agreement with the data in the ASD (open squares). Our calculations provide the correct orders of energy levels, and the deviation from the energy in the ASD is about , except for the configuration. Overall agreement is slightly better than that obtained by Kasen et al. (2013) with the Autostructure code (Badnell 2011). For Nd iii, the energy of the excited level is available only for , and our calculations give excellent agreement except for the Layer 0 calculation of GRASP2K (blue points).
For the case of Er, the agreement with ASD is not as good as in the Nd ions, which reflects the complexity of the Er ions. For Er ii, both HULLAC and GRASP2K calculations do not give the correct order of energy: the derived energy of is too high. On the other hand, for Er iii, the orders of energy are reproduced well although the number of available data in the ASD is small. For the case of Er iii, GRASP2K results give a better agreement than HULLAC. In the next section, we discuss the influence of these results on the opacities.
We calculate the bound–bound opacities by using the energy levels and the transition probabilities obtained from our atomic structure calculations. For bound–bound opacities, we use the formalism of expansion opacity (Karp et al. 1977; Eastman & Pinto 1993; Kasen et al. 2006),
which was also adopted in previous studies of kilonova simulations (Barnes & Kasen 2013; Kasen et al. 2013; Tanaka & Hotokezaka 2013; Tanaka et al. 2014). Here, is the Sobolev optical depth for a transition, and it can be written as
for homologously expanding ejecta. Here, and fl are the wavelength and the oscillator strength of the transition, respectively, and nl is the number density of the lower level of the transition. The summation in the expansion opacity is taken for all of the transitions within a wavelength bin (). The number density of each ion is calculated under the assumption of local thermodynamic equilibrium, and the population of the excited levels is calculated by assuming the Boltzmann distribution.
We confirm that there is an overall trend that the bound–bound opacities of open f-shell (lanthanide) elements are higher than those of open s-shell, p-shell, and d-shell elements over a wide wavelength range. The opacities of open d-shell elements (Fe and Ru) are concentrated on the ultraviolet and optical wavelengths, and those of open s-shell (Ba) and p-shell (Se and Te) elements also have a similar trend with even lower opacities. Figure 3 shows the opacity of each element calculated with , T = 5000 K, and day after the merger. For comparison, Fe opacities calculated with Kurucz's line list (Kurucz & Bell 1995) are also shown. The opacities of Nd and Er (open f-shell) are much higher than that of Fe (open d-shell). The opacity of Ru (open d-shell) is similar to that of Fe, which demonstrates the similarity in the opacity for the elements with the same open shell. The same is true for an open p-shell; the opacities of Se and Te are found to be similar.
The opacities from the two atomic codes agree reasonably well. Figure 4 shows the line expansion opacities of Nd ii, Nd iii, Er ii, and Er iii. As expected from the good agreement in the energy levels (Figure 2), the opacities from HULLAC and GRASP2K are almost indistinguishable for Nd ii and Nd iii. We note that only the Layer 0 calculations with GRASP2K give Nd iii opacities lower than the more realistic (Layer 1 and Layer 2) calculations. This is because the Layer 0 calculations give higher energy levels (Figure 2), which reduce the contribution of bound–bound transitions involving excited levels for a given temperature.
For the Er ions, we find that two atomic codes give larger discrepancies in the energy levels compared with the cases of Nd ions. As in the case of Nd iii, the opacities of Er ii and Er iii from the HULLAC calculations are slightly smaller than those from the GRASP2K calculations because HULLAC calculations give slightly higher energy levels for the excited energy levels. However, the difference in the opacity is only up to a factor of about 2. Therefore, we conclude that the relatively simplified calculations with the HULLAC code gives opacities with sufficient accuracy for astronomical applications.
We also compare our results with those of Kasen et al. (2013) and Fontes et al. (2017), who use different atomic codes. The thin dashed line in Figure 3 shows the opacities of Nd from HULLAC calculations with , T = 4000 K, and day after the merger for direct comparison. Our results are broadly consistent with those presented by Kasen et al. (2013); the opacities from HULLAC calculations are slightly higher but the difference is only a factor of 1.5 or 2 near the peak. Also, our opacities are consistent with the expansion opacities shown by Fontes et al. (2017); our opacities are slightly lower, but the difference is within a factor of 2 in the optical and near-infrared wavelengths.
Finally, we calculate the opacities for a mixture of elements. We use the HULLAC results, which cover more elements and ionization states. Because we have atomic structure calculations for a small number of elements, we assume the same bound–bound transition properties for the elements with the same open shell (see Figure 1). For open f-shell elements, the former and latter halves are replaced with Nd and Er, respectively. For the heavy elements with , we use the data of Ru, Te, Nd, and Er again. For the elements with , we use Kurucz's line list (Kurucz & Bell 1995). We neglect the contribution of open s-shell elements because the total fraction of these elements are small in the ejecta (Figure 1) and the opacities are subdominant (Figure 3).
As a result of the high opacity of lanthanide elements, the opacities for the mixture of elements depend significantly on . Figure 5 shows the line expansion opacity for the element mixture in the dynamical ejecta () and high- ejecta ( and 0.30). If the ejecta is completely lanthanide free as in the case of , the line expansion opacity is smaller than that in the lanthanide-rich ejecta by a factor of near the middle of optical range (~5000 Å). However, a small inclusion of lanthanide elements dramatically enhances the opacities as shown in the case of . This demonstrates the importance of accurate determination in the merger simulations for the accurate prediction of kilonova signals.
4. Radiative Transfer Simulations
We perform radiative transfer simulations by using our new atomic data. We use a three-dimensional, time-dependent, wavelength-dependent Monte Carlo radiative transfer code (Tanaka & Hotokezaka 2013). The code takes into account electron scattering and bound–bound, bound–free, and free–free transitions as sources of opacity. In the previous version of the code (Tanaka & Hotokezaka 2013; Tanaka et al. 2014; Tanaka 2016), we use the VALD database (Piskunov et al. 1995; Ryabchikova et al. 1997; Kupka et al. 1999, 2000) for the bound–bound transitions, while in this paper we use our atomic data presented in Section 2 and treat element mixtures by using representative elements as described in Section 3. To save memory space in the computation, we use a subset of the line list including the transitions whose lower level energy is , 10, and 15 eV for neutral atoms and singly and doubly ionized ions, respectively (Subset 1 in Table 1), and whose oscillator strengths are (Subset 2 in Table 1). We confirm that the use of this subset does not significantly affect the calculated light curves and spectra.
4.1. Simple Models
To study the effect of the element abundances (or ) on the light curves, we calculate the light curves for the three different abundance patterns displayed in Figure 1. For ease in extracting the effect of opacities on the element abundances, we employ a simple model of NS merger ejecta, the parameters of which are set to be the same for three cases. The ejecta mass is taken to be . The density structure of the ejecta is assumed to be spherical with a power-law radial profile of from to (Metzger et al. 2010; Metzger 2017). Light curves with such a simple density structure reasonably agree with those with a realistic structure from numerical relativity simulations (Tanaka & Hotokezaka 2013). This model has a characteristic velocity of , where is the kinetic energy of the ejecta.
For the heating rate of radioactive decays, we adopt ( is the time after the merger in days), as in the simple models adopted in Wollaeger et al. (2017). This gives a reasonable agreement with the nucleosynthesis calculations for a wide range of (Korobkin et al. 2012; Wanajo et al. 2014). We assume that a thermalization factor (), that is, a fraction of the decay energy deposited into the ejecta, is time independent, and adopt , which is a typical value at a few days after the merger (Metzger et al. 2010; Barnes et al. 2016; Rosswog et al. 2017).
Figure 6 shows the bolometric light curves of the simple models. The overall properties of the bolometric light curves can be understood by the dependence on the opacities; the characteristic timescale becomes shorter and the luminosity becomes higher for smaller opacities. As a result, the models with and 0.30 have higher bolometric luminosities at the first few days.
The bolometric light curve with can be compared with that of a Nd-rich model by Barnes & Kasen (2013) with the same ejecta mass and characteristic velocity. The overall properties of the light curves are quite similar as expected from the nice agreement in the opacities as well as similar assumptions adopted in the radiative transfer calculations, e.g., a constant thermalization efficiency. Our light curve gives a lower luminosity by a factor of about 2, but this difference is simply due to the smaller thermalization efficiency assumed in our simulations.
For comparison, we also show the results with gray transfer simulations with and . It should be noted, however, that the gray opacities give a reasonable approximation only for the bolometric luminosity. The properties of multicolor light curves cannot be well described by the gray opacities as expected from the strong wavelength dependence of the opacities (Figure 5).
4.2. Realistic Models
We perform radiative transfer simulations for more realistic models of dynamical and high- post-merger ejecta. For the dynamical ejecta model, we use the density structure from the results of the numerical relativity simulations by Hotokezaka et al. (2013a). As a representative case, we use the APR4-1215 model, which is the merger of two NSs with the gravitational masses of and . The ejecta mass is and the characteristic velocity is .
The dynamical ejecta have a wide range of , which is composed of a low- tidally disrupted component and a high- component by shock heating and/or neutrino absorption (Wanajo et al. 2014; Goriely et al. 2015; Sekiguchi et al. 2015; Foucart et al. 2016; Lehner et al. 2016; Radice et al. 2016; Sekiguchi et al. 2016). Therefore, we approximate the calculated abundance patterns in Wanajo et al. (2014) by a flat distribution from 0.10 to 0.40 as described in Section 2 (orange line in Figure 1). We assume a spatially homogeneous element distribution for simplicity, although merger simulations suggest that the polar regions consist mainly of high- material (e.g., Sekiguchi et al. 2016).
For models of post-merger ejecta, which originate from several possible mechanisms such as viscosity, neutrino heating, and nuclear recombination, we keep using a simple spherical model with a power-law profile as in Section 4.1, instead of using the results of numerical simulations. We set as a representative case. Because the typical velocity of the post-merger ejecta tends to be lower than that of the dynamical ejecta, we set the velocity range from to , which gives a characteristic velocity of . The elemental abundance distributions in the ejecta depend on the details of mass ejection mechanisms as discussed in Section 2. In this paper, to study the effect of lanthanide-free ejecta, we assume relatively high , that is, and 0.30 as shown in Figure 1. The spatial distribution of the elements is assumed to be homogeneous as in the dynamical ejecta.
For both dynamical and post-merger ejecta models, we use the heating rates from nucleosynthesis calculations for the relevant in Wanajo et al. (2014). The radioactive decay and subsequent energy deposition from β-decay, α-decay, and fission are taken into account. For the β-decay, 45% and 20% of energy are assumed to be carried out by γ-rays and β particles, respectively. Then, the thermalization efficiencies of γ-rays, α particles, β particles, and fission fragments are independently evaluated using analytic descriptions in Barnes et al. (2016).
Figure 7 shows the bolometric light curves for all of the models. The overall properties of the bolometric light curves are similar to the simple models in Section 4.1. The luminosity of each model at a few days after the merger is, however, higher than that of the simple model because the thermalization factor is higher than 0.25 at the early time. At days after the merger, the luminosities decline faster than those in the simple models due to decreasing thermalization efficiency. This trend is most notable for the case of because it has the smallest amount of the second peak elements with (see Figure 1) that have dominant contributions to the radioactive heating (Metzger et al. 2010; Wanajo et al. 2014). At late times ( days) when the heating from fission becomes important (Wanajo et al. 2014; Hotokezaka et al. 2016), the dynamical ejecta model (of ) that contains the transuranic elements gives the highest luminosity.
Figures 8 and 9 show the multicolor light curves and spectra for three models, respectively. The spectra of the dynamical ejecta model (APR4-1215) are very red, which peaks in the near-infrared at days. On the other hand, the post-merger ejecta model with has a peak in the optical at days. As a result, the post-merger ejecta model with is much brighter than the dynamical ejecta model in the optical, especially in the u, g, and r bands.
The properties of the light curves of the post-merger ejecta model with are in between the other two models, as expected from the intermediate opacities. Therefore, this model has hybrid properties; the optical brightness is higher than that of the dynamical ejecta model, and the near-infrared brightness is not as faint as that of the post-merger ejecta with (Figure 9).
The multicolor light curves with can be compared with those of a Nd-rich model by Barnes & Kasen (2013). Both models give similar peak magnitudes in the blue optical (−10 to −12 mag), red optical (−12 to −14 mag), and near-infrared (−15 to −14 mag). Due to the time dependence of the thermalization efficiency, our light curves decline more quickly. A closer comparison can be done for the rest-frame J-band magnitude, which was observed in GRB 130603B (Berger et al. 2013; Tanvir et al. 2013). A model with and in Barnes et al. (2016), where a time-dependent thermalization is taken into account, reaches −13.5 mag in the J band at ~7 days. Our model gives a similar brightness at the same epoch. It is noted that our model gives a brighter emission at earlier epochs. This may reflect the difference in the density structure; a flatter distribution is assumed in the outer layer of our models. Such a detail should be refined in the future by using more realistic ejecta models.
Our results confirm the presence of a "blue kilonova" that was previously suggested based on the use of iron opacity for the light r-process elements (Metzger & Fernández 2014; Kasen et al. 2015). For 0.01 of lanthanide-free () ejecta, the optical brightness reaches the absolute magnitude of mag in the g and r bands within a few days after the merger. This corresponds to 21.0 mag and 22.5 mag at 100 Mpc and 200 Mpc, respectively. Thanks to the relatively blue color, this emission is detectable with 1 m class and 2 m class telescopes, respectively.
It should be noted that the observability of blue kilonovae from lanthanide-free post-merger ejecta depends on the properties of the preceding dynamical ejecta as discussed in Kasen et al. (2015). If lanthanide-rich dynamical ejecta are present in all directions, the blue kilonova emission is likely to be absorbed. However, recent relativistic simulations with neutrino interaction show that dynamical ejecta can have relatively high near the polar regions (see, e.g., Sekiguchi et al. 2015; Foucart et al. 2016; Radice et al. 2016). In such cases, the blue emission from the post-merger ejecta can be observable from the polar direction without being absorbed. To test this hypothesis, it is necessary to consistently model the dynamical and post-merger ejecta. It is also noted that our simulations cannot predict the emission within ~1 day after the merger due to the lack of atomic data of more ionized elements. Emission at such early times can peak at optical or even ultraviolet wavelengths (Metzger et al. 2015; Gottlieb et al. 2017), and therefore, it will also be a good target for follow-up observations, especially with small telescopes.
We have newly performed atomic structure calculations for Se (Z = 34), Ru (Z = 44), Te (Z = 52), Ba (Z = 56), Nd (Z = 60), and Er (Z = 68) to construct the atomic data for a wide range of r-process elements. By using two different atomic codes, we confirmed that the atomic structure calculations gave uncertainties in opacities by only a factor of up to about 2. We found that the opacities from the bound–bound transitions of open f-shell elements were the highest from the ultraviolet to near-infrared wavelengths, while those of open s-shell, p-shell, and d-shell elements were lower and concentrated in ultraviolet and optical wavelengths.
Using our new atomic data, we performed multiwavelength radiative transfer simulations to predict a possible variety of kilonova emission. We found that, even for the same ejecta mass, the optical brightness varied by mag depending on the distribution of elemental abundances. If the blue emission from the post-merger, lanthanide-free ejecta with 0.01 is observable without being absorbed by the preceding dynamical ejecta, the brightness will reach the absolute magnitude of mag in the g and r bands within a few days after the merger. This corresponds to 21.0 and 22.5 mag at 100 Mpc and 200 Mpc, which is detectable with 1 m class and 2 m class telescopes, respectively.
We thank Kenta Hotokezaka, Masaru Shibata, Nobuya Nishimura, Kenta Kiuchi, and Koutarou Kyutoku for providing results of simulations and fruitful discussion, and the anonymous referee for constructive comments. M.T. thanks the Institute for Nuclear Theory (INT) at the University of Washington for its hospitality and the Department of Energy for partial support during the completion of this work. M.T. also thanks Rodrigo Fernández, Brian Metzger, Daniel Kasen, and Gabriel Martinez-Pinedo for organizing the workshop and providing the nice research environment at INT. Numerical simulations presented in this paper were carried out with Cray XC30 at the Center for Computational Astrophysics, National Astronomical Observatory of Japan. Computations by G.G. were performed on resources at the High Performance Computing Center "HPC Sauletekis" of the Faculty of Physics at Vilnius University. This research was supported by the NINS program for cross-disciplinary science study, the Inoue Science Research Award from Inoue Foundation for Science, the RIKEN iTHES project, a post-K computer project (Priority issue No. 9) of MEXT, and Grants-in-Aid for Scientific Research from JSPS (15H00782, 15H02075, 15K05077, 16H02183, 16H06341, 16K17706, 17K05391, 26400237) and MEXT (17H06357, 17H06363).