Properties of Kilonovae from Dynamical and Post-merger Ejecta of Neutron Star Mergers

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 ${\varphi }_{{nljm}}$ 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,

$$\begin{eqnarray}&&[c{\boldsymbol{\alpha }}\cdot {\boldsymbol{p}}+(\beta -1){c}^{2}+U(r)]{\varphi }_{{nljm}}={\varepsilon }_{{nlj}}{\varphi }_{{nljm}},\end{eqnarray} \tag{ 1 }$$

where ${\boldsymbol{\alpha }}$ and β are the 4 × 4 Dirac matrices, c the speed of light in atomic units, and ${\varepsilon }_{{nlj}}$ 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,

$$\begin{eqnarray}&&{H}_{\mathrm{DC}}=\displaystyle \sum _{i=1}^{N}(c{{\boldsymbol{\alpha }}}_{i}\cdot {{\boldsymbol{p}}}_{i}+({\beta }_{i}-1){c}^{2}+{V}_{i}^{N})+\displaystyle \sum _{i\gt j}^{N}\displaystyle \frac{1}{{r}_{{ij}}},\end{eqnarray} \tag{ 2 }$$

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), V^N 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 $4{f}^{3}5{d}^{2}6s$ and $4{f}^{3}5{d}^{2}6p$ for Nd i, and $4{f}^{2}5d6p$ and $4{f}^{2}6s6p$ for Nd iii. Note that we do not include $4f5{d}^{2}6s$, which has higher energy levels than $4{f}^{2}5d6p$ and $4{f}^{2}6s6p$.

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

$$\begin{eqnarray}&&\rho (r)=-4\pi {r}^{2}{qA}{[{r}^{l+1}\exp (-\alpha r/2)]}^{2},\end{eqnarray} \tag{ 3 }$$

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 $Z\delta (r)$ seen by an external electron are obtained from the Poisson equation with the boundary condition $U(r){| }_{r\to \infty }=(Z-q)/r$. 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] $4{d}^{5}5{s}^{2}$ gives deeper and quasi-degenerate $4d$ and $5s$ orbital energies, resulting in a better agreement with the energy levels of the ASD. Similarly, the alternative occupancies [Xe] $4{f}^{3}6s$ and [Cd] $5{p}^{5}4{f}^{12}$ 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.

Zoom In Zoom Out Reset image size

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		${N}_{\mathrm{CSFs}}$
	core	conf.	Even	Odd		Even	Odd	Even	Odd
Ba ii	[Kr]	$4{d}^{10}5{s}^{2}5{p}^{6}6s$	$4{d}^{10}5{s}^{2}5{p}^{6}{ns}$	$4{d}^{10}5{s}^{2}5{p}^{6}{np}$	$\{{ns},{np},(n-1)d,$	75	60	838,672	614,880
			$4{d}^{10}5{s}^{2}5{p}^{6}(n-1)d$	$4{d}^{10}5{s}^{2}5{p}^{6}(n-2)f$	$(n-2)f,(n-1)g\}$
			$4{d}^{10}5{s}^{2}5{p}^{6}(n-1)g$
			$n=6\mbox{--}20$
Ba iii	[Kr]	$4{d}^{10}5{s}^{2}5{p}^{6}$	$4{d}^{10}5{s}^{2}5{p}^{6}$ a	$4{d}^{10}5{s}^{2}5{p}^{5}{ns}$	$\{{ns},{np},(n-1)d,$	409	504	70,067	71,388
			$4{d}^{10}5{s}^{2}5{p}^{5}{np}$	$4{d}^{10}5{s}^{2}5{p}^{5}(n-1)d$	$(n-2)f,(n-1)g$
			$4{d}^{10}5{s}^{2}5{p}^{5}(n-2)f$	$4{d}^{10}5{s}^{2}5{p}^{5}(n-1)g$	nha}
			$4{d}^{10}5{s}^{2}5{p}^{5}{nh}$ a
			$n=6\mbox{--}23$
Nd ii	[Xe]	$4{f}^{4}6s$	$4{f}^{4}6s$, $4{f}^{4}5d$	$4{f}^{3}5{d}^{2}$, $4{f}^{4}6p$	$\{7s,7p,6d,5f\}$	3890	2998	24,568	23,966
			$4{f}^{3}5d6p$, $4{f}^{3}6s6p$	$4{f}^{3}5d6s$
Nd iii	[Xe]	4f4	4f4, $4{f}^{3}6p$	$4{f}^{3}5d$, $4{f}^{3}6s$	$\{8s,8p,7d,6g,6h\}$	1020	468	173,816	114,621
			$4{f}^{2}5{d}^{2}$, $4{f}^{2}5d6s$
Er ii	[Xe]	$4{f}^{12}6s$	$4{f}^{12}6s$, $4{f}^{12}5d$	$4{f}^{11}5d6s$, $4{f}^{11}5{d}^{2}$	$\{7s,7p,6d,5f\}$	2836	2497	22,460	21,731
			$4{f}^{11}5d6p$, $4{f}^{11}6s6p$	$4{f}^{11}6{s}^{2}$, $4{f}^{12}6p$
