Can Kilonova Light Curves Be Standardized?

The gravitational wave (GW) event GW170817 marked the birth of a new era in multi-messenger astrophysics (Abbott et al. 2017b). The event was associated to a binary neutron star (BNS) merger located nearly 40 Mpc away in the galaxy NGC 4993. The GW trigger was followed by a nearly coincident short gamma-ray burst 1.7 s after the merger time (Abbott et al. 2017a). The gamma-ray counterpart was subsequently followed up by several ground- and space-based telescopes in the ultraviolet, optical, and near-infrared (UVOIR) bands (Abbott et al. 2017c). The UVOIR emission is largely identified to be from a quasi-thermal transient called a kilonova (KN), powered by radioactive decay of several r-process nuclei (Li & Paczyński 1998; Metzger et al. 2010; Barnes & Kasen 2013), resulting in typical luminosities of ∼10⁴⁰–10⁴² erg s⁻¹.

A detailed description of BNS mergers has been obtained using large-scale numerical relativity (NR) simulations by various groups (see, e.g., Faber & Rasio 2012; Rosswog 2015; Tanaka 2016; Shibata & Hotokezaka 2019 and references therein). These groups have found that the properties (e.g., masses and velocities) of r-process radioactive material ejected from the merger tends to exhibit a non-trivial dependency on the neutron star (NS) masses and tidal deformabilities. Various groups have provided formulae that gives excellent fits to the ejecta properties as a function of the NS masses and equation of state (EOS; Coughlin et al. 2018; Radice et al. 2018). Radiative transfer calculations performed using data from such numerical simulations obtain light curves that qualitatively agree with the observations of GW170817 (Kasen et al. 2017; Tanaka et al. 2017; Tanvir et al. 2017; Miller et al. 2019). Additionally, semi-analytical models based on Arnett–Chatzopoulos framework (Arnett 1982; Chatzopoulos et al. 2012) (originally developed for Type Ia supernovae (SNe Ia)) have been able to explain the light curve fairly well. In such models, the heating rate is calculated using radioactive decay of mostly r-process elements, giving the initial intensity decline rate to be t⁻¹ followed by t⁻³ (Arnett 1982; Li & Paczyński 1998; Metzger et al. 2010; Chatzopoulos et al. 2012; Korobkin et al. 2012; Barnes et al. 2016; Arcavi et al. 2017; Metzger 2017; Villar et al. 2017).

It is interesting to note how much of this picture of the electromagnetic transient described above is similar to that of SNe Ia. SNe Ia are bright optical transients (typical luminosities ∼10⁴³–10⁴⁴ erg s⁻¹) formed from the thermonuclear explosion of white dwarfs (Hoyle & Fowler 1960; Wheeler & Harkness 1990) whose light curves are powered by radioactive decay of Ni⁵⁶. The double degenerate model proposes binary white dwarf mergers as a progenitors of SNe Ia (Iben & Tutukov 1984; Webbink 1984; Nelemans et al. 2001). This model is being supported by several recent observational and theoretical studies, particularly in terms of it being able to explain the Galactic birth rates and delay time distributions (Ruiter et al. 2009; Maoz et al. 2014; Kashyap et al. 2015).

It is well known that SNe Ia light curves show an empirical relationship between the maximum intrinsic luminosity and the decline rate, known as the "Phillips relation" or the "width-luminosity relation" (Phillips 1993). By estimating the maximum intrinsic luminosity from the observed decline rate using the Phillips relation they can be used as a standard candle (Riess et al. 1998). There are several parallels between SNe Ia and KNe: both are believed to be triggered by the merger of compact objects in narrow mass ranges and are powered by radioactive decay of heavy isotopes. Moreover empirical models based on Arnett–Chatzopoulos framework seem to agree with the observations. Hence, it is quite natural to ask this question: could KNe light curves be standardized like SNe Ia light curves?

A definitive answer to this question can only be provided by a large number of KNe observations, as happened in the case of SNe Ia. As we wait for such observations (Yang et al. 2017), we explore whether there is any ground for harboring such expectations. This is done by investigating whether any correlations exist in the synthetic KNe light curves provided by semi-analytical models in conjunction with results from NR simulations of BNSs. Making use of NR fitting formulae for the ejecta properties, we generate synthetic light curves from several putative KNe produced by the merger of several simulated BNS systems with different component masses. We then investigate the correlation between the peak luminosity (${L}_{\mathrm{Bol}}^{\max }$) and the decline in luminosity (${\rm{\Delta }}\mathrm{log}{L}_{\mathrm{Bol}}$) after 5 days following the peak. We find that "Phillips-like" relations exist in these synthetic light curves.

