Tight multi-messenger constraints on the neutron star equation of state from GW170817 and a forward model for kilonova light curve synthesis
Matt Nicholl, Ben Margalit, Patricia Schmidt, Graham P. Smith, Evan J. Ridley, James Nuttall
MMNRAS , 1–15 (2020) Preprint 5 February 2021 Compiled using MNRAS L A TEX style file v3.0
Tight multi-messenger constraints on the neutron star equation of statefrom GW170817 and a forward model for kilonova light curve synthesis
Matt Nicholl , ★ , Ben Margalit † , Patricia Schmidt , , Graham P. Smith , Evan J. Ridley , & James Nuttall School of Physics and Astronomy, University of Birmingham, Birmingham B15 2TT, UK Institute for Gravitational Wave Astronomy and School of Physics and Astronomy, University of Birmingham, Birmingham B15 2TT, UK Astronomy Department and Theoretical Astrophysics Center, University of California, Berkeley, Berkeley, CA 94720, USA
ABSTRACT
We present a rapid analytic framework for predicting kilonova light curves following neutron star (NS) mergers, where the maininput parameters are binary-based properties measurable by gravitational wave detectors (chirp mass and mass ratio, orbitalinclination) and properties dependent on the nuclear equation of state (tidal deformability, maximum NS mass). This enablessynthesis of a kilonova sample for any NS source population, or determination of the observing depth needed to detect a livekilonova given gravitational wave source parameters in low latency. We validate this code, implemented in the public mosfitpackage, by fitting it to GW170817. A Bayes factor analysis overwhelmingly (
𝐵 > ) favours the inclusion of an additionalluminosity source during the first day, well fit by a shock-heated cocoon, alongside lanthanide-poor dynamical ejecta. Theemission thereafter is dominated by a lanthanide-rich viscous wind. We find the mass ratio of the binary is 𝑞 = . ± .
07 (90%credible interval). We place tight constraints on the maximum stable NS mass, 𝑀 TOV = . + . − . M (cid:12) . For a uniform prior intidal deformability, the radius of a 1.4 M (cid:12) NS is 𝑅 . ∼ . 𝑀 TOV , we derive a final measurement 𝑅 . = . + . − . km. Applying our code to the secondgravitationally-detected neutron star merger, GW190425, we estimate that an associated kilonova would have been fainter (by ∼ . Key words: stars:neutron – neutron star mergers – gravitational waves – methods: data analysis
The first binary neutron star (NS) merger detected through its grav-itational wave emission, GW170817 (Abbott et al. 2017a), wasaccompanied by both a short gamma-ray burst (GRB; Goldstein et al.2017; Savchenko et al. 2017) and an optical counterpart (Coulteret al. 2017; Soares-Santos et al. 2017; Valenti et al. 2017; Arcavi et al.2017; Tanvir et al. 2017; Lipunov et al. 2017). Thermal transientsfrom these mergers had long been predicted in the form of a kilonova(Li & Paczyński 1998; Rosswog et al. 1999; Metzger et al. 2010): anon-relativistic outflow heated by the decays of heavy nuclei formedthrough rapid neutron captures (the so-called ‘r-process’; Lattimer &Schramm 1974; Eichler et al. 1989; Davies et al. 1994; Freiburghauset al. 1999). Extensive optical and infrared follow-up of GW170817by many groups showed photometric and spectroscopic propertiesthat were remarkably consistent with the expectations for an r-processkilonova (Andreoni et al. 2017; Arcavi et al. 2017; Chornock et al.2017; Cowperthwaite et al. 2017; Díaz et al. 2017; Drout et al. 2017;Evans et al. 2017; Hu et al. 2017; Kasliwal et al. 2017; McCully et al.2017; Nicholl et al. 2017a; Pian et al. 2017; Pozanenko et al. 2018;Shappee et al. 2017; Smartt et al. 2017; Tanvir et al. 2017; Troja et al.2017; Utsumi et al. 2017; Valenti et al. 2017). ★ Contact e-mail: m.nicholl.1.bham.ac.uk † NASA Einstein Fellow
Studying these sources is crucial for understanding cosmic chemicalevolution. The total ejected mass powering the GW170817 kilonovawas inferred to be ∼ .
05 M (cid:12) , with at least two components ofdifferent compositions and hence opacities (e.g. Villar et al. 2017;Kasen et al. 2017; Margutti & Chornock 2020): one containing light r-process elements (atomic mass number 𝐴 (cid:46) 𝑓 -shells (Kasen et al. 2013; Barnes &Kasen 2013). The large total mass and broad range of mass numberscovering all three r-process abundance peaks (Burbidge et al. 1957)suggested that NS mergers may be the dominant production site ofheavy nuclei in the Universe.Another central goal of multi-messenger astronomy is to determinethe NS equation of state (EoS), i.e. the relation between pressureand density. This provides unique tests of nuclear physics beyond thesaturation density of nuclear matter (e.g. Lattimer & Prakash 2016).Because the EoS determines the NS radius and maximum stable mass,measuring these quantities from astrophysical observation directlyconstrains the EoS. For example, a viable EoS must support theexistence of the most massive known pulsars, with mass (cid:38) (cid:12) (Demorest et al. 2010; Antoniadis et al. 2013; Cromartie et al. 2020).NS mergers provide an excellent laboratory for EoS measurements.The degree of tidal deformation imprinted in the GW signal (Flanagan& Hinderer 2008; Damour et al. 2012) and the strength of shocks atthe contact interface (Bauswein et al. 2013b; Hotokezaka et al. 2013) © a r X i v : . [ a s t r o - ph . H E ] F e b M. Nicholl et al both depend on the radius, typically parameterised as 𝑅 . , the radiusof a 1.4 M (cid:12) NS. Similarly, the lifetime of the massive merger productuntil collapse to a black hole (BH) depends on the maximum stableNS mass, the Tolman-Oppenheimer-Volkoff mass 𝑀 TOV (Shapiro &Teukolsky 1986; Margalit & Metzger 2017).So far, most of the EoS constraints for GW170817 have come fromthe GW data (e.g. Abbott et al. 2018; De et al. 2018; Raithel et al. 2018).However, electromagnetic data have also provided important insight.Nicholl et al. (2017a) pointed out that the inferred combinationof large mass and high velocity in the low-lanthanide ejecta waspredicted only in the case of relatively small radii (Bauswein et al.2013b; Shibata & Hotokezaka 2019). More quantitatively, Margalit& Metzger (2017) used the observed electromagnetic signatures ofGW170817 to infer that the merged remnant must have collapsedto a BH within a short timescale ( ∼
10 ms to ∼ 𝑀 TOV (cid:46) .
17 M (cid:12) . Bauswein et al. (2017) related this lifetime to theradius, finding 𝑅 . > .
68 km.A major caveat in these analyses is our ability to associate ameasured ejecta component with a specific physical origin. NSmergers can eject mass dynamically during the merger, through tidalstripping or shock heating, and after the merger, via winds from anaccretion disk or long-lived remnant (see review by Metzger 2019).Each mechanism is expected to produce a different combination ofejecta mass, velocity and composition, and the relative importanceof each mechanism depends on the properties of the system beforemerger (e.g. Bauswein et al. 2013b; Hotokezaka et al. 2013; Sekiguchiet al. 2016; Radice et al. 2018; Shibata & Hotokezaka 2019; Ciolfi &Kalinani 2020).Most model fits to the electromagnetic data from GW170817 havebeen specified in terms of the ‘post-merger’ parameters: i.e. the indi-vidual masses, velocities and compositions/opacities of some numberof ejecta components. This introduces subjectivity into the originof each component, and neglects physically-expected correlationsbetween the properties of the different components. Specifying in-stead a forward model in terms of ‘pre-merger’ properties – binarymasses and mass ratios, and the EoS – would allow one to traceexplicitly the ejection mechanism of each component, and excludemodels that fit the data well but require unphysical combinations ofejecta parameters. Moreover, analysing the data in terms of pre-mergerproperties provides a more direct link between the electromagneticand GW signals.Progress in combining light curve models with GW inference wasfirst made by Coughlin et al. (2019a), who fit the light curve of theGW170817 kilonova and related the implied ejecta masses to thepredictions of numerical relativity simulations for different systemmass ratios and combinations of 𝑅 . and 𝑀 TOV . Including GWinformation in their priors, they measured 𝑅 . = . − . 𝑅 . = . ± .
86 km. An even more recent analysisof GW170817 was published by Breschi et al. (2021), who found 𝑅 . = . ± . 𝑅 . and 𝑀 TOV .The structure of this paper is as follows. In section 2, we describethe model setup and assumptions. Section 3 shows its application toGW170817 and our main results. We briefly demonstrate the utilityof our code for aiding future kilonova searches, using the example ofGW190425, in section 4, before concluding in section 5.
We have built our model using the Modular Open Source Fitter forTransients (mosfit; Guillochon et al. 2018). Modules to calculatea radioactive kilonova light curve for a given set of ejecta masses,velocities and opacities have already been implemented by Cow-perthwaite et al. (2017) and Villar et al. (2017). Therefore our maintask here is to define modules that provide these masses, velocitiesand opacities for a given binary configuration. We also provide newmodules to take into account the effects of observer viewing angle andshock cooling. We have made these new modules and the parameter
MNRAS , 1–15 (2020) ilonova light curves from binary parameters Observer q = M /M < 1 Λ = f(q, R NS ) ~ M < M TOV θ 𝓜 = (M M ) (M + M ) Observer
GW EM M red ← q, 𝓜 M blue ← q, R NS θ c θ Shock cooling M purp ↑ M disk ← 𝓜 , M TOV , R NS GRB jet (gamma rays)Dynamical ejecta: Post-merger disk/ejecta:v b/r/p ← 𝓜 Figure 1.
Schematic showing model geometry and relations between parameters. A full description of these parameters is given in section 2. Left: Pre-mergerproperties of a NS binary. Bold, green font is used to indicate properties that can be measured from GW waveforms. Right: Ejecta properties probed byelectromagnetic observations. Left arrows are used to mean ‘depends on’, in order to highlight the most important parameter sensitivities. The merged remnantcan be a long-lived or short-lived NS or a BH, depending on the ratio ( 𝑀 + 𝑀 )/ 𝑀 TOV . The blue ejecta are from shocks at the collisional interface and have anopacity 𝜅 blue ∼ . g − , the red ejecta are tidally stripped in the orbital plane and have 𝜅 red ∼
10 cm g − , and the purple ejecta are winds from an accretiondisk and have an intermediate opacity that depends on the lifetime of the remnant. Low opacity ejecta may not be visible for large viewing angle 𝜃 . The shockcooling emission arises from a cocoon of width 𝜃 𝑐 heated at its interface with the GRB jet. files to implement them publicly available via GitHub , as a newmodel called ‘bns’. In this section we will describe the main features.Figure 1 provides a visual overview that summarises the geometryand the relations between key parameters. The connection betweenour modelling framework and the GW observables is discussed indetail in section 2.6. During the merger process, mass is ejected dynamically both by tidalforces in the orbital plane (Rosswog et al. 1999; Sekiguchi et al. 2015)and by shocks at the contact interface between the two NSs (Bausweinet al. 2013b; Hotokezaka et al. 2013). To estimate the mass of thismaterial, we use a fit by Dietrich & Ujevic (2017) to the ejecta massesof 172 numerical relativity simulations, covering a range of binarymasses and equations of state . Their function takes the form 𝑀 dyn − M (cid:12) = (cid:34) 𝑎 (cid:18) 𝑀 𝑀 (cid:19) / (cid:18) − 𝐶 𝐶 (cid:19) + 𝑏 (cid:18) 𝑀 𝑀 (cid:19) 𝑛 + 𝑐 (cid:0) − 𝑀 / 𝑀 ∗ (cid:1) (cid:35) 𝑀 ∗ + ( ↔ ) + 𝑑, (1)where 𝑀 𝑖 is the mass of the 𝑖 th NS and 𝐶 𝑖 ≡ 𝐺 𝑀 𝑖 / 𝑅 𝑖 𝑐 (2) https://github.com/guillochon/MOSFiT Note however that none of these simulations included the effect of NS spin. is its compactness given a radius 𝑅 𝑖 . The term 1 ↔ 𝑎, 𝑏, 𝑐, 𝑑, 𝑛 can be found in their study. Quantities denoted 𝑀 ∗ 𝑖 arethe baryonic masses (differing from the gravitational mass by thebinding energy), which we calculate for a given 𝑀 𝑖 (in solar masses)following Gao et al. (2020): 𝑀 ∗ 𝑖 = 𝑀 𝑖 + . 𝑀 𝑖 . The typical fractionaluncertainty on 𝑀 dyn using this parameterisation was found by Dietrich& Ujevic (2017) to be ≈ − 𝑀 𝑖 / 𝑀 𝑗 . We adopt the conventionthat 𝑀 is the heavier NS (for consistency with GW analyses, e.g.Abbott et al. 2017a, 2018), such that the mass ratio of the systemis 𝑞 = 𝑀 / 𝑀 <
1. The other important sensitivity in equation 1is to the compactness 𝐶 𝑖 , an EoS-dependent quantity. This will bediscussed further in section 2.6.The relationship between binary and ejecta masses for NS mergerswas recently revisited by Krüger & Foucart (2020) and Nedora et al.(2020), who used a different set of numerical data but found broadlyconsistent fits. Coughlin et al. (2019a) provided a fitting formulain terms of log ( 𝑀 dyn ) , which performed similarly well. However,we find this logarithmic form can lead to unrealistically large ejectamasses in some regions of parameter space, which is not optimal forpopulation synthesis.To better capture the multiple physical processes contributing to thedynamical ejecta, we divide it into two components. The tidal ejecta,concentrated in the equatorial plane, are expected to contain a very low MNRAS000
1. The other important sensitivity in equation 1is to the compactness 𝐶 𝑖 , an EoS-dependent quantity. This will bediscussed further in section 2.6.The relationship between binary and ejecta masses for NS mergerswas recently revisited by Krüger & Foucart (2020) and Nedora et al.(2020), who used a different set of numerical data but found broadlyconsistent fits. Coughlin et al. (2019a) provided a fitting formulain terms of log ( 𝑀 dyn ) , which performed similarly well. However,we find this logarithmic form can lead to unrealistically large ejectamasses in some regions of parameter space, which is not optimal forpopulation synthesis.To better capture the multiple physical processes contributing to thedynamical ejecta, we divide it into two components. The tidal ejecta,concentrated in the equatorial plane, are expected to contain a very low MNRAS000 , 1–15 (2020)
M. Nicholl et al electron fraction (the ratio of protons to total nucleons), 𝑌 𝑒 < . 𝑌 𝑒 > .
25, due to the importance of shocksand neutrino heating at higher latitudes (e.g. Bauswein et al. 2013b;Sekiguchi et al. 2016; Metzger 2019), and is correspondingly ‘blue’due to an incomplete r-process only up to mass number 𝐴 ∼ 𝑌 𝑒 distributions in the dynamical ejecta as a function of 𝑞 for twoequations of state, one stiff (large radius) and one soft (compact).We integrate these distributions (in their Figure 5) to determine 𝑓 red ( 𝑞 ) ≡ 𝑀 ( 𝑌 𝑒 < . )/ 𝑀 ( total ) and 𝑓 blue ≡ − 𝑓 red . We find that 𝑓 red is insensitive to the EoS but is a steep function of 𝑞 , rangingfrom 20% for an equal-mass merger ( 𝑞 =
1) to 76% for 𝑞 = .
86. For 𝑞 (cid:46) . 𝑓 red ≈
1, because the less massive NS is tidally disruptedprior to contact, preventing strong shocks (e.g. Bauswein et al. 2013b;Hotokezaka et al. 2013; Lehner et al. 2016; Dietrich & Ujevic 2017).For 𝑞 > .
8, we interpolate 𝑓 red between the available simulationsusing a smooth polynomial fit (shown in appendix Figure A1). Thepredicted ejecta masses in each component, for a range of NS binaries,are shown in Figure 2. Following Radice et al. (2018), we assumea fixed grey opacity for each component: 𝜅 blue = . g − , and 𝜅 red =
10 cm g − (Villar et al. 2017; Metzger 2019; Tanaka et al.2020). For most systems, additional mass will be ejected after the merger,driven by neutrino-heated winds from a viscous accretion disk aroundthe remnant or from the surface of a massive remnant NS (Metzger &Fernández 2014). Simulations show that in many cases the mass ofpost-merger ejecta will exceed that of the dynamical ejecta (Ciolfi et al.2017; Siegel & Metzger 2017, 2018; Radice et al. 2018; Fernándezet al. 2019).The properties of the wind ejecta depend on the lifetime of themerger remnant, which may be a stable NS, a supramassive (hy-permassive) NS temporarily supported by (differential) rotation, ormay promptly collapse to a BH. To first order, longer-lived remnantsproduce a larger disk mass with a higher 𝑌 𝑒 (Metzger & Fernández2014; Lippuner et al. 2017).We use the analytic expression for the disk mass from Coughlinet al. (2019a), who determine the values of their coefficients 𝑎, 𝑏, 𝑐, 𝑑 by fitting to simulations from Radice et al. (2018):log ( 𝑀 disk ) = max (cid:18) − , 𝑎 (cid:18) + 𝑏 tanh (cid:18) 𝑐 − 𝑀 tot / 𝑀 thr 𝑑 (cid:19)(cid:19)(cid:19) . (3)Here 𝑀 tot = 𝑀 + 𝑀 is the total binary mass, and 𝑀 thr is thethreshold mass for prompt collapse to a BH, parameterised as 𝑀 thr = (cid:18) . − . 𝑀 TOV 𝑅 . (cid:19) 𝑀 TOV , (4)i.e. proportional to the maximum stable mass of a non-rotating NS andthe compactness (Bauswein et al. 2013a). An alternative fit by Radiceet al. (2018) was parameterised in terms of the tidal deformabilityrather than 𝑀 TOV . Following Coughlin et al. (2019a), we prefer toretain the explicit dependence on 𝑀 TOV so that we can try to constrainit observationally. Recent fits by Nedora et al. (2020) predict disk -4 -3 -2 -1 E j ec t a m a ss ( M ) Blue dynRed dynDisk( = 0.2)M =0.8M =1.1 M =1.4M =1.7 M =2.00.6 0.8 1.0 1.2 1.4 1.6 1.8Chirp mass (M )0.00.10.20.3 E j ec t a v e l o c it y ( v / c ) Blue dynRed dynDiskM =0.8M =1.1 M =1.4M =1.7 M =2.0 Figure 2.
Top: Predicted ejecta mass in each component (blue and reddynamical, and purple disk ejecta) using relations from Dietrich & Ujevic(2017) and Coughlin et al. (2019a), assuming 𝑅 . =
11 km, 𝑀 TOV = . (cid:12) ,and that 20% of the remnant disk is ejected. Each colour corresponds to afixed primary mass, 𝑀 , with varying mass ratio, 0 . ≤ 𝑞 ≤
1. Verticallines indicate 𝑀 = 𝑀 , with lower mass secondaries (less equal binaries) tothe left. Bottom: Predictions for ejecta velocities. Disk ejecta velocities areestimated using results from Metzger & Fernández (2014). masses broadly consistent with the results of Coughlin et al. (2019a).The parameterisation above is based on simulations of systems withclose to equal mass ratio, 𝑞 ≈ 𝑞 .A fraction 𝜖 disk ∼ . − . 𝜖 disk ≈ . 𝜖 disk maychange as a function of disk mass; De & Siegel 2020). Additionalpost-merger ejecta may arise through neutrino- or magnetically-drivenwinds (Radice et al. 2018; Metzger et al. 2018; Ciolfi & Kalinani2020), which we will discuss in section 2.5.To determine the velocity and composition of the disk ejecta, wefirst estimate the remnant properties as a function of 𝑀 tot , 𝑀 TOV and 𝑀 thr using table 3 of Metzger (2019). For 𝑀 tot < 𝑀 TOV , the
MNRAS , 1–15 (2020) ilonova light curves from binary parameters remnant is indefinitely stable and we assume an ejecta velocity of ≈ . 𝑐 based on the simulations of Metzger & Fernández (2014).For 𝑀 tot > 𝑀 thr (prompt collapse), the same simulations predict avelocity of ≈ . 𝑐 . We linearly interpolate between these velocitiesfor intermediate 𝑀 tot . Figure 2 shows the predicted disk ejecta massesand velocities alongside the dynamical ejecta from section 2.1.To estimate the opacity we use the simulations from Lippuneret al. (2017), who present 𝑌 𝑒 distributions in their Figure 2 for diskswith various remnant lifetimes, ranging from indefinitely stable toprompt collapse. We integrate these distributions to determine themass-weighted 𝑌 𝑒 . We then convert this mean 𝑌 𝑒 into an effectivegrey opacity, 𝜅 purple , using fits to the 𝑌 𝑒 − 𝜅 relations from Tanakaet al. (2020) (shown in appendix Figure A2). The inferred rangeof 0 . < (cid:104) 𝑌 𝑒 (cid:105) < .
38 corresponds to an opacity range of 5 . (cid:38) 𝜅 purple /( cm g − ) (cid:38) . A semi-analytic model that captures the essential light curve physics ofkilonovae was developed by Metzger (2019). This model approximatesthe heating rate at time 𝑡 , averaged over all r-process radioactive decaychains, as (cid:164) 𝑄 ( 𝑡 ) ∝ 𝜖 th 𝑀 r 𝑡 − . (Korobkin et al. 2012), where 𝑀 r is ther-process mass and 𝜖 th is a time-dependent thermalisation efficiency(Barnes et al. 2016). A module to calculate (cid:164) 𝑄 ( 𝑡 ) in this way wasimplemented in mosfit by Cowperthwaite et al. (2017) and Villaret al. (2017).To determine the output luminosity, 𝐿 ( 𝑡 ) , the ejecta are assumedto have a constant grey opacity and expand homologously with scalevelocity 𝑣 ej , allowing for an analytic radiative diffusion solution tothe first law of thermodynamics (Arnett 1982). The temperature ofthe photosphere is initially calculated using the Stefan-Boltzmannlaw, with radius 𝑅 = 𝑣 ej 𝑡 , until such time as the cooling ejecta startto recombine. At this point, the photosphere recedes into the ejectawhile maintaining the recombination temperature, 𝑇 rec . The mosfitmodules to solve for the radiative diffusion and the photospherictemperature and radius are described in Nicholl et al. (2017b).Cowperthwaite et al. (2017) and Villar et al. (2017) generalisedthis model to include the sums of multiple ejecta components: theluminosity scale of each component is proportional to its mass, with adiffusion timescale determined also by its velocity and opacity. Theyallowed 𝑇 rec to vary independently for each of 2-3 ejecta components,finding best fit values in the range 1000 K (cid:46) 𝑇 rec (cid:46) 𝑇 rec = Most analytic light curve models assume spherical symmetry. How-ever, in kilonovae the separate ejecta components occupy distinctspatial regions, due to their different ejection mechanisms (as illus-trated in Figure 1). The low-medium 𝑌 𝑒 (i.e. blue and purple) ejectaare expected to escape mostly orthogonal to the orbital plane, whereshock heating and weak interactions lower the free neutron flux (e.g.Wanajo et al. 2014; Goriely et al. 2015; Foucart et al. 2016). Incontrast, the neutron-rich red ejecta escape parallel to the orbital planesince they form tidally. Simulations show that this geometry can be well approximated as a sphere with conical sections of half-openingangle 𝜃 open ∼ ◦ around the poles. The blue/purple ejecta residewithin the conical sections and the red ejecta outside.This results in a viewing angle dependence in the observed lu-minosity, colour and timescale of the kilonova (Kasen et al. 2015;Kawaguchi et al. 2018; Wollaeger et al. 2018; Bulla 2019; Korobkinet al. 2020, e.g.), which can be understood intuitively: from viewingangles close to the equator, emission from the low-lanthanide polarejecta is obscured by the curtain of lanthanide-rich ejecta at lowerlatitudes, meaning that only a slow-evolving red kilonova is observed.Radiative transfer simulations show that this can reduce the peakluminosity by ≈ − 𝜃 and thepolar ejecta by 1 − cos 𝜃 , with the viewing angle 𝜃 a free parameter.Here we implement a new module that treats this more exactly,using the formalism of Darbha & Kasen (2020). This takes a spherewith conical polar caps of half-opening angle 𝜃 open , and scales theluminosity of each component using the projected area of the caps(blue/purple ejecta) or remaining sphere (red ejecta) as a functionof 𝜃 and 𝜃 open . We have verified that the change in luminosity with 𝜃 is consistent with the radiative transfer results of Bulla (2019) forthe same geometry. We note that this treatment does not accountfor the jet-interaction viewing angle effect discussed by Klion et al.(2021), as this would depend on uncertain properties of the jet. Forour fiducial model, we assume 𝜃 open = ◦ and vary the observerviewing angle 𝜃 . We include in our model two optional final ingredients that cancontribute additional luminosity. Although they are tied less explicitlyin our formalism to the initial binary parameters, neglecting thesecontributions may lead to an underestimate in the model luminosityor to a bias in the inference of other parameters from kilonova data.These components can easily be turned off at run-time, allowing theuser to recover the fully predictive baseline model described in theprevious sections.The first is an enhancement in the blue ejecta due to magnetically-driven winds (Metzger et al. 2018; Ciolfi & Kalinani 2020). Thismechanism is possible only if the remnant avoids prompt collapse,hence we apply this effect only when the merger mass is below 𝑀 thr (equation 4). We adopt a simple parameterisation, inspired byCoughlin et al. (2019a): 𝑀 blue , tot = 𝑀 blue , dyn / 𝛼 , where 0 < 𝛼 < 𝛼 = MNRAS000
38 corresponds to an opacity range of 5 . (cid:38) 𝜅 purple /( cm g − ) (cid:38) . A semi-analytic model that captures the essential light curve physics ofkilonovae was developed by Metzger (2019). This model approximatesthe heating rate at time 𝑡 , averaged over all r-process radioactive decaychains, as (cid:164) 𝑄 ( 𝑡 ) ∝ 𝜖 th 𝑀 r 𝑡 − . (Korobkin et al. 2012), where 𝑀 r is ther-process mass and 𝜖 th is a time-dependent thermalisation efficiency(Barnes et al. 2016). A module to calculate (cid:164) 𝑄 ( 𝑡 ) in this way wasimplemented in mosfit by Cowperthwaite et al. (2017) and Villaret al. (2017).To determine the output luminosity, 𝐿 ( 𝑡 ) , the ejecta are assumedto have a constant grey opacity and expand homologously with scalevelocity 𝑣 ej , allowing for an analytic radiative diffusion solution tothe first law of thermodynamics (Arnett 1982). The temperature ofthe photosphere is initially calculated using the Stefan-Boltzmannlaw, with radius 𝑅 = 𝑣 ej 𝑡 , until such time as the cooling ejecta startto recombine. At this point, the photosphere recedes into the ejectawhile maintaining the recombination temperature, 𝑇 rec . The mosfitmodules to solve for the radiative diffusion and the photospherictemperature and radius are described in Nicholl et al. (2017b).Cowperthwaite et al. (2017) and Villar et al. (2017) generalisedthis model to include the sums of multiple ejecta components: theluminosity scale of each component is proportional to its mass, with adiffusion timescale determined also by its velocity and opacity. Theyallowed 𝑇 rec to vary independently for each of 2-3 ejecta components,finding best fit values in the range 1000 K (cid:46) 𝑇 rec (cid:46) 𝑇 rec = Most analytic light curve models assume spherical symmetry. How-ever, in kilonovae the separate ejecta components occupy distinctspatial regions, due to their different ejection mechanisms (as illus-trated in Figure 1). The low-medium 𝑌 𝑒 (i.e. blue and purple) ejectaare expected to escape mostly orthogonal to the orbital plane, whereshock heating and weak interactions lower the free neutron flux (e.g.Wanajo et al. 2014; Goriely et al. 2015; Foucart et al. 2016). Incontrast, the neutron-rich red ejecta escape parallel to the orbital planesince they form tidally. Simulations show that this geometry can be well approximated as a sphere with conical sections of half-openingangle 𝜃 open ∼ ◦ around the poles. The blue/purple ejecta residewithin the conical sections and the red ejecta outside.This results in a viewing angle dependence in the observed lu-minosity, colour and timescale of the kilonova (Kasen et al. 2015;Kawaguchi et al. 2018; Wollaeger et al. 2018; Bulla 2019; Korobkinet al. 2020, e.g.), which can be understood intuitively: from viewingangles close to the equator, emission from the low-lanthanide polarejecta is obscured by the curtain of lanthanide-rich ejecta at lowerlatitudes, meaning that only a slow-evolving red kilonova is observed.Radiative transfer simulations show that this can reduce the peakluminosity by ≈ − 𝜃 and thepolar ejecta by 1 − cos 𝜃 , with the viewing angle 𝜃 a free parameter.Here we implement a new module that treats this more exactly,using the formalism of Darbha & Kasen (2020). This takes a spherewith conical polar caps of half-opening angle 𝜃 open , and scales theluminosity of each component using the projected area of the caps(blue/purple ejecta) or remaining sphere (red ejecta) as a functionof 𝜃 and 𝜃 open . We have verified that the change in luminosity with 𝜃 is consistent with the radiative transfer results of Bulla (2019) forthe same geometry. We note that this treatment does not accountfor the jet-interaction viewing angle effect discussed by Klion et al.(2021), as this would depend on uncertain properties of the jet. Forour fiducial model, we assume 𝜃 open = ◦ and vary the observerviewing angle 𝜃 . We include in our model two optional final ingredients that cancontribute additional luminosity. Although they are tied less explicitlyin our formalism to the initial binary parameters, neglecting thesecontributions may lead to an underestimate in the model luminosityor to a bias in the inference of other parameters from kilonova data.These components can easily be turned off at run-time, allowing theuser to recover the fully predictive baseline model described in theprevious sections.The first is an enhancement in the blue ejecta due to magnetically-driven winds (Metzger et al. 2018; Ciolfi & Kalinani 2020). Thismechanism is possible only if the remnant avoids prompt collapse,hence we apply this effect only when the merger mass is below 𝑀 thr (equation 4). We adopt a simple parameterisation, inspired byCoughlin et al. (2019a): 𝑀 blue , tot = 𝑀 blue , dyn / 𝛼 , where 0 < 𝛼 < 𝛼 = MNRAS000 , 1–15 (2020)
M. Nicholl et al A b s o l u t e m a gn it ud e r Ku M c = 1.2Mq = 1.0= 0Kasen+ 17 0 10 20 30Time (days)161412108 A b s o l u t e m a gn it ud e q = 0.850 10 20 30Time (days)161412108 A b s o l u t e m a gn it ud e = 60 0 10 20 30Time (days)161412108 A b s o l u t e m a gn it ud e With 30 cocoon
Figure 3.
Example light curves from our model in the near-UV ( 𝑢 ), optical ( 𝑟 ) and near-infrared ( 𝐾 ) bands, for a merger with chirp mass M = . (cid:12) . Top left:Fiducial, equal mass merger viewed along pole. Top right: Merger with same chirp mass but asymmetric mass ratio, 𝑞 = .
85. Lower left: Fiducial model viewedat 60 ◦ off-axis. Lower right: Including shock-heated cocoon with opening angle 30 ◦ . The top two panels show light curves calculated with the models of Kasenet al. (2017) for the closest available ejecta masses (assuming lanthanide fractions of 10 − , 10 − and 10 − for the blue, purple and red components, respectively).These models have not been tuned to match our simulations and are intended only to guide the eye. The viewing angle effect is very roughly approximated byadding − . 𝑟 band and + . 𝐾 band. by the speed of light (Piro & Kollmeier 2018). Following Nakar &Piran (2017), we assume that the shock deposits constant energy perdecade of velocity in the ejecta. The isotropic-equivalent luminosityis proportional to the mass of shocked ejecta, which is a fraction 𝜃 / 𝜃 c equal to zero (or cos 𝜃 c = 𝑀 ej and 𝜃 c .We show examples of light curves calculated for different massratios, inclinations and cocoon opening angles in Figure 3. After implementing the modules described above, our model dependson the masses of the constituent neutron stars, their compactness, themaximum stable mass 𝑀 TOV , and the viewing angle 𝜃 . The fraction ofthe disk mass ejected in the purple component can be varied or fixedto a fiducial value. Two further parameters allow one to increase the early blue flux through magnetic winds and/or shock cooling. Finally,we include line-of-sight extinction, and a white noise parameter 𝜎 inour likelihood function when fitting to data (see Nicholl et al. 2017b;Villar et al. 2017; Guillochon et al. 2018).There is a choice to be made about how best to express some ofthese parameter dependencies. We opt to use the parameterisation thatmost closely matches the quantities constrained by GW observations.This allows the user to quickly predict the electromagnetic lightcurve for a given GW signal, or to use GW information to informtheir priors when fitting to an observed multi-messenger source. Wetherefore parameterise the masses of the system in terms of the ‘chirp’mass M = ( 𝑀 𝑀 ) / ( 𝑀 + 𝑀 ) − / , the quantity most accuratelymeasured from the GW signal, and 𝑞 = 𝑀 / 𝑀 ≤
1, to which thissignal is also somewhat sensitive.The dependence on the NS compactness 𝐶 is more complicated. Forpopulation synthesis, it is simplest to specify a radius 𝑅 , equivalent toimplicitly choosing an equation of state (since all NSs of comparablemass should have approximately the same radius). However, whenfitting to observed data (i.e. when we want to measure, rather thanimpose, a NS radius), it is more useful to express 𝐶 in terms of the NStidal deformability, which measures the responsiveness of the NS to MNRAS , 1–15 (2020) ilonova light curves from binary parameters Table 1.
Parameters in the mosfit bns model and application to GW170817. Some parameters are relevant only when fitting to data. Priors in square brackets areflat over the specified ranges, otherwise Gaussian. Posterior values for tidal deformability are given in terms of ˜ Λ to facilitate comparison to GW analyses. The‘Surface ejecta’ model includes additional blue ejecta from a long-lived remnant; the ‘Shocked cocoon’ model includes emission from a GRB-heated cocoon; the‘Agnostic’ model includes both. The stated best-fit values and uncertainties correspond to the medians and 16th/84th percentiles of the marginalised posteriors.Prior distribution Posterior distribution for fit to GW170817Parameter Default prior GW170817 prior Base model Surface ejecta Shocked cocoon Agnostic M 𝑎 (M (cid:12) ) [ . , . ] . ± .
002 1 . ± .
002 1 . ± .
002 1 . ± .
002 1 . ± . 𝑞 𝑏 [ . , . ] [ . , . ] . ± .
01 0 . + . − . . + . − . . ± . 𝜃 𝑐 [ . , . ] [ . , . ] . ± .
04 0 . + . − . . + . − . . + . − . 𝑀 TOV 𝑑 (M (cid:12) ) [ . , . ] [ . , . ] . ± .
04 2 . ± .
03 2 . + . − . . + . − . Λ s 𝑒 [ , ] [ , ] ∗ ˜ Λ = + − ˜ Λ = + − ˜ Λ = + − ˜ Λ = + − 𝜖 disk 𝑓 [ . , . ] [ . , . ] . ± .
01 0 . ± .
01 0 . ± .
01 0 . ± . 𝛼 𝑔 [ . , . ] [ . , . ] — 0 . + . − . — 0 . + . − . cos 𝜃 c ℎ [ . , . ] [ . , . ] — — 0 . + . − . . + . − . 𝑡 𝑖 (days) [− . , . ] ( − . 𝑁 H 𝑗 [ . , . ] [ . , . ] . + . − . . + . − . . + . − . . + . − . log 𝜎 𝑘 — [− . , . ] − . ± . − . ± . − . ± . − . ± . 𝑍 ): 37.1 56.9 80.3 81.1 𝑎 Chirp mass; 𝑏 Mass ratio; 𝑐 Observer viewing angle; 𝑑 Maximum stable NS mass 𝑒 Symmetric tidal deformability (specified instead as 𝑅 . in generativemode); 𝑓 Fraction of disk ejected; 𝑔 Enhancement (1 / 𝛼 ) of blue ejecta by surface winds; ℎ Opening angle of shocked cocoon; 𝑖 Time of explosion (rest-frame); 𝑗 Hydrogen column density in host galaxy (proportional to extinction); 𝑘 White noise in likelihood function; ∗ Additional constraint that derived parameter˜ Λ < an external tidal field: Λ 𝑖 = ( / ) 𝑘 𝐶 − 𝑖 , where 𝑘 ( 𝑀 ) ≈ . − . Λ and Λ (Flanagan &Hinderer 2008; Wade et al. 2014; Favata 2014),˜ Λ = ( 𝑚 + 𝑀 ) 𝑀 Λ + ( 𝑀 + 𝑀 ) 𝑀 Λ ( 𝑀 + 𝑀 ) , (5)called the binary or effective tidal deformability, which can thereforebe constrained with GW data. The relation between Λ and 𝐶 allowsus to include GW information in our priors to better measure 𝐶 andhence 𝑅 . A similar approach was adopted by Coughlin et al. (2019a),though they used a different relation proposed by De et al. (2018) torelate ˜ Λ to 𝑅 without calculating 𝐶 explicitly.For consistency with GW analyses (Abbott et al. 2018) and toreduce systematic error, we instead sample Λ in terms of the so-called‘symmetric’ tidal deformability: Λ s ≡ ( Λ + Λ )/ . (6)Using ‘quasi-universal relations’ or QURs (largely independent of theEoS, and hence of complicating factors such as 𝑘 ), the antisymmetrictidal deformability Λ a ≡ ( Λ − Λ )/ Λ s to ≈
5% precision (Yagi & Yunes 2016), allowing one to reconstruct Λ , Λ and ˜ Λ directly from Λ s . One may then impose the constraint that˜ Λ be consistent with the GW results, using the constraints class inmosfit (see Guillochon et al. 2018). The compactness of each NS isthen derived from Λ 𝑖 using another QUR: 𝐶 𝑖 = . − . ( Λ 𝑖 ) + . ( Λ 𝑖 ) (7)This relation is accurate to ≈ .
5% across a variety of equations ofstate (Yagi & Yunes 2017). These 𝐶 𝑖 are then used in equation 1to estimate the dynamical ejecta mass. Following a fit to observeddata, the NS radii can be derived from the posterior distribution of Λ s following the same procedure and using the constituent NS masses inequation 2.Our final parameter set therefore consists of 5-11 parameters, whichdiffer slightly between the generative and fitting versions of the model.We list all of these parameters and suggest priors in Table 1. For comparison, the parameterisation of the original mosfit kilonovamodel by Villar et al. (2017) also had 11 free parameters. As with any model of NS mergers, GW170817 provides an idealtesting ground for the formalism described above. In this section, wefit the observed optical data to demonstrate (i) how using a binary(rather than ejecta) based parameter set provides additional insightinto the physical origin of each luminosity component; (ii) the utilityof this model in probing the NS EoS with kilonova data; and (iii)the extent to which the use (or not) of GW information in the modelpriors affects our posteriors.We use data from the Open Kilonova Catalog (Guillochon et al.2017), compiled by Villar et al. (2017). Given the high level of overlapbetween the observations from different groups, we fit the followingsubset, chosen to cover the full range of bands with a comparable(approximately nightly) density of sampling in each: 𝑢𝑔𝑟𝑖𝑧𝑦𝐻 fromCowperthwaite et al. (2017),
𝐵𝑉𝑔𝑟𝐽𝐻 from Drout et al. (2017), 𝑟𝑖𝑧𝑦 from Smartt et al. (2017), 𝑟𝑧𝐾 from Tanvir et al. (2017),
𝐽𝐻𝐾 fromKasliwal et al. (2017), 𝐵𝑉 from Troja et al. (2017), and the 𝑈 andUV bands from Evans et al. (2017).Our priors are given in Table 1. The time of merger is fixedby the GW detection. The priors on M and 𝑞 are from Abbottet al. (2017a), while the prior on Λ s is chosen to match the analysisemployed by Abbott et al. (2018). The range of viewing angles,0 . < cos 𝜃 < .
97, is from Nakar & Piran (2020), who compiledresults from GRB afterglow modelling (Alexander et al. 2017, 2018;Haggard et al. 2017; Margutti et al. 2017, 2018b; Troja et al. 2017,2018, 2019; D’Avanzo et al. 2018; Dobie et al. 2018; Gill & Granot2018; Granot et al. 2018; Lazzati et al. 2018; Lyman et al. 2018;Mooley et al. 2018; Fong et al. 2019; Hajela et al. 2019; Lamb et al.2019; Wu & MacFadyen 2019; Ryan et al. 2020) and from very longbaseline interferometry (VLBI; Ghirlanda et al. 2019; Hotokezakaet al. 2019). We assume a distance of 40.7 Mpc (Cantiello et al. 2018)to the host galaxy, NGC 4993. We run our fits on the University ofBirmingham bluebear cluster, and use dynesty (Speagle 2020) tosample the posteriors of the model parameter space.
MNRAS000
MNRAS000 , 1–15 (2020)
M. Nicholl et al -1 Days since merger12151821242730 A pp a r e n t m a gn it ud e + c on s t a n t W2+7M2+5W1+5 U+5u+4B+3 g+2V+1r+0 i-1z-2y-3 J-4H-5K-6 -1 Days since merger12151821242730 A pp a r e n t m a gn it ud e + c on s t a n t W2+7M2+5W1+5 U+5u+4B+3 g+2V+1r+0 i-1z-2y-3 J-4H-5K-6
Figure 4.
Left: Best-fit kilonova model for GW170817 without surface ejecta enhancement or GRB shock cooling. Right: Best-fit model with these effectsincluded. Both models provide a good match to the data beyond 𝑡 (cid:38) . ≈ . 𝐵 > . Light curve fits for the surface ejectaand shock-only models are shown in the appendix. We begin by fitting GW170817 with four variants of the model todetermine the relative importance of different effects. These are (i)the baseline model ( 𝛼 = cos 𝜃 c = 𝛼 ≤ 𝜃 c ≤ 𝛼, cos 𝜃 c ≤ 𝐵 ≡ 𝑍 / 𝑍 where 𝑍 𝑖 is the Bayesianevidence for model 𝑖 , to compare these models. Speagle (2020) givean extensive account of how 𝑍 is evaluated by dynesty (and a goodoverview of Bayes’ Theorem); we refer the reader there for details.Generally, 𝐵 >
10 is taken to indicate a strong preference, based onthe available data, for one model over another, and
𝐵 >
100 to indicatea decisive preference. Compared to the base model, the model withenhanced surface wind ejecta is preferred with
𝐵 > .However, our fits overwhelmingly prefer models that include shockcooling compared to either of the other models, with Bayes factors 𝐵 > compared to the surface ejecta model and 𝐵 > compared to the base model. The agnostic model (with both surfacewinds and shocks) is weakly favoured over the shock-only model,with 𝐵 = .
2. The decisive preference for shock cooling can beunderstood from Figure 4. Models without shocks, generating theirearly luminosity only from radioactive decay in the blue ejecta, aretoo faint by ≈ . 𝑡 ∼ . 𝛼 improves the fit at 𝑡 ∼ . 𝜃 c ∼ − ◦ provides the excess luminosity at 𝑡 = . Taking the agnostic version as our preferred model, we show thetwo-dimensional posteriors for this fit in Figure 5. As with all modelvariants considered, the chirp mass posterior is essentially the prior,since this parameter is constrained to such high precision by the GWdata. The electromagnetic data tighten the mass ratio of the systemto 𝑞 > .
84 (about 50% of the GW-inferred prior volume) due tothe requirement for blue ejecta. Viewing angles close to ∼ ◦ arepreferred, more in line with afterglow modelling than with VLBI(Hajela et al. 2019; Nakar & Piran 2020). This is also consistent withthe GW-inferred viewing angle at the distance of the GW170817 host(Finstad et al. 2018).The largest degeneracies in the model parameters are between themass ratio and the surface ejecta enhancement 1 / 𝛼 , since equal massbinaries also produce more blue ejecta; between 𝑀 TOV and 𝜖 disk ,since both affect the mass of purple disk ejecta; and between 𝜃 c and˜ Λ , with a larger shock heating contribution to the luminosity relaxingthe preference for low tidal deformability (compact EoS). We willdiscuss the EoS-dependent quantities, 𝑀 TOV and 𝑅 . ( ˜ Λ ), in moredetail in section 3.3.The best-fitting models produce typical ejecta masses 𝑀 blue ≈ MNRAS , 1–15 (2020) ilonova light curves from binary parameters M chirp (M ) = 1.19 +0.00-0.00 . . . . q = M / M q=M /M = 0.92 +0.05-0.05 . . . . . c o s cos = 0.84 +0.03-0.02 . . . d i s k disk = 0.12 +0.01-0.01 . . . . = 0.63 +0.18-0.19 = 231.29 +92.49-73.63 . . . . M T O V ( M ) M TOV (M ) = 2.17 +0.05-0.06 . . . . l og A h o s t V ( m ag ) log A hostV ( mag ) = -1.59 +0.39-0.41 . . . . . M chirp (M ) . . . . . l og .
84 0 .
88 0 .
92 0 . q = M /M . . . . . cos .
105 0 .
120 0 . disk . . . . .
08 2 .
16 2 .
24 2 . M TOV (M ) . . . . log A hostV ( mag ) .
52 0 .
48 0 .
44 0 .
40 0 . log log = -0.44 +0.02-0.02 Figure 5.
Posterior distributions of free parameters for the preferred agnostic model (including shock cooling and enhanced surface wind ejecta) fit to GW170817,using multi-messenger priors (Table 1). .
005 M (cid:12) , 𝑀 red ≈ .
001 M (cid:12) , and 𝑀 purple ≈ .
02 M (cid:12) . As we haveestablished, the blue ejecta can only account for ∼
50% of the flux at 𝑡 ∼ . 𝑡 ∼ −
10 days, with the low mass of red ejecta contributingsignificantly only in the extended tail. These results are in agreementwith Villar et al. (2017); however in our forward model we canuniquely associate the purple component with remnant disk winds,the blue component with shock-driven dynamical ejecta (possiblybut not necessarily boosted by magnetic surface winds), and the red component with tidal dynamical ejecta, while identifying theimportance of the cocoon contribution to the luminosity at early times.
In this section, we examine the posteriors for EoS-dependent quantitiesin our fit, including estimates of the systematic uncertainty. While 𝑀 TOV is treated explicitly as a model parameter, the NS radii mustbe derived from the posterior of the symmetric tidal deformability,
MNRAS000
MNRAS000 , 1–15 (2020) M. Nicholl et al p ( R a d i u s ) -only: R = 10.72 +1.67-1.48 km PrimarySecondaryA20188 10 12 14 16Radius (km) p ( R a d i u s ) + M TOV : R = 11.06 +1.01-0.98 km PrimarySecondaryA2018-EoSCo2019Ca2020D2020B20218 10 12 14 16 18Radius (km)0.00.51.01.52.02.53.0 M a ss ( M ) Radius prior weights from M
TOV
10 15Radius (M=1.4) (km)=11.55, =0.75
Figure 6.
Posterior distributions for the NS radii using the agnostic model,including systematic errors. Vertical lines mark the 5th, 50th and 95th per-centiles for the primary (more massive) NS. Plotted for comparison are themeasurements by Abbott et al. (2018) (with and without imposing constraintson the EoS), Coughlin et al. (2019a), Capano et al. (2020), Dietrich et al.(2020) and Breschi et al. (2021). Top: NS radii from the posterior of Λ s .Middle: NS radii after re-weighting using equations of state that support an 𝑀 TOV within our 90% confidence interval. Bottom: construction of theseprior weights using mass-radius curves from Dietrich et al. (2020). The redcurves satisfy our constraint on 𝑀 TOV . Inset: Gaussian fit to allowed radii. Λ s = + − , and the masses of the two NSs, obtained from M and 𝑞 . The major source of systematic uncertainty in 𝑅 . and 𝑀 TOV isthe scatter in relations 1 and 3, which are respectively ∼
70% for thedynamical ejecta mass (Dietrich & Ujevic 2017) and ∼
50% of ourinferred disk mass (Coughlin et al. 2019a). We test two methods tosimultaneously model these systematic errors. First we run a suite offits where equations 1 and 3 are modified each time by a random factordrawn from a Gaussian distribution with mean equal to 1 and widthequal to the calibration uncertainties. Comparing the posteriors of 50such runs, we find a scatter of ±
40 in Λ s , and negligible scatter in 𝑀 TOV . As a complementary method, we run a single realisation of thefit but with two additional free parameters with Gaussian priors, 𝜎 dyn and 𝜎 disk (where the dynamical ejecta mass is 𝑀 (cid:48) dyn = 𝜎 dyn 𝑀 dyn andthe disk ejecta is treated analogously). In the former case, we broadenthe posterior of Λ s using random draws from a Gaussian of width40, whereas in the latter case we can use the posterior distribution of Λ s directly. We find that the two approaches yield indistinguishableposteriors for Λ s (and hence 𝑅 . ) and 𝑀 TOV .Since the systematic error is negligible for 𝑀 TOV , our estimate ofthis quantity is equal to the posterior distribution shown in Figure 5.We measure a 90% credible interval 𝑀 TOV = . + . − . M (cid:12) . Modelswithout shock cooling (strongly disfavoured by the Bayes factor anal-ysis) prefer a more massive 𝑀 TOV ∼ . (cid:12) . For comparison, theanalysis of Margalit & Metzger (2017) found that 𝑀 TOV (cid:46) . (cid:12) ,which is consistent with our preferred shock cooling models butin potential tension with the shock-free models (although Shibataet al. 2019 have extended the Margalit & Metzger 2017 analysisand found 𝑀 TOV (cid:46) . (cid:12) ). It is worth noting that the 𝑀 TOV con-straints obtained here are complementary to the Margalit & Metzger(2017) constraints, which were based on energetic arguments that areindependent from the kilonova modeling approach adopted in ourpresent work. It is therefore interesting (and non-trivial) that the twoapproaches yield consistent results.We next derive 𝑅 . from the posteriors of Λ s , M and 𝑞 , usingequations 2, 6 and 7. During this conversion we include the systematicuncertainty introduced by the QURs (Yagi & Yunes 2016, 2017) byadding random scatter to the derived values for 𝐶 𝑖 (6.5%) and Λ a (5%), drawn from Gaussian distributions. The masses of the primaryand secondary, which enter equation 2, are 𝑀 = . ± .
04 M (cid:12) and 𝑀 = . ± .
03 M (cid:12) ; therefore their measured radii (especially forthe primary) can be taken as very close to 𝑅 . .The top panel of Figure 6 shows the derived probability densityfunctions for the NS radii evaluated in this way. We find the primaryhas a radius 𝑅 = . + . − . km (90% credible interval). The modelswithout shock cooling prefer a smaller radius, ∼
10 km, due to theirlower values for Λ s . The most direct comparison for our radius is withthe GW-inferred value; Abbott et al. (2018) find 𝑅 GW = . + . − . km(without imposing any EoS constraints a priori ). Our measurement isin excellent agreement with this value, and has a 90% credible boundthat is smaller by ∼ . 𝑅 GW = . ± . 𝑀 TOV to tightenour posteriors on the NS radii. This assumes that the posteriors of˜ Λ and 𝑀 TOV are uncorrelated (Figure 5 shows this to be a goodapproximation). We download the mass-radius curves for a sample of4000 equations of state from Dietrich et al. (2020), constructed usingchiral effective theory (up to 1 . MNRAS , 1–15 (2020) ilonova light curves from binary parameters that is motivated by nuclear theory (Tews et al. 2018; Capano et al.2020). As shown in the bottom panel of Figure 6, we select only thosethat support a TOV mass within our 90% confidence interval. Theinset histogram shows the distribution of 𝑅 . for this EoS subset.We take this as the a priori probability of a given radius in order tore-weight our posteriors for 𝑅 and 𝑅 (after first transforming to aprior flat in 𝑅 . following Raithel et al. 2020), resulting in the middlepanel of Figure 6. Our final measurement using Λ s and ensuringconsistency with 𝑀 TOV is 𝑅 = . + . − . km.Similar to Abbott et al. (2018), imposing a constraint on the EoSresults in an increase in 𝑅 . and a tightening of the posterior. Our finalvalue is consistent with theirs at around the 1 𝜎 level, and narrowerby ∼ . . + . − . km found byCapano et al. (2020). It is also consistent with the 11 . + . − . km foundby Dietrich et al. (2020), at around the 1 𝜎 level. Breschi et al. (2021)report a slightly larger value for 𝑅 . , finding it to be 12 . + . − . km,which is consistent with our results only at the ∼ 𝜎 level. In this section we emphasise the importance of the multi-messenger constraints we used to fit GW170817. Table 1 lists two sets of priors:one set of suggested default priors for fitting an arbitrary NS merger,and one specifically tuned for GW170817. This latter set of priors,which we employed in the previous section, makes use of the GWconstraints on M , 𝑞 , 𝑡 and ˜ Λ , GRB afterglow constraints on 𝜃 (theGRB also constrains 𝑡 ), and constraints on extinction from analysisof the host galaxy (Blanchard et al. 2017; Levan et al. 2017; Pan et al.2017).Taking our preferred model for GW170817 (the agnostic modelincluding both shock cooling and an enhanced surface wind ejecta),we change the priors to the default set and re-run the fit to GW170817.Figure 7 shows the difference in posteriors when using the broaddefault priors, compared to the results using the multi-messengerpriors. We find a fit of similar quality (ln 𝑍 = . M obtained using the broad priorshas a median M broad = . (cid:12) , close to the known value from theGW data. However, the credible region is ∼
30 times wider than whenusing the very tight GW constraints. Another notable difference is inthe viewing angle, where much larger viewing angles 𝜃 broad ∼ ◦ are preferred without the GRB constraints. The distribution of massratio spans a similar range as in the multi-messenger case, but isskewed towards more equal systems ( 𝑞 broad ∼ . 𝛼 and cos 𝜃 c are very similar in the twocases, as are the posteriors of ˜ Λ . However, the weaker constraints onthe total system mass lead to a broadening of the posterior distributionof 𝑀 TOV , and a skew towards larger values.The takeaway message is that our model can do a reasonable job ofmeasuring important system properties for a GW170817-like mergereven without a GW detection. However, having information fromelectromagnetic, GW, GRB and host galaxy analyses leads to tighterconstraints, and this effect is perhaps most important when attemptingto probe the NS EoS, as in section 3.3.
Having validated our model using GW170817, we now (briefly)demonstrate its utility for generating simulated kilonova data. Whilethe broader applications in terms of population synthesis will beexplored in future works, an important use case that we discuss hereis to predict the electromagnetic observables of a particular kilonovagiven a detected GW signal (see e.g. Margalit & Metzger 2019), inorder to make follow-up observations more efficient. This is especiallyrelevant given the possibility of target-of-opportunity GW follow-upwith the upcoming Vera Rubin Observatory (Margutti et al. 2018a;Cowperthwaite et al. 2019; Smith et al. 2019; Chen et al. 2020).During the LIGO-Virgo O3 run, the GW detectors identified severalbinary NS candidates. One event is contained in the latest GW sourcecatalog, GWTC-2 (Abbott et al. 2020a), with the others discoveredin the second half of O3 such that source parameters have not yetbeen released. The published event is GW190425, which was notablefor its large total mass 𝑀 + 𝑀 = . (cid:12) (Abbott et al. 2020c),implying that (unlike GW170817) at least one of the constituent NSswas more massive than the canonical ≈ . (cid:12) .Various groups of astronomers followed up the GW discovery ofGW190425 with electromagnetic observations in an attempt to detecta counterpart, though without success on this occasion (Hosseinzadehet al. 2019; Coughlin et al. 2019b; Lundquist et al. 2019; Antier et al.2020). This may have been due in part to the large distance to thesource, 159 + − Mpc, and in part because the signal was significantonly in one GW detector, leading to a wide sky localization of > (Abbott et al. 2020c).Hosseinzadeh et al. (2019) compiled detection limits from bothgalaxy-targeted and publicly-reported wide-field searches, and com-pared to the light curves of GW170817 (shifted to the distance ofGW190425). In Figure 8, we show these limits compared to ourbest-fit model of GW170817 in the left panel. In agreement withHosseinzadeh et al. (2019), we see that the deep galaxy-targetedsearches rule out a GW170817-like counterpart (in those particulargalaxies) over most of the 90% plausible distance range.With our new model, we can use the chirp mass and mass ratioconstraints obtained by LIGO and Virgo for GW190425 to predict alight curve more directly applicable to this source. A similar approachwas adopted by Barbieri et al. (2020) for both NS-NS and NS-BHbinaries. This prediction is not only tuned specifically to GW190425,but also folds in the uncertainties on the binary parameters as wellas the distance uncertainty. We use a Gaussian distribution for thechirp mass M = . ± . . < 𝑞 < < cos 𝜃 <
1, 0 . < cos 𝜃 c <
1, and a fiducial 𝜖 disk = .
15 basedon GW170817. The result is shown in the right panel of Figure 8.Two features of this prediction are noteworthy. One is the increasedwidth of the 90% distribution of credible light curves, allowing morespace for a faint kilonova to go undetected at the observed depths.More significantly, the median light curve is fainter (by ∼ . 𝑟 band at 1 day post-merger) and declines more rapidly than GW170817.The main reason for this is that the more massive remnant in thissystem is expected to collapse promptly to a BH, reducing the massof the purple ejecta component by an order of magnitude compared toGW170817. With this faster decline, even many of the galaxy-targetedobservations obtained (cid:38) ∼
50% ofthe plausible light curves.Our goal here is not to calculate precisely the fraction of thelocalization volume probed by optical observations of GW190425,but simply to make the general point that how constraining a detectionlimit is depends on the choice of comparison model, and that this in
MNRAS000
MNRAS000 , 1–15 (2020) M. Nicholl et al M chirp (M ) = 1.20 +0.06-0.07 . . . q = M / M q=M /M = 0.94 +0.04-0.05 . . . . c o s cos = 0.71 +0.06-0.05 . . . . d i s k disk = 0.15 +0.01-0.01 . . . . = 0.64 +0.18-0.19 . . . . c o s c o c oo n cos cocoon = 0.93 +0.04-0.07 = 210.14 +89.52-59.73 . . . . M T O V ( M ) M TOV (M ) = 2.32 +0.11-0.15 . . . . . l og A h o s t V ( m ag ) log A hostV ( mag ) = -1.74 +0.40-0.33 .
08 1 .
14 1 .
20 1 .
26 1 . M chirp (M ) . . . . l og .
88 0 .
92 0 . q=M /M .
64 0 .
72 0 .
80 0 . cos .
12 0 .
14 0 .
16 0 . disk . . . . .
78 0 .
84 0 .
90 0 . cos cocoon
100 200 300 400 500 2 . . . . M TOV (M ) . . . . . log A hostV ( mag ) .
48 0 .
45 0 .
42 0 . log log = -0.45 +0.02-0.02 Multimessenger-priorsBroad-priors
Figure 7.
Posterior distributions for model fit to GW170817 when using the broad default priors (blue), compared to the posteriors obtained using themulti-messenger priors (black; same as Figure 5). Both sets of priors are listed in Table 1. turn depends on the binary source parameters. A better strategy infuture would be to tune the observed depths to a level necessary tocover the 90% credible range of kilonova luminosities for a given GWdetection . At present, this is not possible, as the GW-inferred binaryparameters needed to simulate the light curve are not released at thetime of merger. We suggest that greater cooperation between the GWand electromagnetic astronomy communities could be very useful inthis regard: if the GW source parameters ( M , 𝑞 , ˜ Λ ) were releasedin low latency, we could use our model to immediately predict the required observing depth at any wavelength to find or rule out akilonova counterpart at high significance. We have presented a simple forward model that combines a range ofliterature results on ejecta masses, r-process opacities, and the NS EoS,to predict kilonova light curves directly from GW-accessible binaryparameters. If such parameters (which may be estimated even in lowlatency; Biscoveanu et al. 2019; Finstad & Brown 2020; Krastev et al.
MNRAS , 1–15 (2020) ilonova light curves from binary parameters A pp a r e n t m a gn it ud e u r K0.0 0.5 1.0 1.5 2.0Time (days)161820222426 A pp a r e n t m a gn it ud e Using GW parameters u r K
Figure 8.
Model light curves for a GW190425 kilonova compared to limitscompiled by Hosseinzadeh et al. (2019). The shaded regions correspondto the 90% credible ranges in magnitude for each band. Top: assuming aGW170817-like kilonova shifted to the GW-inferred distance of GW190425(following Hosseinzadeh et al. 2019). Bottom: predicting the light curvedirectly from the GW constraints on the chirp mass, mass ratio, distance andinclination. The wider uncertainty bands and faster decline in this case showthat assuming a GW170817-like event leads to unrealistically tight constraintsand observations that may be too shallow to confidently detect this source. ∼ −
10 days after merger is a viscous disk wind. Thevery lanthanide-rich tidal ejecta makes a significant contribution onlyat very late times.The posteriors of our model fit contain important information onthe GW170817 progenitor system. We find that this system likelyhad a small mass ratio, 𝑞 ≈ .
9, and was viewed ≈ ◦ off-axis. Byincluding the maximum stable NS mass 𝑀 TOV and the symmetrictidal deformability Λ s in our free parameters, we constrained theseEoS-dependent parameters. We measure 𝑀 TOV = . + . − . M (cid:12) (90%credible interval) and derive 𝑅 . =11 . + . − . km using the posteriorsof Λ , M , 𝑞 and 𝑀 TOV , carefully taking into account the systematicerrors on all relations used. This radius constraint is consistent andcompetitive with others in the literature.We have not discussed in this work the signature of a mergerbetween a NS and a companion BH; these are also expected toproduce kilonovae as long as the mass ratio is modest (thoughwithout shock-driven, low-lanthanide ejecta). Our next step will be toimplement an analogous model to simulate light curves for this classof events.Finding a larger sample of kilonovae in the coming years, as GWdetectors and optical telescopes increase their survey power, will becrucial to improve our understanding of the NS EoS and the sourcepopulation of NS binaries. Hopefully, the code we present here willbe useful to the community both in helping to plan observations thatincrease this sample and in deriving important physics from the dataobtained.
ACKNOWLEDGEMENTS
MN thanks Jenni Barnes for providing kilonova thermalisation data,and Edo Berger, Laurence Datrier, Michael Fulton, Martin Hendry, Al-bert Kong, Gavin Lamb, Brian Metzger, Christopher Moore, GeraintPratten, Surojit Saha, Stephen Smartt, Alberto Vecchio and AshleyVillar for their insights and discussions. MN is supported by a RoyalAstronomical Society Research Fellowship. This project has receivedfunding from the European Research Council (ERC) under the Eu-ropean Union’s Horizon 2020 research and innovation programme(grant agreement No. 948381). BM is supported by NASA through theNASA Hubble Fellowship grant
DATA AVAILABILITY
This paper makes use of existing public data from kilonova.space .The code has been made available at https://github.com/guillochon/MOSFiT . REFERENCES
Abbott B. P., et al., 2017a, Physical Review Letters, 119, 161101Abbott B. P., et al., 2017b, ApJ, 848, L12 MNRAS000
Abbott B. P., et al., 2017a, Physical Review Letters, 119, 161101Abbott B. P., et al., 2017b, ApJ, 848, L12 MNRAS000 , 1–15 (2020) M. Nicholl et al
Abbott B. P., et al., 2018, Physical Review Letters, 121, 161101Abbott R., et al., 2020a, arXiv e-prints, p. arXiv:2010.14527Abbott B. P., et al., 2020b, Living Reviews in Relativity, 23, 3Abbott B. P., et al., 2020c, ApJ, 892, L3Abbott R., et al., 2020d, ApJ, 896, L44Ackley K., et al., 2020, A&A, 643, A113Alexander K. D., et al., 2017, ApJ, 848, L21Alexander K. D., et al., 2018, ApJ, 863, L18Anand S., et al., 2021, Nature Astronomy, 5, 46Andreoni I., et al., 2017, Publ. Astron. Soc. Australia, 34, e069Andreoni I., et al., 2020a, arXiv e-prints, p. arXiv:2008.00008Andreoni I., et al., 2020b, ApJ, 890, 131Antier S., et al., 2020, MNRAS, 492, 3904Antoniadis J., et al., 2013, Science, 340, 448Arcavi I., 2018, ApJ, 855, L23Arcavi I., et al., 2017, ApJ, 848, L33Arnett W. D., 1982, ApJ, 253, 785Barbieri C., Salafia O. S., Colpi M., Ghirlanda G., Perego A., 2020, arXive-prints, p. arXiv:2002.09395Barnes J., Kasen D., 2013, ApJ, 775, 18Barnes J., Kasen D., Wu M.-R., Martínez-Pinedo G., 2016, ApJ, 829, 110Bauswein A., Baumgarte T. W., Janka H. T., 2013a, Phys. Rev. Lett., 111,131101Bauswein A., Goriely S., Janka H.-T., 2013b, ApJ, 773, 78Bauswein A., Just O., Janka H.-T., Stergioulas N., 2017, ApJ, 850, L34Berger E., 2014, ARA&A, 52, 43Berger E., Fong W., Chornock R., 2013, ApJ, 774, L23Biscoveanu S., Vitale S., Haster C.-J., 2019, ApJ, 884, L32Blanchard P. K., et al., 2017, ApJ, 848, L22Breschi M., Perego A., Bernuzzi S., Del Pozzo W., Nedora V., Radice D.,Vescovi D., 2021, arXiv e-prints, p. arXiv:2101.01201Bulla M., 2019, MNRAS, 489, 5037Burbidge E. M., Burbidge G. R., Fowler W. A., Hoyle F., 1957, Reviews ofModern Physics, 29, 547Cantiello M., et al., 2018, ApJ, 854, L31Capano C. D., et al., 2020, Nature Astronomy, 4, 625Chen H.-Y., Cowperthwaite P. S., Metzger B. D., Berger E., 2020, arXive-prints, p. arXiv:2011.01211Chornock R., et al., 2017, ApJ, 848, L19Ciolfi R., Kalinani J. V., 2020, ApJ, 900, L35Ciolfi R., Kastaun W., Giacomazzo B., Endrizzi A., Siegel D. M., Perna R.,2017, Phys. Rev. D, 95, 063016Coughlin M. W., Dietrich T., Margalit B., Metzger B. D., 2019a, MNRAS,489, L91Coughlin M. W., et al., 2019b, ApJ, 885, L19Coulter D. A., et al., 2017, Science, 358, 1556Cowperthwaite P. S., et al., 2017, ApJ, 848, L17Cowperthwaite P. S., Chen H.-Y., Margalit B., Margutti R., May M., MetzgerB., Pankow C., 2019, arXiv e-prints, p. arXiv:1904.02718Cromartie H. T., et al., 2020, Nature Astronomy, 4, 72D’Avanzo P., et al., 2018, A&A, 613, L1Damour T., Nagar A., Villain L., 2012, Phys. Rev. D, 85, 123007Darbha S., Kasen D., 2020, ApJ, 897, 150Davies M. B., Benz W., Piran T., Thielemann F. K., 1994, ApJ, 431, 742De S., Siegel D., 2020, arXiv e-prints, p. arXiv:2011.07176De S., Finstad D., Lattimer J. M., Brown D. A., Berger E., Biwer C. M., 2018,Physical Review Letters, 121, 091102Demorest P. B., Pennucci T., Ransom S. M., Roberts M. S. E., Hessels J. W. T.,2010, Nature, 467, 1081Díaz M. C., et al., 2017, ApJ, 848, L29Dietrich T., Ujevic M., 2017, Classical and Quantum Gravity, 34, 105014Dietrich T., Coughlin M. W., Pang P. T. H., Bulla M., Heinzel J., Issa L., TewsI., Antier S., 2020, arXiv e-prints, p. arXiv:2002.11355Dobie D., et al., 2018, ApJ, 858, L15Drout M. R., et al., 2017, Science, 358, 1570Eichler D., Livio M., Piran T., Schramm D. N., 1989, Nature, 340, 126Evans P. A., et al., 2017, Science, 358, 1565Favata M., 2014, Phys. Rev. Lett., 112, 101101 Fernández R., Tchekhovskoy A., Quataert E., Foucart F., Kasen D., 2019,MNRAS, 482, 3373Finstad D., Brown D. A., 2020, ApJ, 905, L9Finstad D., De S., Brown D. A., Berger E., Biwer C. M., 2018, ApJ, 860, L2Flanagan É. É., Hinderer T., 2008, Phys. Rev. D, 77, 021502Fong W., et al., 2019, ApJ, 883, L1Fong W., et al., 2020, arXiv e-prints, p. arXiv:2008.08593Foucart F., O’Connor E., Roberts L., Kidder L. E., Pfeiffer H. P., Scheel M. A.,2016, Phys. Rev. D, 94, 123016Freiburghaus C., Rosswog S., Thielemann F. K., 1999, ApJ, 525, L121Gao H., Ai S.-K., Cao Z.-J., Zhang B., Zhu Z.-Y., Li A., Zhang N.-B., BausweinA., 2020, Frontiers of Physics, 15, 24603Ghirlanda G., et al., 2019, Science, 363, 968Gill R., Granot J., 2018, MNRAS, 478, 4128Goldstein A., et al., 2017, ApJ, 848, L14Gomez S., et al., 2019, ApJ, 884, L55Gompertz B. P., et al., 2018, ApJ, 860, 62Gompertz B. P., et al., 2020, MNRAS, 497, 726Goriely S., Bauswein A., Just O., Pllumbi E., Janka H. T., 2015, MNRAS,452, 3894Gottlieb O., Nakar E., Piran T., 2018, MNRAS, 473, 576Graham M. J., et al., 2020, Phys. Rev. Lett., 124, 251102Granot J., Gill R., Guetta D., De Colle F., 2018, MNRAS, 481, 1597Guillochon J., Parrent J., Kelley L. Z., Margutti R., 2017, ApJ, 835, 64Guillochon J., Nicholl M., Villar V. A., Mockler B., Narayan G., Mandel K. S.,Berger E., Williams P. K. G., 2018, ApJS, 236, 6Haggard D., Nynka M., Ruan J. J., Kalogera V., Cenko S. B., Evans P., KenneaJ. A., 2017, ApJ, 848, L25Hajela A., et al., 2019, ApJ, 886, L17Hinderer T., 2008, ApJ, 677, 1216Hosseinzadeh G., et al., 2019, ApJ, 880, L4Hotokezaka K., Kiuchi K., Kyutoku K., Okawa H., Sekiguchi Y.-i., ShibataM., Taniguchi K., 2013, Phys. Rev. D, 87, 024001Hotokezaka K., Nakar E., Gottlieb O., Nissanke S., Masuda K., Hallinan G.,Mooley K. P., Deller A. T., 2019, Nature Astronomy, 3, 940Hu L., et al., 2017, Science Bulletin, Vol.~62, No.21, p.1433-1438, 2017, 62,1433Jin Z.-P., et al., 2016, Nature Communications, 7, 12898Kasen D., Badnell N. R., Barnes J., 2013, ApJ, 774, 25Kasen D., Fernández R., Metzger B. D., 2015, MNRAS, 450, 1777Kasen D., Metzger B., Barnes J., Quataert E., Ramirez-Ruiz E., 2017, Nature,551, 80Kasliwal M. M., et al., 2017, Science, 358, 1559Kasliwal M. M., et al., 2020, arXiv e-prints, p. arXiv:2006.11306Kawaguchi K., Shibata M., Tanaka M., 2018, ApJ, 865, L21Kiuchi K., Kyutoku K., Shibata M., Taniguchi K., 2019, ApJ, 876, L31Klion H., Duffell P. C., Kasen D., Quataert E., 2021, MNRAS,Korobkin O., Rosswog S., Arcones A., Winteler C., 2012, MNRAS, 426, 1940Korobkin O., et al., 2020, arXiv e-prints, p. arXiv:2004.00102Krastev P. G., Gill K., Villar V. A., Berger E., 2020, arXiv e-prints, p.arXiv:2012.13101Krüger C. J., Foucart F., 2020, Phys. Rev. D, 101, 103002Kulkarni S. R., 2005, arXiv e-prints, pp astro–ph/0510256LSST Science Collaboration et al., 2009, arXiv e-prints, p. arXiv:0912.0201Lamb G. P., et al., 2019, ApJ, 883, 48Lattimer J. M., Prakash M., 2016, Phys. Rep., 621, 127Lattimer J. M., Schramm D. N., 1974, ApJ, 192, L145Lazzati D., Perna R., Morsony B. J., Lopez-Camara D., Cantiello M., CiolfiR., Giacomazzo B., Workman J. C., 2018, Phys. Rev. Lett., 120, 241103Lehner L., Liebling S. L., Palenzuela C., Caballero O. L., O’Connor E.,Anderson M., Neilsen D., 2016, Classical and Quantum Gravity, 33,184002Levan A. J., et al., 2017, ApJ, 848, L28Li L.-X., Paczyński B., 1998, ApJ, 507, L59Lippuner J., Fernández R., Roberts L. F., Foucart F., Kasen D., Metzger B. D.,Ott C. D., 2017, MNRAS, 472, 904Lipunov V. M., et al., 2017, ApJ, 850, L1Lundquist M. J., et al., 2019, ApJ, 881, L26MNRAS , 1–15 (2020) ilonova light curves from binary parameters Lyman J. D., et al., 2018, Nature Astronomy, 2, 751Margalit B., Metzger B. D., 2017, ApJ, 850, L19Margalit B., Metzger B. D., 2019, ApJ, 880, L15Margutti R., Chornock R., 2020, arXiv e-prints, p. arXiv:2012.04810Margutti R., et al., 2017, ApJ, 848, L20Margutti R., et al., 2018a, arXiv e-prints, p. arXiv:1812.04051Margutti R., et al., 2018b, ApJ, 856, L18McBrien O. R., et al., 2020, MNRAS,McCully C., et al., 2017, ApJ, 848, L32Metzger B. D., 2019, Living Reviews in Relativity, 23, 1Metzger B. D., Fernández R., 2014, MNRAS, 441, 3444Metzger B. D., et al., 2010, MNRAS, 406, 2650Metzger B. D., Bauswein A., Goriely S., Kasen D., 2015, MNRAS, 446, 1115Metzger B. D., Thompson T. A., Quataert E., 2018, ApJ, 856, 101Mooley K. P., et al., 2018, Nature, 554, 207Nakar E., Piran T., 2017, ApJ, 834, 28Nakar E., Piran T., 2020, arXiv e-prints, p. arXiv:2005.01754Narayan R., Paczynski B., Piran T., 1992, ApJ, 395, L83Nedora V., et al., 2020, arXiv e-prints, p. arXiv:2011.11110Nicholl M., et al., 2017a, ApJ, 848, L18Nicholl M., Guillochon J., Berger E., 2017b, ApJ, 850, 55O’Connor B., et al., 2021, MNRAS,Pan Y.-C., et al., 2017, ApJ, 848, L30Pian E., et al., 2017, Nature, 551, 67Piro A. L., Kollmeier J. A., 2018, ApJ, 855, 103Postnikov S., Prakash M., Lattimer J. M., 2010, Phys. Rev. D, 82, 024016Pozanenko A. S., et al., 2018, ApJ, 852, L30Radice D., Perego A., Hotokezaka K., Fromm S. A., Bernuzzi S., RobertsL. F., 2018, ApJ, 869, 130Raithel C. A., Özel F., Psaltis D., 2018, ApJ, 857, L23Raithel C., Ozel F., Psaltis D., 2020, arXiv e-prints, p. arXiv:2004.00656Rastinejad J. C., et al., 2021, arXiv e-prints, p. arXiv:2101.03175Rosswog S., Liebendörfer M., Thielemann F.-K., Davies M. B., Benz W.,Piran T., 1999, A&A, 341, 499Ryan G., van Eerten H., Piro L., Troja E., 2020, ApJ, 896, 166Salafia O. S., Colpi M., Branchesi M., Chassande-Mottin E., Ghirlanda G.,Ghisellini G., Vergani S. D., 2017, ApJ, 846, 62Savchenko V., et al., 2017, ApJ, 848, L15Sekiguchi Y., Kiuchi K., Kyutoku K., Shibata M., 2015, Phys. Rev. D, 91,064059Sekiguchi Y., Kiuchi K., Kyutoku K., Shibata M., Taniguchi K., 2016, Phys.Rev. D, 93, 124046Shapiro S. L., Teukolsky S. A., 1986, Black Holes, White Dwarfs and NeutronStars: The Physics of Compact ObjectsShappee B. J., et al., 2017, Science, 358, 1574Shibata M., Hotokezaka K., 2019, Annual Review of Nuclear and ParticleScience, 69, 41Shibata M., Zhou E., Kiuchi K., Fujibayashi S., 2019, Phys. Rev. D, 100,023015Siegel D. M., Metzger B. D., 2017, Phys. Rev. Lett., 119, 231102Siegel D. M., Metzger B. D., 2018, ApJ, 858, 52Smartt S. J., et al., 2017, Nature, 551, 75Smith G. P., Robertson A., Bianconi M., Jauzac M., 2019, arXiv e-prints, p.arXiv:1902.05140Soares-Santos M., et al., 2017, ApJ, 848, L16Speagle J. S., 2020, MNRAS, 493, 3132Tanaka M., Kato D., Gaigalas G., Kawaguchi K., 2020, MNRAS, 496, 1369Tanvir N. R., Levan A. J., Fruchter A. S., Hjorth J., Hounsell R. A., WiersemaK., Tunnicliffe R. L., 2013, Nature, 500, 547Tanvir N. R., et al., 2017, ApJ, 848, L27Tews I., Margueron J., Reddy S., 2018, Phys. Rev. C, 98, 045804Troja E., et al., 2017, Nature, 551, 71Troja E., et al., 2018, Nature Communications, 9, 4089Troja E., et al., 2019, MNRAS, 489, 2104Utsumi Y., et al., 2017, PASJ, 69, 101Valenti S., et al., 2017, ApJ, 848, L24Vieira N., et al., 2020, ApJ, 895, 96Villar V. A., et al., 2017, ApJ, 851, L21 Wade L., Creighton J. D. E., Ochsner E., Lackey B. D., Farr B. F., LittenbergT. B., Raymond V., 2014, Phys. Rev. D, 89, 103012Wanajo S., Sekiguchi Y., Nishimura N., Kiuchi K., Kyutoku K., Shibata M.,2014, ApJ, 789, L39Wollaeger R. T., et al., 2018, MNRAS, 478, 3298Wu Y., MacFadyen A., 2019, ApJ, 880, L23Yagi K., Yunes N., 2016, Classical and Quantum Gravity, 33, 13LT01Yagi K., Yunes N., 2017, Phys. Rep., 681, 1
APPENDIX A: NEW FITS TO SIMULATION DATA
Figure A1 shows our quadratic fit to determine 𝑓 red and 𝑓 blue ≡ − 𝑓 red ,the fractions of red (high opacity) and blue (low opacity) dynamicalejecta, using the simulated 𝑌 𝑒 distributions from Sekiguchi et al.(2016). Ejecta with 𝑌 𝑒 ≤ .
25 are assumed to be red.Figure A2 shows our broken polynomial fit to the simulations ofTanaka et al. (2020), which we use to estimate the opacity of thepurple (post-merger) ejecta as a function of 𝑌 𝑒 . APPENDIX B: PLOTS OF ADDITIONAL LIGHT CURVEFITS
Figure B1 shows the light curve fits for the surface-wind-only andshock-only versions of the model.
This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000
This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000 , 1–15 (2020) M. Nicholl et al /M M a ss fr ac ti on w it h Y e < . SFHoDD2
Figure A1.
Mass fraction of lanthanide-poor dynamical ejecta as a function of the binary mass ratio, 𝑞 = 𝑀 / 𝑀 , using 𝑌 𝑒 distributions from Sekiguchi et al.(2016) for the SFHo (soft) and DD2 (stiff) equations of state. Our polynomial fit predicts essentially no lanthanide-poor ejecta for 𝑞 (cid:46) .
8, in agreement withDietrich & Ujevic (2017). e O p ac it y ( c m g - ) Figure A2.
Opacity of r-process ejecta as a function of electron fraction and temperature from Tanaka et al. (2020). The broken polynomial fit is used to estimatethe opacity of the disk ejecta, using the average 𝑌 𝑒 predicted by Lippuner et al. (2017) for stable, short-lived or prompt-collapse merger remnants. In practice, theaverage 𝑌 𝑒 > .
25 for all reasonable input parameters.MNRAS , 1–15 (2020) ilonova light curves from binary parameters -1 Days since merger12151821242730 A pp a r e n t m a gn it ud e + c on s t a n t W2+7M2+5W1+5 U+5u+4B+3 g+2V+1r+0 i-1z-2y-3 J-4H-5K-6 -1 Days since merger12151821242730 A pp a r e n t m a gn it ud e + c on s t a n t W2+7M2+5W1+5 U+5u+4B+3 g+2V+1r+0 i-1z-2y-3 J-4H-5K-6
Figure B1.
Left: Best-fit kilonova model for GW170817 with surface ejecta enhancement. Right: Best-fit model with shock cooling. The shock cooling model ispreferred with a Bayes factor
𝐵 > . This figure is similar to Figure 4 in the main text, which shows the best fit light curves including either both or neither ofthese effects. MNRAS000