Can Kilonova Light curves be Standardized?
DD RAFT VERSION F EBRUARY
6, 2020
Preprint typeset using L A TEX style emulateapj v. 12/16/11
CAN KILONOVA LIGHT CURVES BE STANDARDIZED? R AHUL K ASHYAP , G AYATHRI R AMAN , P ARAMESWARAN A JITH , International Centre for Theoretical Sciences, Tata Institute of Fundamental Research, Bengaluru 560089, India and Canadian Institute for Advanced Research, CIFAR Azrieli Global Scholar, MaRS Centre, West Tower, 661 University Ave, Toronto, ON M5G 1M1, Canada
Draft version February 6, 2020
ABSTRACTBinary neutron star mergers have been recently confirmed to be the progenitors of the optical transients kilo-novae (KNe). KNe are powered by the radioactive decay of neutron-rich elements (r-process elements) whichare believed to be the product of disruption of neutron stars during their merger. KNe exhibit interesting par-allels with type Ia supernovae (SNe), whose light curves show specific correlations which allow them to beused as standardizable candles. In this paper, we investigate whether KNe light curves could exhibit similarcorrelations. While a satisfactory answer to this question can only be provided by future KNe observations,employing theoretical models we explore whether there is any ground for harboring such expectations. Usingsemi-analytic models of KNe light curves in conjunction with results from numerical relativity simulations ofbinary neutron star mergers, we obtain the maximum bolometric luminosity ( L maxBol ) and decline from peak lu-minosity ( ∆ log L Bol ) for a simulated population of mergers. We find that theoretical light curves of KNe showremarkable correlations despite the complex physics governing their behavior. This presents a possibility offuture observations to uncover such correlations in the observed light curves, eventually allowing observers tostandardize these light curves and to use them for local distance measurements.
Keywords: kilonovae, supernovae: general — neutron stars, compact binary coalescence, standardizable candle INTRODUCTION
The gravitational wave (GW) event GW170817 marked thebirth of a new era in multi-messenger astrophysics (Abbottet al. 2017b). The event was associated to a binary neutronstar (BNS) merger located nearly 40 Mpc away in the galaxyNGC 4993. The GW trigger was followed by a nearly coin-cident short gamma ray burst (sGRB) 1.7 s after the mergertime (Abbott et al. 2017a). The gamma-ray counterpart wassubsequently followed up by several ground and space basedtelescopes in the ultraviolet, optical and near-infrared (here-after, UVOIR) bands (Abbott et al. 2017c). The UVOIR emis-sion is largely identified to be from a quasi-thermal transientcalled a kilonova (KN), powered by radioactive decay of sev-eral r-process nuclei (Li & Paczy´nski 1998; Metzger et al.2010; Barnes & Kasen 2013), resulting in typical luminosi-ties of ∼ - 10 erg s − .A detailed description of BNS mergers has been ob-tained using large scale numerical relativity (NR) simulationsby various groups (see, e.g., Shibata & Hotokezaka 2019;Tanaka 2016; Rosswog 2015; Faber & Rasio 2012 and ref-erences therein). These groups have found that the proper-ties (e.g., masses and velocities) of r-process radioactive ma-terial ejected from the merger tends to exhibit a non-trivialdependency on the neutron star (NS) masses and tidal de-formabilities. Various groups have provided formulae thatgives excellent fits to the ejecta properties as a function ofthe NS masses and equation of state (EOS) (Radice et al.2018; Coughlin et al. 2018). Radiative transfer calculationsperformed using data from such numerical simulations ob-tain light curves that qualitatively agree with the observa-tions 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 devel-oped for Type Ia supernovae) have been able to explain thelight curve fairly well. In such models, the heating rate iscalculated using radioactive decay of mostly r-process ele- ments, giving the initial intensity decline rate to be t − fol-lowed by t − (Arnett 1982; Li & Paczy´nski 1998; Metzgeret al. 2010; Chatzopoulos et al. 2012; Korobkin et al. 2012;Metzger 2017; Barnes et al. 2016; Villar et al. 2017; Arcaviet al. 2017).It is interesting to note how much of this picture of the elec-tromagnetic transient described above is similar to that of typeIa supernovae (SNe Ia). SNe Ia are bright optical transients(typical luminosities ∼ - 10 erg s − ) formed from ther-monuclear explosion of white dwarfs (Hoyle & Fowler 1960;Wheeler & Harkness 1990) whose light curves are poweredby radioactive decay of Ni . The double degenerate modelproposes binary white dwarf mergers as a progenitors of SNeIa (Iben & Tutukov 1984; Webbink 1984; Nelemans et al.2001). This model is being supported by several recent ob-servational and theoretical studies, particularly in terms of itbeing able to explain the Galactic birth rates and delay timedistributions (Ruiter et al. 2009; Maoz et al. 2014; Kashyapet al. 2015).It is well known that SNe Ia light curves show an empir-ical relationship between the maximum intrinsic luminosityand the decline rate, known as the “Phillips relation” or the“width-luminosity relation” (Phillips 1993). By estimatingthe maximum intrinsic luminosity from the observed declinerate using the Phillips relation they can be used as a standardcandle (Riess et al. 1998). There are several parallels betweenSNe Ia and KNe: Both are believed to be triggered by themerger of compact objects in narrow mass ranges and arepowered by radioactive decay of heavy isotopes. Moreoverempirical models based on Arnett-Chatzopoulos frameworkseem to agree with the observations. Hence, it is quite naturalto ask this question: Could KNe light curves be standardizedlike SNe Ia light curves?A definitive answer to this question can only be providedby a large number of KNe observations, as happened in thecase of SNe Ia. As we wait for such observations (Yang et al.2017), we explore whether there is any ground for harboring a r X i v : . [ a s t r o - ph . S R ] F e b . . . . . . M (M (cid:12) ) . . . . . . M ( M (cid:12) ) ejecta mass log ( M ej / M (cid:12) ) − . − . − . − . − . − . − . − . − . . . . . . . M (M (cid:12) ) . . . . . . M ( M (cid:12) ) ejecta velocity v ej / c .
15 0 .
17 0 .
19 0 .
21 0 .
23 0 .
25 1 . . . . . . M (M (cid:12) ) . . . . . . M ( M (cid:12) ) ejecta mass log ( M ej / M (cid:12) ) − . − . − . − . − . − . − . − . − . . . . . . . M (M (cid:12) ) . . . . . . M ( M (cid:12) ) ejecta velocity v ej / c .
15 0 .
17 0 .
19 0 .
21 0 .
23 0 . Figure 1.
Total ejecta mass log ( M ej / M (cid:12) ) (top panels) and ejecta velocity v ej / c (bottom panels) as a function of the NS masses M and M in the bi-nary, computed using the fits given in Radice et al. (2018) (left panels) andCoughlin et al. (2018) (right panels). We assume here that 30% of the diskmass contributes to the unbound r-process ejecta that powers the light curve. such expectations. This is done by investigating whether anycorrelations exist in the synthetic KNe light curves providedby semi-analytical models in conjunction with results fromNR simulations of BNSs. Making use of NR fitting formulaefor the ejecta properties, we generate synthetic light curvesfrom several putative KNe produced by the merger of severalsimulated BNS systems with different component masses. Wethen investigate the correlation between the peak luminosity( L maxBol ) and the decline in luminosity ( ∆ log L Bol ) after 5 daysfollowing the peak. We find that “Phillips-like” relations existin these synthetic light curves.Indeed, the current semi-analytical models are unlikely tocapture the complex physics and the rich phenomenology ofKNe in entirety. Hence, the specific relationship that we findusing the current semi-analytic models are unlikely to hold upagainst actual observations. However, they hint a possibilityof the existence of such relationships in real light curves. Thispaper is organized as follows: Sec. 2 provides a summary ofsynthetic light curve models that we are employing along withthe NR fitting formulas for ejecta properties. In Sec. 3 wediscuss our main results while Sec. 4 presents a summaryand outlook. SEMI-ANALYTICAL MODELING OF KILONOVA LIGHT CURVES
In the absence of a large enough number of KNe obser-vations, here we explore the possibility of the existence ofa Phillips-like relation in the synthetic light curves predictedby semi-analytical KNe models. In this paper we adopt theArnett-Chatzopoulos model to obtain the light curves for apopulation of binary NS mergers using the NR fitting for-mulae for ejecta mass and velocity provided by Radice et al. time [days] L bo l [ e r g / s ] Drout et al (2017)Cowperthwaite et al (2017)Vilar (2017)Abbott (2017b) + Radice (2018a)Abbott (2017b) + Dietrich (2019)
Figure 2.
The solid traces show the evolution of bolometric luminosity ofthe KNe associated with GW170817 as predicted by the semi-analytical KNemodels. The Different markers show the observed luminosity evolution fromthe same event calculated by Drout et al. (2017) and Cowperthwaite et al.(2017). (2018) and Coughlin et al. (2018).The light curve modeling justifiably assumes that the phys-ical processes responsible for heating that produces UVOIRare well separated in time from the processes responsible for γ -rays, X-rays and radio (Hotokezaka et al. 2016). In additionit is assumed that there is a homologously expanding isotropicejecta of neutron-rich radioactive isotopes. This expandingejecta will follow the same evolution in an ambient mediumas seen for SNe, for example. The model constructs the lightcurve by taking into account the work done in expansion, theheating done by radioactive decay along with the knowledgeof velocity and opacity of the expanding ejecta (Arnett 1982;Chatzopoulos et al. 2012). The properties of the ejecta, inturn, depend on the mass of the NSs in the binary, for a givenEOS as shown by NR simulations (Dietrich & Ujevic 2017;Radice et al. 2018; Coughlin et al. 2018).Ejecta masses and velocities are functions of gravitationaland baryonic masses, mass-ratio and weighted average ˜ Λ ofindividual tidal deformabilities (Flanagan & Hinderer 2008).From the NS masses in the binary, assuming the DD2 EOS(Typel et al. 2010; Hempel & Schaffner-Bielich 2010), weobtain the masses and velocities of the dynamical ejecta anddisk mass using NR fits provided by Radice et al. (2018) (Eqn18-25) and Coughlin et al. (2018) (Eqn D1-D5). We assumethat 30% of the disk mass becomes unbound and contributesas the disk ejecta. The total ejecta mass, M ej is then taken tobe the sum of the dynamical and disk ejecta. Figure 1 showsthe ejecta mass and velocity as a function of the NS massesin 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 esti-mates. The NS radius decreases with the mass; hence, thelower-mass companion is more prone to tidal deformation,producing larger ejecta mass. Also, since more massive NSsare more compact, they would also produce larger ejecta ve-locities. We observe these trends in figure 1.In our model, the total ejecta mass is then decomposed into . . . . . . . . ∆ log L bol . . . . . . l og L bo l , m a x GW170817 0 . . . . . . . . ∆ log L bol . . . . . . Radice + DD2 + f e Coughlin + DD2 + f e Radice + WFF2 + f e Coughlin + WFF2 + f e Radice + DD2 + f e Coughlin + DD2 + f e Radice + WFF2 + f e Coughlin + WFF2 + f e Radice + DD2 + f e ( q ) Figure 3.
The relation between maximum luminosity L maxbol and decline in luminosity in few days ∆ log L bol from simulated light curves. We consider masses inthe range 1 . − . M (cid:12) , two different NR fitting formulae for ejecta mass and velocity (Radice et al. (2018) and Coughlin et al. (2018); shown in different colors),two different EoSs (DD2 and WFF2; shown by markers with and without black edges), and two different choices of f e ( f e0 ≡ [0 . , . , . f e1 ≡ [0 . , . , . f e ( q ) as discussed in the text; shown by different markers). The filled and unfilled stars correspond to the GW170809 KNeobservations by Cowperthwaite et al. (2017) and Drout et al. (2017). The left panel corresponds to a delay time of 5 days and the right panel to 7 days. “blue”, “purple” and “red” components, which differ in theirelectron fraction Y e and hence the opacity κ . Following Vil-lar et al. (2017), we assume κ m = 0 . , ,
10 cm g − for theblue, purple and red components, respectively. Finally, weadd the light curve due to the three components to obtain bolo-metric luminosity. We use the symbol f e to denote array ofthe fraction of ejecta mass distributed in to the blue, purpleand red components ( f e m being m th element of this array):for example, f e = [0 . , . , .
2] means that total ejecta massis decomposed into 20% blue ( κ = 0 . g − , Y e > . κ = 3 cm g − , . < Y e < .
4) and 20% red( κ = 10 cm g − , Y e < .
1) components. Note that f e of dy-namical and disk ejecta can be, in general, different (see, e.g.,Radice et al. (2018) and Lippuner et al. (2017)). Leaving amore careful treatment of this for a future work, we considerthree simple choices of f e (see section 3). In the absence of ac-curate predictions from NR simulations, we assume the samevelocity v ej for all components of the ejecta.The analytical modelling of light curve depends on the in-put heating rate and thermal efficiency of each component ofthe ejecta, given by Korobkin et al. (2012) and Barnes et al.(2016), respectively. L in , m ( t ) = 4 × M rp , m [0 . − π − arctan (cid:18) t − t o σ (cid:19) ] . erg s − (cid:15) ( t ) = 0 . (cid:20) exp( − at ) + ln(1 + bt d )2 bt d (cid:21) (1)where M rp , m is the total mass (in M (cid:12) ) of the r-process el-ements synthesized for each component m (blue, purple orred), i.e., M rp , m = f e m M ej , and t is time in days. We use 2-dinterpolation for each of the ejecta components to obtain thevalues for the fit parameters a , b , d for different ejecta massesand 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 mL m ( t ) = exp (cid:18) − t τ m (cid:19) (cid:90) t L in , m ( t (cid:48) ) (cid:15) ( t (cid:48) ) exp (cid:18) t (cid:48) τ m (cid:19) t (cid:48) τ m dt (cid:48) (2)where, τ m = (cid:112) κ m M rp , m /β v ej c , with κ m being is the grayopacity of the ejecta component, and β = 13 .
4, a dimen-sionless constant associated with the geometric profile of theejecta (Villar et al. 2017).Figure 2 shows the synthetic light curves computed forthe KN associated with GW170817, along with the observedones. Villar et al. (2017)’s best fit model uses M rp , m =[0 . , . , . (cid:12) and v ej , m = [0 . , . , . c forthe blue, purple and red components of the ejecta. We alsoplot the light curves computed using the estimated componentmasses from GW170817 ( M = 1 . − . (cid:12) ; M = 1 . − . (cid:12) , with the constraint M + M = 2 . + . − . M (cid:12) Abbottet al. 2017b), where the ejecta mass and velocity estimatedusing NR fitting formulae of Radice et al. (2018) and Cough-lin et al. (2018). Here, as discussed earlier, we assume that f e = [0 . , . , . RESULTS
We generate a population of BNS mergers whose massesare uniformly distributed in the mass range 1 . − . (cid:12) andcompute the synthetic light curves produced by each merger,using the procedure outlined in section 2. We use these the-oretical light curves to find the relation between maximumluminosity L maxbol and decrease in luminosity ∆ log L bol in 5days from the time of the peak luminosity, where ∆ log L bol ≡ log( L maxbol / L ). The choice of 5 days is arbitrary, but ismotivated by the fact that UVOIR observations of KNe canbe typically performed over a few days. We vary the pa-rameters in the model and discuss the possible variations inthe correlation. In particular, we vary the choice of NRbased fitting formula for the ejecta mass and velocity (pro-vided by Radice et al. 2018; Coughlin et al. 2018), the nuclearEOS (DD2 by Typel et al. 2010 and WFF2 by Wiringa et al.1988), and distribution of the electron fraction Y e of the ejecta( f e = [0 . , . , .
2] and f e = [0 . , . , .
14] used by Villaret al. 2017). We also consider a case where f e depends onthe NS masses . We also examine the same correlations com-puted assuming a time delay of 7 days from peak luminosity.These results are plotted in figure 3, suggesting clear correla-tions between L maxbol and ∆ log L bol .It should be mentioned here that the relation found in thiswork factors in the full non-linearity in the NR simulationsand the non-trivial relationship between NS masses and bolo-metric light curve. The fact that such a correlation has beenobserved in the synthetic light curves suggests that a similarcorrelation should exist in the actual light curves, even thoughthe actual observed correlation may turn out to be differentthan what are presented here. If this is vindicated by futureKNe observations, this will provide an independent distanceladder. For example, in figure 3, the decline in luminosity in5 days (or any other suitably chosen time) is an independentobservable which can be used to find the maximum intrinsicluminosity using the correlation. Thus the luminosity distancecan be estimated by comparing the intrinsic and apparent lu-minosities. SUMMARY AND OUTLOOK
Motivated by the similarities of KNe with SNe Ia, we haveexplored the possibility of KNe providing a set of standardiz-able candles analogous to SNe Ia. Indeed, such a possibilitycan be confirmed or refuted only by a large number KNe ob-servations. As we await such observations, we studied simplesemi-analytical KNe models (in conjunction with NR fittingformulae for ejecta mass and velocity) and discovered cor-relations that exist between the peak bolometric luminosity L maxbol and the decline in the luminosity ∆ log L bol after a fewdays (figure 3). This is performed by computing L maxbol and ∆ log L bol from synthetic light curves generated from ejectaproduced by a population of BNS mergers. We employ dif-ferent NR fitting formulae, NS EOS and electron fraction dis-tribution of the ejecta to study the robustness of our results.We note that Coughlin et al. (2019) has taken a different ap-proach to the standardization of KNe.The light curve calculation presented in this work has mul-tiple simplifying assumptions, and are subject to errors in theNR simulations and the KNe models. Hence they are onlycrude estimates. In particular, anisotropy of the ejecta com-ponents and time-variation of the ejecta opacities needs futureinvestigation. However, there is preliminary evidence (com-ing 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 particularcorrelation, which has to be left to future observations. Sucha correlation in the observed light curves could have potentialusage in distance measurement, which will have applicationsin fundamental physics, astrophysics and cosmology. Possi-ble applications include constraining the number of spacetime We use figure 5 from Dietrich & Ujevic (2017) to get the average electronfraction of the ejecta, ¯ Y e ( q ), as a function of the mass ratio q := M / M of thebinary. Then we solve two constraint equations, (cid:80) m f e m = 1 and (cid:80) f e m Y e m = ¯ Y e ( q ), where m denotes the blue, purple and red components. The system isunder-determined as there are three unknowns ( f e m ) and only two equations.We choose f eblue = 0 . f eblue = 0 . dimensions by comparing distance estimates from GW andKN observation from a BNS merger (along the lines of Pardoet al. 2018; Abbott et al. 2019), the calibration of distanceladders (along the lines of Gupta et al. 2019), and potentiallythe estimation of cosmological parameters (along the lines ofCoughlin et al. 2019).We admit that there are however key differences betweenSN and KNe light curves. KNe have peak luminosities ( ∼ - 10 erg s − ) much lower than SNe Ia peak luminosities( ∼ - 10 erg s − ) making KNe observable only in thelocal universe. The SNe Ia B-band light curve usually peaks ∼
20 days post explosion (Riess et al. 1999) where the spec-tral peak shifts from optical to infrared in about 2-3 months.In contrast, the KNe light curve reaches the maximum valuewithin few hours and shifts from optical to infrared within 10days. Because of these differences, KNe light curves are rel-atively more difficult to standardize and also the counterpartof the “Phillips relation” would be expected to be different forthem. Only future observations can tell the full story.
Acknowledgments: — We thank Kenta Hotokezaka, Bala Iyer,Ramesh Narayan, Sumit Kumar, K. Haris for their patientfeedback to the work. We also thank Bharat Kumar for pro-viding polytropic fits to the EOS DD2, and David Radice andTim Dietrich for clarifications on their NR fits. RK would liketo thank G. Srinivasan for his encouragement to pursue theproblem during a meeting in 2017 at ICTS, Bengaluru. We aregrateful to the anonymous referee for many useful commentsand suggestions on an earlier version of the manuscript. Thisresearch was supported by the Max Planck Society through aMax Planck Partner Group at ICTS-TIFR and by the CanadianInstitute for Advanced Research through the CIFAR AzrieliGlobal Scholars program. Computations were performed atthe ICTS cluster Alice.REFERENCES