Indeed, the current semi-analytical models are unlikely to capture the complex physics and the rich phenomenology of KNe in entirety. Hence, the specific relationship that we find using the current semi-analytic models are unlikely to hold up against actual observations. However, they hint a possibility of the existence of such relationships in real light curves. This Letter is organized as follows. Section 2 provides a summary of synthetic light curve models that we are employing along with the NR fitting formulas for ejecta properties. In Section 3 we discuss our main results, while Section 4 presents a summary and outlook.

In the absence of a large enough number of KNe observations, here we explore the possibility of the existence of a Phillips-like relation in the synthetic light curves predicted by semi-analytical KNe models. In this Letter we adopt the Arnett–Chatzopoulos model to obtain the light curves for a population of binary NS mergers using the NR fitting formulae for ejecta mass and velocity provided by Radice et al. (2018) and Coughlin et al. (2018).

The light curve modeling justifiably assumes that the physical processes responsible for heating that produces UVOIR are well separated in time from the processes responsible for γ-rays, X-rays, and radio (Hotokezaka et al. 2016). In addition it is assumed that there is a homologously expanding isotropic ejecta of neutron-rich radioactive isotopes. This expanding ejecta will follow the same evolution in an ambient medium as seen for SNe, for example. The model constructs the light curve by taking into account the work done in expansion, the heating done by radioactive decay along with the knowledge of velocity and opacity of the expanding ejecta (Arnett 1982; Chatzopoulos et al. 2012). The properties of the ejecta, in turn, depend on the mass of the NSs in the binary, for a given EOS as shown by NR simulations (Dietrich & Ujevic 2017; Coughlin et al. 2018; Radice et al. 2018).

Ejecta masses and velocities are functions of gravitational and baryonic masses, mass ratio, and weighted average $\tilde{{\rm{\Lambda }}}$ of individual tidal deformabilities (Flanagan & Hinderer 2008). From the NS masses in the binary, assuming the DD2 EOS (Hempel & Schaffner-Bielich 2010; Typel et al. 2010), we obtain the masses and velocities of the dynamical ejecta and disk mass using NR fits provided by Radice et al. (2018; Equations (18)–(25)) and Coughlin et al. (2018; Equations (D1)–(D5)). We assume that 30% of the disk mass becomes unbound and contributes as the disk ejecta. The total ejecta mass, ${M}_{\mathrm{ej}}$ is then taken to be the sum of the dynamical and disk ejecta. Figure 1 shows the ejecta mass and velocity as a function of the NS masses in the binary. We find that the numerical fits of the eject mass (velocity) provided by both groups differ, on average, by ∼29% (12%), which is a reflection of the error in these estimates. The NS radius decreases with the mass; hence, the lower-mass companion is more prone to tidal deformation, producing larger ejecta mass. Also, since more massive NSs are more compact, they would also produce larger ejecta velocities. We observe these trends in Figure 1.

Zoom In Zoom Out Reset image size

In our model, the total ejecta mass is then decomposed into "blue," "purple," and "red" components, which differ in their electron fraction ${Y}_{{\rm{e}}}$ and hence the opacity κ. Following Villar et al. (2017), we assume κ_m = 0.5, 3, 10 cm² g⁻¹ for the blue, purple, and red components, respectively. Finally, we add the light curve due to the three components to obtain bolometric luminosity. We use the symbol ${{f}_{{\rm{e}}}}_{m}$ to denote array of the fraction of ejecta mass distributed in to the blue, purple, and red components: for example, ${{f}_{{\rm{e}}}}_{m}=[0.2,0.6,0.2]$ means that total ejecta mass is decomposed into 20% blue ($\kappa =0.5\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1},{Y}_{{\rm{e}}}\gt 0.4$), 60% purple ($\kappa =3\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1},0.1\lt {Y}_{{\rm{e}}}\lt 0.4$), and 20% red ($\kappa \,=10\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1},{Y}_{{\rm{e}}}\lt 0.1$) components. Note that ${{f}_{{\rm{e}}}}_{m}$ of dynamical and disk ejecta can be, in general, different (see, e.g., Radice et al. 2018 and Lippuner et al. 2017). Leaving a more careful treatment of this for a future work, we consider three simple choices of ${{f}_{{\rm{e}}}}_{m}$ (see Section 3). In the absence of accurate predictions from NR simulations, we assume the same velocity ${v}_{\mathrm{ej}}$ for all components of the ejecta.

The analytical modeling of light curve depends on the input heating rate and thermal efficiency of each component of the ejecta, given by Korobkin et al. (2012) and Barnes et al. (2016), respectively.

