Can Kilonova Light curves be Standardized?

TL;DR

This study explores whether kilonova light curves from neutron star mergers can be standardized for distance measurement, using theoretical models to identify potential correlations similar to those in supernovae.

Contribution

The paper demonstrates, through theoretical modeling, that kilonova light curves exhibit significant correlations, suggesting they could be standardized in future observations.

Findings

Theoretical models show strong correlations in kilonova light curves.

Potential for standardizing kilonovae for cosmological distance measurements.

Future observations could confirm these correlations in real data.

Abstract

Binary neutron star mergers have been recently confirmed to be the progenitors of the optical transients kilonovae (KNe). KNe are powered by the radioactive decay of neutron-rich elements (r-process elements) which are believed to be the product of disruption of neutron stars during their merger. KNe exhibit interesting parallels with type Ia supernovae (SNe), whose light curves show specific correlations which allow them to be used as standardizable candles. In this paper, we investigate the possibility of the KN light curves exhibiting similar correlations. While a satisfactory answer to this question can only be provided by future KN observations, employing theoretical models we explore whether there is any ground for harboring such expectations. Using semi-analytic models of KN light curves in conjunction with results from numerical relativity simulations of binary neutron star mergers, we obtain the maximum bolometric luminosity ($L_{\mathrm{Bol}}^{\mathrm{max}}$) and decline in luminosity ($\Delta L_{\mathrm{Bol}}$) for a simulated population of mergers. We find that theoretical light curves of KNe show remarkable correlations despite the complex physics governing their behavior. This presents a possibility of future observations to uncover such correlations in the observed light curves, eventually allowing observers to standardize these light curves and to use them for local distance measurements.

Equations

$$L_{\mathrm{in},m}(t)$$

$$ϵ(t)$$

$$L_m(t)=\exp \left(\frac{-t^2}{{\tau }_m^2}\right){\int }_0^{t}2L_{\mathrm{in},m}(t^{\prime })ϵ(t^{\prime })\exp \left(\frac{t^{\prime 2}}{{\tau }_m^2}\right)\frac{t^{\prime }}{{\tau }_m^2}d{t}^{\prime }$$

Full text

\address

1 International Centre for Theoretical Sciences, Tata Institute of Fundamental Research, Bengaluru 560089, India \address2 Canadian Institute for Advanced Research, CIFAR Azrieli Global Scholar, MaRS Centre, West Tower, 661 University Ave, Toronto, ON M5G 1M1, Canada

1. Introduction

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 (sGRB) 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 (hereafter, 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 $\sim$1040 - 1042 erg s*-1*.

A detailed description of BNS mergers has been obtained using large scale numerical relativity (NR) simulations by various groups (see, e.g., Shibata & Hotokezaka 2019; Tanaka 2016; Rosswog 2015; Faber & Rasio 2012 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) (Radice et al., 2018; Coughlin et al., 2018). Radiative transfer calculations performed using data from such numerical simulations obtain light curves that qualitatively agree with the observations of GW170817 (Tanaka et al., 2017; Kasen 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) 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^{-1}$ followed by $t^{-3}$ (Arnett, 1982; Li & Paczyński, 1998; Metzger et al., 2010; Chatzopoulos et al., 2012; Korobkin et al., 2012; Metzger, 2017; Barnes et al., 2016; Villar et al., 2017; Arcavi et al., 2017).

It is interesting to note how much of this picture of the electromagnetic transient described above is similar to that of type Ia supernovae (SNe Ia). SNe Ia are bright optical transients (typical luminosities $\sim$1043 - 1044erg s*-1*) formed from thermonuclear explosion of white dwarfs (Hoyle & Fowler, 1960; Wheeler & Harkness, 1990) whose light curves are powered by radioactive decay of Ni56. 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}}^{\mathrm{max}}$) and the decline in luminosity ($\Delta\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 paper is organized as follows: Sec. 2 provides a summary of synthetic light curve models that we are employing along with the NR fitting formulas for ejecta properties. In Sec. 3 we discuss our main results while Sec.  4 presents a summary and outlook.

2. Semi-analytical modeling of Kilonova light curves

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 paper 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 $\gamma$-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; Radice et al., 2018; Coughlin et al., 2018).

