Modeling Kilonova Light Curves: Dependence on Nuclear Inputs

The astrophysical parameters that describe the initial conditions of r-process outflows launched by compact object mergers giving rise to KNe have been studied extensively, for example, in Korobkin et al. (2012), Martin et al. (2015), Lippuner & Roberts (2015), Wanajo (2018), and Miller et al. (2019). The nuclei synthesized in NSM ejecta are determined by a number of factors, including the entropy, the gas rate of expansion and initial neutron richness of the ejecta. Here we choose astrophysical conditions typical of merger accretion disk winds and use variations in the initial electron fraction, Y_e , of the ejecta as a proxy for all astrophysical variations.

The base astrophysical trajectory is a standard parameterized wind (Panov & Janka 2009), with initial entropy per baryon in units of a Boltzmann constant of s/k = 40 and an expansion timescale of 20 ms, conditions similar to those found in the multidimensional NSM simulations of, e.g., Just et al. (2015), Radice et al. (2018a), and as applied in Zhu et al. (2018). Note that nucleosynthesis proceeds slightly differently at different entropy, with a higher entropy working similarly to a lower Y_e , as both can increase the r-process reach (Meyer & Brown 1997). Simulations begin in nuclear statistical equilibrium at a temperature of T = 10 GK with seed nuclear abundances generated using the SFHo equation of state (Steiner et al. 2013). For all simulations, we use the same evolution of density as a function of time (expansion rate) but choose a variety of initial electron fractions. With different initial electron fractions, nucleosynthesis processes and their corresponding energy generation as a function of time are different. Therefore, the temperature of each simulation must be determined separately; and we use an adiabatic expansion law that is modified assuming a 10% nuclear reheating efficiency as in Holmbeck et al. (2018) and Vassh et al. (2019). As the energy generation is different for each set of nuclear inputs, each simulation we consider has a unique thermodynamic history.

Based on an initial survey of Y_e from 0.01–0.30 (with increments of 0.01) we find that a subset of these is sufficient to capture most of the effects on the light curve. Therefore, to keep simulation expense low without losing a reasonable coverage of electron fractions, we use Y_e of 0.02, 0.12, 0.14, 0.16, 0.18, 0.21, 0.24, and 0.28 in our nucleosynthesis calculations, which we describe next.

For the nucleosynthesis simulations, we use the Portable Routines for Integrated nucleoSynthesis Modeling (PRISM) (T. M. Sprouse et al. 2020, in preparation) reaction network developed at the University of Notre Dame and Los Alamos National Laboratory. PRISM is designed to work with prepared input files of astrophysical conditions, nuclear mass models, and nuclear reaction channels (such as charged-particle reactions, neutron capture, photodissociation, β-decay, β-delayed neutron emission, neutron-induced fission, β-delayed fission, and spontaneous fission as in Zhu et al. 2018; Vassh et al. 2019). For experimentally measured nuclei, we use measured masses from AME2016 from Wang et al. (2017) and experimental data from NUBASE2016 from Audi et al. (2017), including α-decay, β-decay, electron capture, neutron emission, and spontaneous fission data. Since we are interested in uncertainties from theoretical nuclear inputs for nuclei heavier than iron, we use the JINA Reaclib nuclear reaction database (Cyburt et al. 2010) for charged-particle and light-nuclei reactions, but otherwise explore the effect of a wide range of predictions for as yet unmeasured masses and reaction rates for heavier nuclei.

Where experimental data is not available, we use nuclear physics inputs starting from the eight theoretical nuclear models listed in Table 1, and all the reaction rates are determined consistently with the associated masses (Mumpower et al. 2015b). The neutron capture and neutron-induced fission rates are recalculated with each chosen mass model using Los Alamos National Laboratory statistical Hauser–Feshbach code CoH (Kawano et al. 2016). The β-decay strength functions are from Möller et al. (2019), and the relative probabilities for β-decay, β-delayed neutron emission, and β-delayed fission are calculated using the QRPA+HF framework (Mumpower et al. 2016a, 2018). Theoretical α-decay rates are obtained with a Viola–Seaborg relation using Q-values calculated from each chosen mass model and parameters fit to known data (Vassh et al. 2019). For neutron-induced fission, β-delayed fission, and spontaneous fission, we use either a symmetric (${A}_{\mathrm{fragment}}=\tfrac{1}{2}{A}_{\mathrm{fission}}$) or a double Gaussian fragment distribution from Kodama & Takahashi (1975); most modern fission yield predictions fall somewhere between these extremes. We make an exception for the fission yields for ²⁵⁴Cf, in which case we use the yield distribution described by Zhu et al. (2018). For spontaneous fission half-lives, we apply a barrier-height-dependent prescription from Karpov et al. (2012) and Zagrebaev et al. (2011) (KZ) or the parameterization from Xu & Ren (2005) (XR). Fission rates are updated with the fission barrier predictions most closely associated with a given mass model as in Vassh et al. (2019). That is, for the Extended Thomas–Fermi plus Strutinsky integral (ETFSI) model we apply ETFSI barriers, for HFB22 and HFB27 we apply HFB-14 barriers, for the Thomas–Fermi (TF) model we apply TF barriers, and for all other cases, including FRDM2012, we apply finite-range liquid drop model (FRLDM) barriers.