$$\begin{eqnarray}&&\begin{array}{l}{L}_{\mathrm{in},m}(t)\\ =\,4\times {10}^{18}\ {M}_{\mathrm{rp},m}\ {\left[0.5-{\pi }^{-1}\arctan \left.\left(\displaystyle \frac{t-{t}_{o}}{\sigma }\right)\right]\right.}^{1.3}\ \mathrm{erg}\ {{\rm{s}}}^{-1}\\ \epsilon (t)=0.36\left[\exp (-a\,t)+\displaystyle \frac{\mathrm{ln}(1+2\,{{bt}}^{d})}{2\,{{bt}}^{d}}\right]\end{array}\end{eqnarray} \tag{ 1 }$$

where ${M}_{\mathrm{rp},m}$ is the total mass (in ${\text{}}{M}_{\odot }$) of the r-process elements synthesized for each component m (blue, purple, or red), i.e., ${M}_{\mathrm{rp},m}={{f}_{{\rm{e}}}}_{m}{M}_{\mathrm{ej}}$, and t is time in days. We use 2D interpolation for each of the ejecta components to obtain the values for the fit parameters a, b, d for different ejecta masses and velocities using the Table 1 of Barnes et al. (2016).

Assuming homologous expansion as described in Arnett (1982), we use the prescription outlined in Chatzopoulos et al. (2012) and Villar et al. (2017) to compute the luminosity for each component m

$$\begin{eqnarray}&&{L}_{m}(t)=\exp \left(\displaystyle \frac{-{t}^{2}}{{\tau }_{m}^{2}}\right){\int }_{0}^{t}2\,{L}_{\mathrm{in},m}(t^{\prime} )\,\epsilon (t^{\prime} )\,\exp \left(\displaystyle \frac{t{{\prime} }^{2}}{{\tau }_{m}^{2}}\right)\,\displaystyle \frac{t^{\prime} }{{\tau }_{m}^{2}}\,{dt}^{\prime} \end{eqnarray} \tag{ 2 }$$

where, ${\tau }_{m}=\sqrt{2{\kappa }_{m}{M}_{\mathrm{rp},m}/\beta {v}_{\mathrm{ej}}c}$, with κ_m being is the gray opacity of the ejecta component, and β = 13.4, a dimensionless constant associated with the geometric profile of the ejecta (Villar et al. 2017).

Figure 2 shows the synthetic light curves computed for the KN associated with GW170817, along with the observed ones. Villar et al. (2017)'s best-fit model uses M_rp,m = [0.02, 0.047, 0.011] M_⊙ and ${v}_{\mathrm{ej},m}=[0.266,0.152,0.137]\,c$ for the blue, purple, and red components of the ejecta. We also plot the light curves computed using the estimated component masses from GW170817 (M₁ = 1.36–1.6 M_⊙; M₂ = 1.16–1.36 M_⊙, with the constraint ${M}_{1}+{M}_{2}\,={2.73}_{-0.01}^{+0.04}\,{M}_{\odot }$ Abbott et al. 2017b), where the ejecta mass and velocity estimated using NR fitting formulae of Radice et al. (2018) and Coughlin et al. (2018). Here, as discussed earlier, we assume that ${{f}_{{\rm{e}}}}_{m}=[0.2,0.6,0.2]$. For comparison, we also plot the observed light curves as presented by Drout et al. (2017) and Cowperthwaite et al. (2017). The general agreement between the theoretical models and observations is encouraging.

Zoom In Zoom Out Reset image size

We generate a population of BNS mergers whose masses are uniformly distributed in the mass range 1.2–1.7 M_⊙ and compute the synthetic light curves produced by each merger, using the procedure outlined in Section 2. We use these theoretical light curves to find the relation between maximum luminosity ${L}_{\mathrm{bol}}^{\max }$ and decrease in luminosity ${\rm{\Delta }}\mathrm{log}{L}_{\mathrm{bol}}$ in 5 days from the time of the peak luminosity, where ${\rm{\Delta }}\mathrm{log}{L}_{\mathrm{bol}}\equiv \mathrm{log}({L}_{\mathrm{bol}}^{\max }/{L}_{\mathrm{bol}}^{5\,\mathrm{days}})$. The choice of 5 days is arbitrary, but is motivated by the fact that UVOIR observations of KNe can be typically performed over a few days. We vary the parameters in the model and discuss the possible variations in the correlation. In particular, we vary the choice of NR based fitting formula for the ejecta mass and velocity (provided by Coughlin et al. 2018; Radice et al. 2018), the nuclear EOS (DD2 by Typel et al. 2010 and WFF2 by Wiringa et al. 1988), and distribution of the electron fraction ${Y}_{{\rm{e}}}$ of the ejecta (${{f}_{{\rm{e}}}}_{m}=[0.2,0.6,0.2]$ and ${{f}_{{\rm{e}}}}_{m}=[0.26,0.6,0.14]$ used by Villar et al. 2017). We also consider a case where ${{f}_{{\rm{e}}}}_{m}$ depends on the NS masses.³ We also examine the same correlations computed assuming a time delay of 7 days from peak luminosity. These results are plotted in Figure 3, suggesting clear correlations between ${L}_{\mathrm{bol}}^{\max }$ and ${\rm{\Delta }}\mathrm{log}{L}_{\mathrm{bol}}$.