Er iii	[Xe]	4f12	4f12, $4{f}^{11}6p$	$4{f}^{11}5d$, $4{f}^{11}6s$	$\{7s,7p,6d,5f\}$	255	468	503,824	842,643

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 (${N}_{\mathrm{CSFs}}$) in the final list for each parity.

^a The ground state and states of the configuration $4{d}^{10}5{s}^{2}5{p}^{5}{nh}$ as the nh orbital were excluded from the computations for $n=7\mbox{--}23$.

Download table as:  ASCIITypeset image

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 $4f5{d}^{2}6s$ for Nd iii, which has high energy levels. Also, the GRASP2K calculations do not include $4{f}^{2}5d6p$ and $4{f}^{2}6s6p$, 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 $\lesssim 30 \%$, except for the $4{f}^{3}6s6p$ 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 $4{f}^{3}5d$, 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 $4{f}^{11}6{s}^{2}$ 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.

To study the effect of the element abundances (or ${Y}_{{\rm{e}}}$) 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 ${M}_{\mathrm{ej}}=0.01\,{M}_{\odot }$. The density structure of the ejecta is assumed to be spherical with a power-law radial profile of $\rho \propto {r}^{-3}$ from $v=0.05c$ to $0.2c$ (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 ${v}_{\mathrm{ch}}=\sqrt{2{E}_{{\rm{K}}}/{M}_{\mathrm{ej}}}=0.1c$, where ${E}_{{\rm{K}}}$ is the kinetic energy of the ejecta.

For the heating rate of radioactive decays, we adopt $2\times {10}^{10}\,{t}_{{\rm{d}}}^{-1.3}$ $\mathrm{erg}\,{{\rm{s}}}^{-1}\ {{\rm{g}}}^{-1}$ (${t}_{{\rm{d}}}$ 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 ${Y}_{{\rm{e}}}$ (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 $\epsilon =0.25$, 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 ${Y}_{{\rm{e}}}=0.25$ and 0.30 have higher bolometric luminosities at the first few days.

Zoom In Zoom Out Reset image size

The bolometric light curve with ${Y}_{{\rm{e}}}=0.10\mbox{--}0.40$ 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 $\kappa =1.0$ and $10\,{\mathrm{cm}}^{2}\ {{\rm{g}}}^{-1}$. 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).

We perform radiative transfer simulations for more realistic models of dynamical and high-${Y}_{{\rm{e}}}$ 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 $1.2\,{M}_{\odot }$ and $1.5\,{M}_{\odot }$. The ejecta mass is ${M}_{\mathrm{ej}}\,=0.01\,{M}_{\odot }$ and the characteristic velocity is ${v}_{\mathrm{ch}}=0.24c$.

The dynamical ejecta have a wide range of ${Y}_{{\rm{e}}}$, which is composed of a low-${Y}_{{\rm{e}}}$ tidally disrupted component and a high-${Y}_{{\rm{e}}}$ 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 ${Y}_{{\rm{e}}}$ 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-${Y}_{{\rm{e}}}$ 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 ${M}_{\mathrm{ej}}=0.01\,{M}_{\odot }$ 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 $v=0.02c$ to $0.1c$, which gives a characteristic velocity of ${v}_{\mathrm{ch}}=0.05c$. 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 ${Y}_{{\rm{e}}}$, that is, ${Y}_{{\rm{e}}}=0.25$ 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 ${Y}_{{\rm{e}}}$ 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 $t\,\gtrsim \,5$ 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 ${Y}_{{\rm{e}}}=0.3$ because it has the smallest amount of the second peak elements with $Z\sim 50$ (see Figure 1) that have dominant contributions to the radioactive heating (Metzger et al. 2010; Wanajo et al. 2014). At late times ($\gt 14$ days) when the heating from fission becomes important (Wanajo et al. 2014; Hotokezaka et al. 2016), the dynamical ejecta model (of ${Y}_{{\rm{e}}}=0.1\mbox{--}0.4$) that contains the transuranic elements gives the highest luminosity.

Zoom In Zoom Out Reset image size

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 $t=1\mbox{--}20$ days. On the other hand, the post-merger ejecta model with ${Y}_{{\rm{e}}}=0.3$ has a peak in the optical at $t\lesssim 5$ days. As a result, the post-merger ejecta model with ${Y}_{{\rm{e}}}=0.3$ is much brighter than the dynamical ejecta model in the optical, especially in the u, g, and r bands.

Zoom In Zoom Out Reset image size

The properties of the light curves of the post-merger ejecta model with ${Y}_{{\rm{e}}}=0.25$ 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 ${Y}_{{\rm{e}}}=0.3$ (Figure 9).

Zoom In Zoom Out Reset image size

The multicolor light curves with ${Y}_{{\rm{e}}}=0.10\mbox{--}0.40$ 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 ${M}_{\mathrm{ej}}=0.01\,{M}_{\odot }$ and ${v}_{\mathrm{ch}}=0.1c$ 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 ${M}_{\odot }$ of lanthanide-free (${Y}_{{\rm{e}}}=0.3$) ejecta, the optical brightness reaches the absolute magnitude of $M=-14$ 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 ${Y}_{{\rm{e}}}$ 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.