The total heat supplied to the light curve is a function of the energy emitted by radioactive decay and the efficiency with which that energy thermalizes, i.e., is converted into the thermal energy that ultimately forms the light curve. The distribution of radioactively decaying nuclei evolves rapidly in the days and weeks following the merger. We keep track of the energy emitted through each channel (β-decay, α-decay, and fission) as a function of time for each nucleosynthesis simulation. This is important because there is not a one-to-one correspondence between a final abundance pattern and the history of heating and nuclear decays (Kawaguchi et al. 2020).

Different decay channels have different overall levels of thermalization efficiency. Energy released in β-decays heats the ejecta less effectively than that released in α-decay, which in turn is less efficient than fission, as shown in Barnes et al. (2016). This is partly because a substantial fraction of β-decay energy is released as neutrinos, which free-stream out of the diffuse KN ejecta without interacting. In contrast, in α-decays and fission, most of the energy goes to massive particles, which thermalize more effectively in the ejecta. To capture these effects, we calculate the rate at which energy is produced by β-decay, α-decay, and fission as a function of time. We then apply analytic thermalization efficiencies from Kasen & Barnes (2019) (hereafter referred to as KB19) for each channel to estimate the total rate at which energy is thermalized, $\dot{Q}(t)$, as a function of time,

$$\begin{eqnarray}&&\dot{Q}(t)=\sum _{}{{if}}_{i}({M}_{\mathrm{ej}},{v}_{\mathrm{ej}},t)\ {\dot{q}}_{i}(t)\ {M}_{\mathrm{ej}}.\end{eqnarray} \tag{ 1 }$$

In the above, f_i is the thermalization efficiency for the reaction channel i, which depends on the mass, M_ej, and characteristic velocity, v_ej, of the ejecta in addition to time, and ${\dot{q}}_{i}(t)$ is the nuclear heating released by channel i per unit mass of ejecta, as determined from the reaction rates and Q-values of each nuclear reaction in each nucleosynthesis simulation. We determine the thermalization efficiencies, f_i (t), following the procedure outlined in KB19 in which the thermalization of massive particles is assumed to have the form f = (1 + τ)⁻ⁿ . Here, τ is the time, scaled to a characteristic "thermalization time," t_th, at which thermalization starts to become inefficient. Both t_th and the power-law index, n, depend on the parameters of the ejecta, the energetics of the decaying nuclei, the magnitude of the energy-loss cross-sections, and their dependence on particle energy.

The total thermalization efficiency for β-decay depends on the efficiencies for electrons and γ-rays, as well as the fraction of energy lost to neutrinos. We assume that 25% of the energy goes to electrons, 25% to γ-rays, and the remainder to neutrinos (Barnes et al. (2016); Hotokezaka et al. (2016)), which is consistent with the decays that are occurring in the calculations presented in the next sections. The electron efficiency follows the pattern defined above, with n = 1, and t_th defined by (KB19, Equation (43))

$$\begin{eqnarray}&&{t}_{\mathrm{th},\beta }=12.9\ {M}_{0.01}^{2/3}\ {v}_{0.2}^{-2}\ {\zeta }^{2/3}\ \mathrm{days},\end{eqnarray} \tag{ 2 }$$

where M_0.01 is the ejecta mass in units of 0.01 M_⊙, v_0.2 = v_ej/0.2c, and ζ is a constant defined to be close to unity (we adopt ζ = 1). For γ-rays, the thermalization efficiency is assumed to be comparable to the probability of absorption or scattering in the ejecta (KB19, Equations (48) and (53)),

$$\begin{eqnarray}&&{f}_{\gamma }(t)=1-\exp [-{t}_{\gamma }^{2}/{t}^{2}],\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{t}_{\gamma }=0.3\ {M}_{0.01}^{1/2}\ {v}_{0.2}^{-1}.\end{eqnarray} \tag{ 4 }$$

The efficiencies for α-decay and fission are the efficiencies of the corresponding particles. For α-decay,

$$\begin{eqnarray}&&{f}_{\alpha }{(t)=(1+t/{t}_{\mathrm{th},\alpha })}^{-1.5},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{t}_{\mathrm{th},\alpha }=2\ {t}_{\mathrm{th},\beta },\end{eqnarray} \tag{ 6 }$$

while for fission fragments,