Ejecta masses and velocities are functions of gravitational and baryonic masses, mass-ratio and weighted average $\tilde{\Lambda}$ of individual tidal deformabilities (Flanagan & Hinderer, 2008). From the NS masses in the binary, assuming the DD2 EOS (Typel et al., 2010; Hempel & Schaffner-Bielich, 2010), we obtain the masses and velocities of the dynamical ejecta and disk mass using NR fits provided by Radice et al. (2018) (Eqn 18-25) and Coughlin et al. (2018) (Eqn 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 $\sim$ 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.

In our model, the total ejecta mass is then decomposed into “blue”, “purple” and “red” components, which differ in their electron fraction $Y_{\mathrm{e}}$ and hence the opacity $\kappa$. Following Villar et al. (2017), we assume $\kappa_{m}=0.5,3,10~{}\mathrm{cm}^{2}~{}\mathrm{g}^{-1}$ 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}}$ to denote array of the fraction of ejecta mass distributed in to the blue, purple and red components (${f_{\rm e}}_{m}$ being $m$th element of this array): for example, ${f_{\rm e}}=[0.2,0.6,0.2]$ means that total ejecta mass is decomposed into 20% blue ($\kappa=0.5~{}\mathrm{cm}^{2}~{}\mathrm{g}^{-1},~{}Y_{\mathrm{e}}>0.4$), 60% purple ($\kappa=3~{}\mathrm{cm}^{2}~{}\mathrm{g}^{-1},~{}0.1<Y_{\mathrm{e}}<0.4$) and 20% red ($\kappa=10~{}\mathrm{cm}^{2}~{}\mathrm{g}^{-1},Y_{\mathrm{e}}<0.1$) components. Note that ${f_{\rm e}}$ 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}}$ (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 modelling 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.

[TABLE]

where $M_{\mathrm{rp},\,m}$ is the total mass (in $\mathrm{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 2-d 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$

[TABLE]

where, $\tau_{m}=\sqrt{{2\kappa_{m}M_{\mathrm{rp},\,m}}/{\beta v_{\mathrm{ej}}c}}$, with $\kappa_{m}$ being is the gray opacity of the ejecta component, and $\beta=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_{\mathrm{rp},~{}m}=[0.02,0.047,0.011]\mathrm{M}_{\odot}$ 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}=1.36-1.6\mathrm{M}_{\odot}$; $M_{2}=1.16-1.36\mathrm{M}_{\odot}$, with the constraint $M_{1}+M_{2}=2.73^{+0.04}_{-0.01}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}}=[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.

3. Results

We generate a population of BNS mergers whose masses are uniformly distributed in the mass range $1.2-1.7~{}\mathrm{M}_{\odot}$ 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}}^{\mathrm{max}}$ and decrease in luminosity $\Delta\log L_{\mathrm{bol}}$ in 5 days from the time of the peak luminosity, where $\Delta\log L_{\mathrm{bol}}\equiv\log(L_{\mathrm{bol}}^{\mathrm{max}}/L_{\mathrm{bol}}^{\mathrm{5~{}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 Radice et al. 2018; Coughlin 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_{\mathrm{e}}$ of the ejecta (${f_{\rm e}}=[0.2,0.6,0.2]$ and ${f_{\rm e}}=[0.26,0.6,0.14]$ used by  Villar et al. 2017). We also consider a case where ${f_{\rm e}}$ depends on the NS masses111We use figure 5 from Dietrich & Ujevic (2017) to get the average electron fraction of the ejecta, $\bar{Y}_{e}(q)$, as a function of the mass ratio $q:=M_{1}/M_{2}$ of the binary. Then we solve two constraint equations, $\sum_{m}{f_{\rm e}}_{m}=1$ and $\sum{f_{\rm e}}_{m}\,{Y_{\mathrm{e}}}_{m}=\bar{Y}_{e}(q)$, where $m$ denotes the blue, purple and red components. The system is under-determined as there are three unknowns (${f_{\rm e}}_{m}$) and only two equations. We choose ${f_{\rm e}}_{\mathrm{blue}}=0.1$ for the blue component to obtain the same for the red and purple components. Corresponding results are plotted in figure 3. A choice ${f_{\rm e}}_{\mathrm{blue}}=0.2$ does not make a big difference.. 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}}^{\mathrm{max}}$ and $\Delta\log L_{\mathrm{bol}}$.

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 which 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.

4. Summary and outlook

Motivated by the similarities of KNe with SNe Ia, we have explored the possibility of KNe providing a set of standardizable candles 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}}^{\mathrm{max}}$ and the decline in the luminosity $\Delta\log L_{\mathrm{bol}}$ after a few days (figure 3). This is performed by computing $L_{\mathrm{bol}}^{\mathrm{max}}$ and $\Delta\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 that there are however key differences between SN and KNe light curves. KNe have peak luminosities ($\sim$1040 - 1042 erg s*-1*) much lower than SNe Ia peak luminosities ($\sim$1043 - 1044 erg s*-1*) making KNe observable only in the local universe. The SNe Ia B-band light curve usually peaks $\sim$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.

Acknowledgments:

We thank Kenta Hotokezaka, Bala Iyer, Ramesh Narayan, Sumit Kumar, 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. RK 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.