Zoom In Zoom Out Reset image size

It should be mentioned here that the relation found in this work factors in the full non-linearity in the NR simulations and the non-trivial relationship between NS masses and bolometric light curve. The fact that such a correlation has been observed in the synthetic light curves suggests that a similar correlation should exist in the actual light curves, even though the actual observed correlation may turn out to be different than what are presented here. If this is vindicated by future KNe observations, this will provide an independent distance ladder. For example, in Figure 3, the decline in luminosity in 5 days (or any other suitably chosen time) is an independent observable that can be used to find the maximum intrinsic luminosity using the correlation. Thus the luminosity distance can be estimated by comparing the intrinsic and apparent luminosities.

Motivated by the similarities of KNe with SNe Ia, we have explored the possibility of KNe providing a set of standardizable candles that are analogous to SNe Ia. Indeed, such a possibility can be confirmed or refuted only by a large number KNe observations. As we await such observations, we studied simple semi-analytical KNe models (in conjunction with NR fitting formulae for ejecta mass and velocity) and discovered correlations that exist between the peak bolometric luminosity ${L}_{\mathrm{bol}}^{\max }$ and the decline in the luminosity ${\rm{\Delta }}\mathrm{log}{L}_{\mathrm{bol}}$ after a few days (Figure 3). This is performed by computing ${L}_{\mathrm{bol}}^{\max }$ and ${\rm{\Delta }}\mathrm{log}{L}_{\mathrm{bol}}$ from synthetic light curves generated from ejecta produced by a population of BNS mergers. We employ different NR fitting formulae, NS EOS, and electron fraction distribution of the ejecta to study the robustness of our results. We note that Coughlin et al. (2019) has taken a different approach to the standardization of KNe.

The light curve calculation presented in this work has multiple simplifying assumptions, and are subject to errors in the NR simulations and the KNe models. Hence they are only crude estimates. In particular, anisotropy of the ejecta components and time-variation of the ejecta opacities needs future investigation. However, there is preliminary evidence (coming from the observation of KNe associated with GW170817) that they capture the essential features of KNe light curves. We stress the fact that we are not proposing any particular correlation, which has to be left to future observations. Such a correlation in the observed light curves could have potential usage in distance measurement, which will have applications in fundamental physics, astrophysics, and cosmology. Possible applications include constraining the number of spacetime dimensions by comparing distance estimates from GW and KN observation from a BNS merger (along the lines of Pardo et al. 2018; Abbott et al. 2019), the calibration of distance ladders (along the lines of Gupta et al. 2019), and potentially the estimation of cosmological parameters (along the lines of Coughlin et al. 2019).

We admit, however, that there are key differences between SN and KNe light curves. KNe have peak luminosities (∼10⁴⁰–10⁴² erg s⁻¹) that are much lower than SNe Ia peak luminosities (∼10⁴³–10⁴⁴ erg s⁻¹) making KNe observable only in the local universe. The SNe Ia B-band light curve usually peaks ∼20 days post explosion (Riess et al. 1999) where the spectral peak shifts from optical to infrared in about 2–3 months. In contrast, the KNe light curve reaches the maximum value within few hours and shifts from optical to infrared within 10 days. Because of these differences, KNe light curves are relatively more difficult to standardize and also the counterpart of the "Phillips relation" would be expected to be different for them. Only future observations can tell the full story.

We thank Kenta Hotokezaka, Bala Iyer, Ramesh Narayan, Sumit Kumar, and K. Haris for their patient feedback to the work. We also thank Bharat Kumar for providing polytropic fits to the EOS DD2, and David Radice and Tim Dietrich for clarifications on their NR fits. R.K. would like to thank G. Srinivasan for his encouragement to pursue the problem during a meeting in 2017 at ICTS, Bengaluru. We are grateful to the anonymous referee for many useful comments and suggestions on an earlier version of the manuscript. This research was supported by the Max Planck Society through a Max Planck Partner Group at ICTS-TIFR and by the Canadian Institute for Advanced Research through the CIFAR Azrieli Global Scholars program. Computations were performed at the ICTS cluster Alice.