$$\begin{eqnarray}&&{f}_{{\rm{f}}}{(t)=(1+t/{t}_{\mathrm{th},{\rm{f}}})}^{-1},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{t}_{\mathrm{th},{\rm{f}}}=4\ {t}_{\mathrm{th},\beta }.\end{eqnarray} \tag{ 8 }$$

The scalings and choice of n come from comparing the typical energies and energy-loss rates of the different particle types (Barnes et al. 2016). The thermalization efficiencies f_i (t) for each decay channel are plotted in Figure 1.

Zoom In Zoom Out Reset image size

We estimate the opacity of the KN ejecta from the composition of each simulation at a representative time t_κ  = 1 day after the merger. The assumption of constant composition is standard in radiation transport models of KNe (Kasen et al. 2017; Tanaka et al. 2020; Even et al. 2020, among many others), and is justified by the much slower evolution of the composition on these longer timescales relative to the first minutes and hours after free neutron depletion. In particular, the mass fractions of the highest opacity species are not found to change dramatically over the course of the KN, and the effects of temperature are expected in this epoch to have a greater effect on the opacity than any lingering compositional drift (Kasen et al. 2013; Tanaka et al. 2020).

Opacity has been shown to depend critically on bound-bound structure, and because we lack atomic data for all atomic species, we carry out our opacity calculation for a carefully constructed modified composition. We take the mass fraction in all s- and p-block elements and, for the purpose of calculating opacity only, assume that this mass resides in calcium (Z = 20). Since the opacities from d- and f-block dominate the total opacity of the gas, we modify the full composition in a manner that preserves the total mass fractions in d- and f-blocks. In the d-block, we evenly divide the mass fractions in d-block elements among Z = 21–28. In the f-block (lanthanides and actinides) we take care to preserve the mass fraction in each column of the periodic table, since atomic opacities vary both among and within the blocks of the periodic table (Kasen et al. 2013; Tanaka et al. 2020), and because f-block abundances have the strongest impact on the total opacity. In our modified composition, lanthanide and actinide mass fractions are distributed as described above among Z = 59–70.

We calculate the Planck mean of the expansion opacity (Karp et al. 1977) for the modified compositions described above at a range of temperatures, assuming local thermodynamic equilibrium and adopting a mass density ρ = 10⁻¹⁴ cm² g⁻¹, typical of conditions in the ejecta on KN timescales. The results of these calculations inform the temperature-dependent opacity function we use to compute semi-analytic light curves:

$$\begin{eqnarray}\kappa =\left\{\begin{array}{ll}{\kappa }_{\max }{\left(\tfrac{T}{4000{\rm{K}}}\right)}^{5.5}, & T\lt 4000\,{\rm{K}}\\ {\kappa }_{\max } & \mathrm{otherwise}.\end{array}\right.\end{eqnarray} \tag{ 9 }$$

In the equations above, κ_max is the maximum value from the Planck mean opacity calculation, and is unique for each composition. The steep decline in κ for T < 4000 K accounts for the opacity lost as the higher-energy electronic states of high-opacity lanthanide and actinide atoms and ions are depopulated due to the decreasing amount of available thermal energy in a cooling gas.

While opacity does have an effect on our calculations, since we are looking at lanthanide-rich scenarios, it tends to be less significant than other uncertainties such as the nuclear heating. We explore the effects of different lanthanide and actinide mass fractions, and therefore changing opacity, in Section 4.2.

The effective heating, $\dot{Q}(t)$ and opacity, κ(T(t)) are input into a semi-analytic light-curve model based on that of Metzger (2017), in which the KN ejecta is discretized as a set of concentric shells of mass M_v . We choose a power-law mass density profile

$$\begin{eqnarray}&&\rho (v,t)=\left(\displaystyle \frac{3}{4\pi }\right)\displaystyle \frac{{M}_{\mathrm{ej}}}{{v}^{3}{t}^{3}}{\left(\displaystyle \frac{v}{{v}_{0}}\right)}^{-3},\end{eqnarray} \tag{ 10 }$$

where v is the velocity coordinate, which extends from v₀ = 0.1c to v_max = 0.4c. Since we are interested in the red KN component containing the most massive r-process nuclei, we adopt as our fiducial model an ejecta with mass M_ej = 0.05 M_⊙ and v_ej = 0.15 c, typical of the values inferred for the outflow powering the red emission of the KN associated with GW170817 (Drout et al. 2017; Kasen et al. 2017; Kasliwal et al. 2017; Perego et al. 2017; Tanaka et al. 2017). The effect of varying ejecta mass is explored in Section 4. The characteristic ejecta velocity, defined in terms of the ejecta's total kinetic energy, E_kin, through ${v}_{\mathrm{ej}}={(2{E}_{\mathrm{kin}}/{M}_{\mathrm{ej}})}^{1/2}$, is used to calculate thermalization efficiencies as expressed in Equations (2) and (4). In our fiducial model, v_ej = 0.15 c. We consider 100 shells and evolve the thermal energy, E_v , of a shell of mass M_v at velocity v as

$$\begin{eqnarray}&&\displaystyle \frac{{{dE}}_{v}}{{dt}}=\displaystyle \frac{{M}_{v}}{{M}_{\mathrm{ej}}}\dot{Q}(t,v)-\displaystyle \frac{{E}_{v}}{t}-{L}_{v},\end{eqnarray} \tag{ 11 }$$

where the first term on the right-hand side is the effective heating from (thermalized) radioactivity, the second term accounts for loss of radiation energy due to adiabatic expansion, and the final term represents the outgoing luminosity from the mass shell. The luminosity comprises photons that diffuse and free-stream out of the shell, and is approximated as

$$\begin{eqnarray}&&{L}_{v}=\displaystyle \frac{{E}_{v}}{{t}_{{\rm{d}},{\rm{v}}}+{t}_{\mathrm{lc}}}.\end{eqnarray} \tag{ 12 }$$

Here, t_lc is the light-crossing time and t_d,v the diffusion timescale at velocity v. These quantities are given by

$$\begin{eqnarray}&&{t}_{\mathrm{lc}}={vt}/c,\ \mathrm{and}\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{t}_{{\rm{d}},{\rm{v}}}={M}_{\gt v}\kappa ({T}_{v}(t))/4\pi {vtc}.\end{eqnarray} \tag{ 14 }$$

In Equation (14), M_>v is the amount of mass exterior to the mass shell with velocity v, i.e., the amount of mass which has velocity greater than v. The opacity of a given layer, as described above, depends on the temperature, T_v , which we calculate from E_v assuming the ejecta is radiation dominated. The total luminosity at time t is the sum of the luminosity escaping from each layer.

In Figure 8, we show the overall effective heating rate for the 13 simulations, calculated using the thermalization efficiencies described in Section 2.3 with M_ej = 0.05 M_⊙ and v = 0.15 c. In this figure, the gray band represents the range of effective heating rates from all simulations with Y_e of 0.02–0.28, as well as the full set of nuclear physics inputs described in Section 2.2. The combined mass fraction of lanthanides and actinides, ${X}_{\mathrm{Ln}+\mathrm{An}}({t}_{\kappa }=1\ \mathrm{day})$, varies from 0.020–0.481, as shown in the corner subplot of Figure 8. In Figure 9, we show the fractional effective heating rates for the individual simulations 1–11. Meanwhile, Figure 10 shows the fractional effective heating rates for simulations 12 and 13, in addition to the individual simulations of which they are composed.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Prior to 1 day post-merger, we see that most simulations are dominated by β-decay. The exceptions to this are simulations 3 and 4, which use the Hartree–Fock–Bogoliubov (HFB) theoretical nuclear model. In these cases, Figure 9 shows the dominance of spontaneous fission at this time, and we see a correspondingly high effective heating rate, shown in Figure 8. Beginning at 8–10 days post-merger, we start to see α-decay and spontaneous fission playing an increasingly important role in the overall heating for simulations with Y_e  < 0.21, as well as those fit to the solar pattern. We note that these cases are weighted such that they include substantial contribution from simulations with Y_e  < 0.21 in order to produce solar actinide yields. Around this time, we see that the effective heating rates for individual simulations decay with slopes roughly proportional to t⁻¹ − t^−1.1. For comparison purposes, these timescales are shown by the red lines in Figure 8. Simulations 1, 8, and 10 (Y_e  ≥ 0.21) lack significant contributions from spontaneous fission and α-decay, and show lower heating rates than their lower Y_e counterparts. However, we also see that for a given Y_e , the effective heating rate is also sensitive to the theoretical nuclear model.

Around 100 days, simulations 1, 8, and 10 continue to be dominated by β-decay. The effective heating rate for these simulations decays rapidly at this time, roughly proportional to t^−2.3, as indicated by the solid red line in Figure 8. Meanwhile, simulations with low Y_e , as well as simulations 12 and 13, are dominated by spontaneous fission and in most cases show a "bump" in the overall heating compared to the power-law decays. This "extra" effective heating mainly comes from the spontaneous fission of ²⁵⁴Cf, which is consistent with the findings of Zhu et al. (2018). After hundreds of days, we tend to see a decreased importance of spontaneous fission and a corresponding rapid decay in the effective heating rate. We further explore the sources and roles of these different processes and timescales in energy generation in Section 5.

The relationship between r-process radioactivity and the emerging light curve is governed by multiple interacting factors. While the total effective heating sets the overall magnitude of the light curve, the maximum luminosity also depends on the time at which the light-curve peaks, which is a function of the ejecta opacity. Changes in the slope of the effective heating curve can also affect the shape of the light curve near its peak, while the decline of the light curve is a function of late-time effective heating, which is directly tied to nucleosynthetic outcomes, particularly the prevalence of ${\unicode{x00251}}$-decay and fission. We now take the 13 models given in Table 2 and calculate the bolometric luminosity for each of them as described in Section 2. Figure 11 shows the results of this calculation assuming an ejecta mass of M_ej = 0.05 M_⊙ and an ejecta velocity of v_ej = 0.15 c. For comparison, the observed bolometric luminosity from the KN associated with GW170817 (Drout et al. 2017) is also shown. This includes emission from a blue component, which dominates for the first ∼2 days, and a red component, which is the main contributor after roughly 4 days. Since our model is concerned only with the red component, the relevant time period for comparison with data is t ≳ 4 days. As can be seen in Figure 11, the peak bolometric luminosity is on the scale of 10⁴¹ erg s⁻¹, and all simulations reach their maximum values by times ranging from 2 days to about 12 days. We note that this behavior is for an ejecta mass close to that generally inferred from GW170817. We explore the effects of varying the ejecta mass in the next section.

Zoom In Zoom Out Reset image size

Much of the variability in the luminosity is a reflection of the variability in the effective nuclear heating seen in Figure 8. For the selected simulations, the range of uncertainty in the peak bolometric luminosity is from 5 × 10⁴⁰ to 3 × 10⁴¹ erg s⁻¹. The general trend is that a larger effective heating rate at the peak time leads to an overall more luminous light curve. For example, the effective heating rates of simulation 3 and 4 are close to the upper bound of the effective heating rate range shown in Figure 8 and the peak bolometric luminosities of simulations 3 and 4 shown in Figure 11 are also the most luminous. At a few days after the merger, simulations 9 and 10 have the lowest effective heating rate as shown in Figure 8 and their peak bolometric luminosities (brown lines in Figure 11) are also on the low side.

Another way in which differences in nuclear heating manifest themselves is in the shape of the luminosity curve at late times. With lower Y_e simulations, it is more likely that nucleosynthesis will proceed to heavy nuclei that can fission. As discussed in Sections 2 and 3, this boosts the effective heating rate, both because of the large Q-values involved in fission and also due to the enhanced thermalization of fission products. In contrast, the element synthesis in simulation 1 never reaches fissioning nuclei, and only β-decays contribute significantly to the heating. Therefore, the luminosity in simulation 1 (dark purple line in Figure 11) drops much more steeply at late times. We explore the landscape of decaying nuclei more closely in Section 5.

If the effective heating rate is the same, the higher opacity that stems from a larger ${X}_{\mathrm{Ln}+\mathrm{An}}$ results in an elongated light curve and a later, dimmer peak, as compared to a simulation with a smaller ${X}_{\mathrm{Ln}+\mathrm{An}}$. For example, simulations 1 and 12 have similar effective heating rates (dark purple and green lines in Figure 8) at several days, but the peak bolometric luminosity of simulation 1 (dark purple line in Figure 11) is noticeably earlier than the peak bolometric luminosity of simulation 12 (dark green line in Figure 11). This is consistent with the subplot of Figure 8, where simulation 1 (dark green bar) has an ${X}_{\mathrm{Ln}+\mathrm{An}}$ that is less than 10%, while simulation 12 (dark green bar) has an ${X}_{\mathrm{Ln}+\mathrm{An}}$ of about 30%. We caution that while lanthanides and actinides produce a high-opacity component, the relationship is one of diminishing returns. Increasing the lanthanide mass fraction from 0.1–1.0 raises the mean opacity by a factor of less than 2 and even an increase from 0.01–1.0 results in a change of a factor of ≲5 (Kasen et al. 2013). These differences tend to be less significant than that of the variation in nuclear heating across models. If we were working in the context of a fixed heating rate which is insensitive to nuclear physics variations, then the ejecta mass would be tightly constrained by a comparison with observed luminosity, and with this information about the mass, then the opacity could be constrained by matching the evolution of the light curve. Since in our study the heating rate is calculated to be consistent with the nuclear model, we can see from Figure 11 that the (smaller) differences in opacity can be offset by a combination of the heating rate and ejecta mass, or vice versa.

The combination of the effective heating and ${X}_{\mathrm{Ln}+\mathrm{An}}$ leads to substantial differences in the duration of the bolometric light-curve peaks shown in Figure 11. Simulations 3 and 4 (blue lines in Figure 11) have a fast rise at about 1 day, and their peak luminosity begins to decline only around 12 days. This correlates with the shallow slope in the effective heating rate which comes from an early time fission contribution, as well as a substantial ${X}_{\mathrm{Ln}+\mathrm{An}}$, which increases the diffusion time of photons through the ejecta. These long plateaus are interesting as a potential signal of early time fission and need to be verified with a more sophisticated radiation transport model (presented in Barnes et al. 2020).

Figure 11 demonstrates that the single-component (red) model of a given mass and velocity, as outlined in Section 2, could have a wide range of predicted luminosities, depending on the nuclear physics and astrophysics inputs chosen for the r-process calculation. We now turn to an exploration of the uncertainty in the inferred ejecta mass given uncertainties in these inputs in the case where there is a component dominated by lanthanide/actinide rich material.

In Figure 12(a), for simulations 12 and 13 (which are the linear combinations that fit the solar abundance pattern), we show the calculated light curve for various ejecta masses ranging from 0.01–0.08 M_⊙. In all models, the average ejecta velocity was v_ej = 0.15 c. As can be seen from Figure 12(a), this grid of masses covers a substantial range of bolometric luminosities. Additionally, the peak luminosity, width of the peak, and decay of the light curve also depend on the ejecta mass. In the context of the methods outlined in Section 2, if we assume the merger event GW170817 to have synthesized nuclear species in ratios consistent with the solar abundance pattern (for A > 125) in conditions like simulations 12 or 13, 0.02 M_⊙ is the best inferred ejecta mass to match the light curve from 8–12 days. Additional work with a more complete set of solar abundance-tuned models is needed to determine the full extent of the uncertainty on the bolometric luminosity from models that match the solar abundance pattern but have different nuclear inputs.

Zoom In Zoom Out Reset image size

Next, we take simulations 1–11 and calculate a light curve for four different values of the ejecta mass: 0.01 M_⊙, 0.02 M_⊙, 0.05 M_⊙, and 0.08 M_⊙. In Figure 12(b), for each simulation, we plot the light curve that best matches the data at 4–12 days and show the corresponding ejecta mass in the key. We see from this figure that a single set of data can be fairly well fit by a variety of ejecta masses. The inferred ejecta mass can be relatively small with the simulations with higher effective heating rates and ${X}_{\mathrm{Ln}+\mathrm{An}}$. For example, the best match for simulations 3, 4, and 5 (HFB and ETFSI mass models) occurs with an ejecta mass of 0.01 M_⊙. The simulations with lower effective heating rates at the time their light-curve peaks and lower mass fractions of lanthanides and actinides, such as simulations 9 and 10 with the SLY4 mass model, have higher inferred ejecta masses, 0.08 M_⊙. The ejecta mass inferred from simulations 1, 2, 6, 7, 8, edit112, and 13, which have light curves that reflect an interplay between the amount of lanthanides and actinides and the effective heating rates, is in between these two extremes.

This analysis indicates the extent to which the inferred ejecta mass can vary with different nuclear physics and different initial neutron richness, while keeping the rest of the parameters in the KN model constant. In particular, it is clear that if the model is not required to be consistent with solar abundance ratios, the inferred ejecta mass may suffer from large uncertainties. Note that this analysis is intended to show the range of uncertainty that would exist if one considered material that is synthesized with Y_e  < 0.24 in the 8–12 day range. Our analysis considers scenarios where the ejecta reaches rather high ${X}_{\mathrm{Ln}+\mathrm{An}}$. Previous estimates of an ${X}_{\mathrm{Ln}+\mathrm{An}}\sim {10}^{-2}$ (e.g., Chornock et al. 2017 calculated based on observations are derived under specific assumptions about the r-process heating rate, which we vary. If a smaller mass fraction is preferred than is obtained from neutron-rich element synthesis, then these curves could be combined with lanthanide-free material to reach the desired lanthanide mass fraction, and the previous analysis could be repeated along with complete fits to the spectra and time evolution of the luminosity.

For some practical applications, it may be useful to approximate the light curve as a power law. For the power law t^−α we find this range to be 1 ≤ α ≤ 2.7, depending on the nuclear physics inputs, for the time period we consider, t < 40 days. The broad range of power-law indices we find reflects the fact that our analysis includes contributions from all decay channels, including α-decay and fission, which can contribute significantly to the heating. For example, simulation 1 (purple line in Figure 12(b)), is not neutron rich enough to significantly produce heavy nuclei, which undergo spontaneous fission or α-decay. Therefore, the late-time decay is not supplemented by these extra processes, and the bolometric luminosity decays rapidly, proportional to t^−2.7. Meanwhile, simulations 3–8 have more significant fissioning and/or α-decaying nuclei, and show a much slower decay rate that is closer to t⁻¹.

Other estimates of this power law at a different times can be found in KB19 , Waxman et al. (2019), and Hotokezaka & Nakar (2020).

Given the potential importance of spontaneous fission, as well as the large uncertainties associated with spontaneous fission heating, we examine the consequences of spontaneous fission rate prescriptions as well as the influence of daughter product distributions. In Figure 14, we show fractional effective heating rates in order to compare simulations for which we altered only these aspects of spontaneous fission, while holding all other nuclear and astrophysical inputs constant. The most striking takeaway from Figure 14(a) is that prior to tens of days, there can be a clear difference between the results with KZ and XR for cases that permit nuclei with high mass numbers to be significantly populated, as is the case for HFB with KZ rates. This is because XR spontaneous fission rates are relatively high starting around Z ∼ 94 regardless of the predicted fission barriers since this prescription has no dependence on this input. Therefore, in the XR case, high mass, high fission Q-value nuclei are never significantly populated. We see also from both panels of this figure that changing the daughter product distribution can influence the fraction of β-decay heating.

Zoom In Zoom Out Reset image size

Returning to all 224 simulations, we now identify those isotopes responsible for most of the spontaneous fission heating. In general, we determine key contributors at each selected time by ranking the heaters within each channel. We then sum this list until 80% of the total heating in that channel is accounted for, as we find that 80% adequately represents the diversity of significant contributors. In Figure 15, we show a section of the chart of the nuclides, graphically representing the number of simulations for which each isotope was responsible for all or part of the top 80% of the spontaneous fission heating. We find that in most cases, even at early times, ²⁵⁴Cf (t_1/2 = 60.5 ± 0.2 days) (Phillips et al. 1963) is among the top isotopes contributing to spontaneous fission heating. This is reflected in the behavior of the spontaneous fission heating lines in Figure 13, showing a consistent underlying plateau and decay after about 50 days, indicating a population and evolution of ²⁵⁴Cf on its spontaneous fission timescale. However, in some cases, we find other isotopes to be responsible for all or part of 80% of the spontaneous fission heating in that model, either because they dominate over ²⁵⁴Cf or because ²⁵⁴Cf itself does not make up the entire 80%. We list the simulations in which an isotope other than ²⁵⁴Cf is responsible for all or part of 80% of the total spontaneous fission heating in Table 3. The dominant trend is that models that permit the synthesis of heavy species near or beyond the N = 184 predicted shell closure have the greatest diversity in nuclear species that contribute to spontaneous fission heating, as is the case with HFB + KZ combinations. We note, additionally, that very low fission barriers will result in limited population of nuclei near ²⁵⁴Cf (Vassh et al. 2019; Giuliani et al. 2020), and thus little to no contribution of spontaneous fission to the nuclear heating. At times much earlier than the timescale for spontaneous fission of ²⁵⁴Cf, e.g., at 1 day, our simulations have a more diverse set of important spontaneous fission heaters than at later times, especially from nuclei that decay on shorter timescales.

Zoom In Zoom Out Reset image size

Table 3. Top Contributing Nuclei: Spontaneous Fission (80%)

Time (days)	Isotope	Model	Fission Yield	Fission Rate	Ye
1	${}_{96}^{254}\mathrm{Cm}$	HFB22	Symmetric	XR	0.24*
...			Kodama	XR	0.24*
...	${}_{100}^{258}\mathrm{Fm}$	HFB27	Symmetric	KZ	(0.02–0.24)*
...			Kodama	KZ	(0.02–0.24)*
...	${}_{104}^{267}\mathrm{Rf}$	FRDM2012	Kodama	KZ	0.12, 0.14
...		UNEDF1	Symmetric	KZ	0.14, 0.18
...			Kodama	KZ	(0.02–0.18)
...		HFB22	Symmetric	KZ	0.02*, 0.12, 0.14, 0.16*, 0.18*, 0.21, 0.24
...			Kodama	KZ	(0.02–0.18)*, 0.21, 0.24
...	${}_{104}^{270}\mathrm{Rf}$	HFB22	Symmetric	KZ	0.02*
...	${}_{104}^{271}\mathrm{Rf}$	HFB22	Symmetric	KZ	0.02*, 0.12, 0.14, 0.16*
...			Kodama	KZ	(0.02–0.16)*
...		HFB27	Symmetric	KZ	(0.02–0.18)*
...			Kodama	KZ	(0.02–0.18)*
...	${}_{105}^{273}\mathrm{Db}$	UNEDF1	Kodama	KZ	0.12
...	${}_{108}^{288}\mathrm{Hs}$	ETFSI	Symmetric	KZ	0.02, 0.12, 0.14*
...			Kodama	KZ	0.14
8	${}_{100}^{259}\mathrm{Fm}$	HFB22	Symmetric	KZ	0.18
...			Kodama	KZ	0.18
...	${}_{104}^{269}\mathrm{Rf}$	HFB22	Symmetric	KZ	0.16, 0.18
...			Kodama	KZ	0.02, (0.14–0.18)
...		HFB27	Symmetric	KZ	0.02, 0.12, 0.16, 0.18
...			Kodama	KZ	0.02, 0.12, 0.16, 0.18
...	${}_{104}^{270}\mathrm{Rf}$	HFB27	Symmetric	KZ	0.02, 0.18
...			Kodama	KZ	0.02, 0.14
50	${}_{100}^{259}\mathrm{Fm}$	HFB22	Symmetric	KZ	0.18
...			Kodama	KZ	0.18
...	${}_{104}^{269}\mathrm{Rf}$	HFB22	Symmetric	KZ	0.16
...			Kodama	KZ	0.02, 0.16, 0.18
...		HFB27	Symmetric	KZ	0.02, 0.12, 0.16, 0.18
...			Kodama	KZ	0.02, 0.12, 0.14, 0.16, 0.18

Note. Simulations in which an isotope other than ²⁵⁴Cf was responsible for all or part of 80% of the total spontaneous fission heating. The right columns list out the specific simulations in which each isotope was seen, labeled first by nuclear model, then in parentheses by spontaneous fission yield, spontaneous fission rate, and Y_e . An asterisk indicates that for the corresponding simulations, ²⁵⁴Cf did not appear in the list of top 80%.

Download table as:  ASCIITypeset image

As Figure 14 shows, varying the nuclear inputs, such as the mass model or the spontaneous fission rate and yield distribution, can also impact the α-decay heating prediction. Depending on the specific model inputs, we find that α-decay heating can range from 10% of the total heating to up to 50% prior to 100 days. Following the same procedure described in Section 5.1, we identify isotopes responsible for the top 80% of α-decay heating. In Figure 16, we show the evolution of the top contributing nuclei between 1, 8, and 30 days, with a more complete list included in the Appendix. At each of the selected times, there is generally good agreement across simulations regarding which isotopes contribute the most.

Zoom In Zoom Out Reset image size

We distinguish two regions of importance: 83 ≤ Z ≤ 89 and 98 ≤ Z ≤ 105. The first region, shown in the bottom left corner of each subplot of Figure 16, is accessible to most trajectories, and is a region where the preferential decay mode tends to be α-decay. The magnitude of the α-heating in this region is largely determined by the amount of material that reaches these high-mass regions, thus we find that it exhibits a relatively high sensitivity to Y_e . Comparing the important isotopes in this region as a function of time (across panels), we see the expected evolution toward stability.

Meanwhile, in the upper right region of the nuclear chart (higher atomic numbers), decay modes display a stronger competition between α-decay and spontaneous fission, e.g., one species might primarily decay via α-decay with a secondary or tertiary spontaneous fission decay, while its immediate neighbor might primarily decay via spontaneous fission. Therefore, while the α-heating in this region is sensitive to the amount of material that reaches the region, it is also sensitive to the relative rates of α-decay and spontaneous fission along the relevant decay chains. By examining the ordered list of contributors to α-decay heating, we find that in simulations where isotopes from this second region are important, they tend to supply an increasing fraction of the total α-decay heating as time progresses. We find that the two most consistent isotopes from the second region, ²⁵³Es and ²⁵⁵Fm, are present in most nuclear model simulations from this sample, but are restricted to simulations using lower Y_e values (only below 0.18). Finally, while many nuclei in this region have measured α-decay rates, some do not, making both theoretical predictions and measurements for these α-decays particularly important.

We find the largest variation in fractional effective heating in the contribution from β-decay heating. Depending on the time post-merger, varying the fission prescription, nuclear models, and astrophysical conditions can result in β-decay heating comprising anywhere between 10% to more than 90% of the total effective heating. The spontaneous fission reactions and α-decays discussed in the previous sections have direct impacts on the β-decay heating, as these processes are a significant factor in determining the population of nuclei that eventually β-decay. Figure 17 shows the frequency of isotopes responsible for all or part of 80% of the β-decay heating at 8 (left) and 50 (right) days: a more complete list is found in the Appendix. Close to stability, there are a few important isotopes that show up in nearly every simulation. Further away from stability, and as early as 30 days, we start to see the significant impact of ²⁵⁴Cf spontaneous fission daughter products in the A = 130 region. Given that the simulations we compare in this section span Y_e values from 0.02–0.24, it is expected that the abundance of ²⁵⁴Cf varies from simulation to simulation, and therefore the fraction of β-heating that comes from decays of daughter products of ²⁵⁴Cf is also not uniform; this is reflected in Figure 17. Surrounding the two distinct green patches seen in the bottom right panel of Figure 17 is a more diffuse distribution of isotopes (yellow patches) that appear to be significant in few simulations. These are primarily found in simulations using HFB nuclear models with Kodama fission yields and KZ spontaneous fission rates. This highlights the potential importance of isotopes with unmeasured characteristics in the high-Z region (100 ≤ Z ≤ 106), both directly via spontaneous fission or α-decay, as well as indirectly via the decay products of these processes. The amount of heating from spontaneous fission, α-decay, and subsequently, β-decay, is ultimately sensitive to whether the synthesis of heavy nuclei reaches this high-Z region, and if it does, the process by which they preferentially decay.

Zoom In Zoom Out Reset image size