Multi-messenger Bayesian parameter inference of a binary neutron-star merger
The combined detection of a binary neutron-star merger in both gravitational waves (GWs) and electromagnetic (EM) radiation spanning the entire spectrum – GW170817 / AT2017gfo / GRB170817A – marks a breakthrough in the field of multi-messenger astronomy. Between the plethora of modeling and observations, the rich synergy that exists among the available data sets creates a unique opportunity to constrain the binary parameters, the equation of state of supranuclear density matter, and the physical processes at work during the kilonova and gamma-ray burst. However, previous works use simplified lightcurve models and fits to numerical relativity simulation that do not account for all of the relevant physical processes. We report, for the first time, Bayesian parameter estimation combining information from GW170817, AT2017gfo, GRB170817 to obtain truly multi-messenger constraints on the tidal deformability , total binary mass , the radius of a solar mass neutron star (with additional systematic uncertainty), and an upper bound on the mass ratio of , all at 90% confidence. Our joint novel analysis makes use of new phenomenological descriptions of the dynamical ejecta, debris disk mass, and remnant black hole properties, all derived from a large suite of numerical relativity simulations.
The combined detection of a GW event, GW170817 TheLIGOScientific:2017qsa (), a gamma ray burst (GRB) of short duration, GRB170817A Monitor:2017mdv () accompanied by a non-thermal afterglow, and thermal emission (“kilonova”) at optical, near-infrared, and ultraviolet wavelengths, AT2017gfo GBM:2017lvd (); Arcavi:2017xiz (); Coulter:2017wya (); Lipunov:2017dwd (); Soares-Santos:2017lru (); Tanvir:2017pws (); Valenti:2017ngx () from a binary neutron star (BNS) merger has enabled major leaps forward in several research areas. The latter include new limits on the equation of state (EOS) of cold matter at supranuclear densities (e.g. De:2018uhw (); Abbott:2018exr (); Radice:2017lry (); Radice:2018ozg (); Coughlin:2018miv (); Bauswein:2017vtn (); Annala:2017llu (); Most:2018hfd (); Ruiz:2017due (); Margalit:2017dij (); Rezzolla:2017aly (); Shibata:2017xdx ()). One of the main goals of the nascent field of “multi-messenger astronomy” is to obtain complementary observations of the same object or event. These observations, potentially across a variety of wavelengths and particle types, probe different aspects of the system. In the case of GW170817, GW detectors such as LIGO and Virgo provide a highly accurate measurement of the binary chirp mass , but leave the mass ratio, , poorly constrained.
A variety of studies over the last year focused on the properties of this first detection of a BNS system, including detailed analyses of the GW signal by the LVC (e.g. TheLIGOScientific:2017qsa (); Abbott:2018wiz (); Abbott:2018exr (); LIGOScientific:2018mvr ()) and external groups (e.g. De:2018uhw (); Finstad:2018wid (); Dai:2018dca ()), relying on different parameter estimation techniques and a variety of GW models. Despite this diversity of methods, all of the published works predict small tidal deformabilities, favoring relatively soft EOSs and placing upper limits on the radii of NSs. For this first BNS system the GW analyses broadly agree, and studies indicate that systematic errors are below the statistical errors Abbott:2018wiz (); Dudi:2018jzn (); Samajdar:2018dcx (). However, this might not be the case for future GW observations with larger signal-to-noise ratios, thus emphasizing the need for further improvements in the current infrastructure and GW modeling.
Fortunately, deficiencies in the available GW information can sometimes be supplemented with EM observations, potentially improving the measurements of key parameters. For instance, the results of numerical relativity simulations were used to argue against the EOS being too soft, as the mass of the remnant accretion disk and its associated wind ejecta would be insufficient to account for the luminosity of the observed kilonova, e.g., Radice:2018xqa (); Bauswein:2017vtn (); Coughlin:2018miv (). Combining GW and EM observations thus provides an opportunity to independently constrain the binary parameters, place tighter bounds on the EOS, and obtain a better understanding of the physical processes and outcomes of BNS mergers.
One of the first multi-messenger constraints on the tidal deformability and consequently the supranuclear EOS was presented in Radice:2017lry (). Based on state-of-the art numerical relativity (NR) simulations, they proposed that the tidal deformability needs to be to ensure that a significant fraction of matter was either ejected from the system or contained within a debris disk around the BH remnant to explain the bright EM counterpart. Recently, Radice:2018ozg () updated this first analysis and obtained constraints on the tidal deformability of and on the corresponding radius of a 1.4 neutron star of km. To the best of our knowledge, Coughlin:2018miv () presented the first analysis of the lightcurves and spectra of AT2017gfo linking with a Bayesian analysis the kilonova properties to the source properties of the binary. We used the kilonova model of Kasen:2017sxr () combined with methods of Gaussian Process Regression (GPR), Doctor:2017csx (); Purrer:2014fza (); Coughlin:2018miv (), and related a fraction of the ejected material to dynamical ejecta. Based on the analysis, the tidal deformability was limited to .
In addition, there have been studies placing limits on the maximum NS mass of a stable TOV star, . Those studies are orthogonal to the works constraining the tidal deformability since both quantities test different parts of the NS EOS. Ref. Margalit:2017dij () places a upper limit on the mass of a non-rotating NS of , Rezzolla:2017aly () report a maximum TOV mass of , and Shibata:2017xdx () provide an estimate for the maximum mass of . All these constraints have been derived by assuming the formation of a BH after the merger of GW170817 and incorporating the measured chirp mass inferred from the GW analysis. We employ for our analysis the maximum mass constraint derived in Margalit:2017dij ().
In this article, we will perform an analysis combining information from three separate sources: GW170817, GRB170817A, and the kilonova AT2017gfo to perform a multi-messenger Bayesian parameter analysis of a neutron star merger. The flowchart in Fig. 1 highlights the interplay between the different observable signatures and presents the joint posteriors obtained on the tidal deformability , the binary mass ratio , and the maximum mass of a stable non-rotating neutron star .
We begin by analyzing GW170817 (blue shaded region of Fig. 1) and use the publicly available “low spin” posterior samples (https://dcc.ligo.org/LIGO-P1800061) published by the LIGO Scientific and Virgo Collaborations (LVC) TheLIGOScientific:2017qsa (); Abbott:2018wiz (). As these sample use the sky localization obtained from EM observations, they already incorporate EM information. Under the assumption that the merging objects are two NSs described by the same EOS De:2018uhw (); Abbott:2018exr (), we can further restrict the posterior distribution. Finally, we discard from our sample those systems with viewing angles which are inconsistent with those allowed by the analysis of the GRB afterglow by vanEerten:2018vgj (). These restrictions more tightly constrain the tidal deformability Abbott:2018exr (), in particular showing that the GW data disfavors values of , motivating us to conservatively restrict in the following EM analysis.
In the second phase of our work, we analyze the light curves of AT2017gfo (red shaded region in Fig. 1). We fit the observational data Coughlin:2018miv (); Smartt:2017fuw (); GBM:2017lvd () with the 2-component radiative transfer model of Kasen:2017sxr (). The usage of multiple components, proposed prior to the discovery of GW170817 Metzger:2014ila (), is motivated by different ejecta mechanisms contributing to the total -process yields of BNS mergers. The first type of mass ejection are “dynamical ejecta” generated during the merger process itself. Dynamical ejecta are typically characterized by a low-electron fraction when they are created by tidal torque, but the electron fraction can extend to higher values (and thus the lanthanide abundance be reduced) in the case of shock-driven ejecta. In addition to dynamical ejecta, disk winds driven by neutrino energy, magnetic fields, viscous evolution and/or nuclear recombination (e.g. Kohri:2005tq (); Surman:2005kf (); Metzger:2008av (); Dessart:2008zd (); Fernandez:2013tya (); Perego:2014fma (); Siegel:2014ita (); Just:2014fka (); Rezzolla:2014nva (); Ciolfi:2014yla (); Siegel:2017nub ()) leads to a large quantity of ejecta, which in many cases exceeds that of the dynamical component. The ejecta components employed in our kilonova light curve analysis are related to these different physical ejecta mechanisms: the first ejecta component is assumed to be proportional to dynamical ejecta, , while the second ejecta component arises from the disk wind and is assumed to be proportional to the mass of the remnant disk, . By fitting the observed lightcurves with the kilonova models Kasen:2017sxr () withing a Gaussian Process Regression framework Coughlin:2018miv (), we obtain for each component posterior distributions for the ejecta mass , the lanthanide mass fraction (related to the initial electron fraction), and the ejecta velocity .
The values of , , and obtained from our kilonova analysis are related to the properties of the binary and EOS using new phenomenological fits to numerical relativity simulations, which we briefly described below. First, we revisit the phenomenological fit presented in Radice:2018pdn () between the disk mass and tidal deformability to correlate the disk mass, , to the properties of the merging binary. Simulations following the merger aftermath suggest that the disk mass is accumulated primarily through radial redistribution of matter in the post-merger remnant. Thus, the lifetime of the remnant prior to its collapse is related to its stability and found to strongly correlate with the disk mass Radice:2018xqa (). We find that the lifetime in turn is governed to a large degree by the ratio of , where is the total binary mass and is the threshold mass Bauswein:2013jpa () above which the merger results in prompt (dynamical timescale) collapse to a black hole, which depends on the NS compactness and thus . Therefore, , rather than alone, provides a better measure of the stability of the post-merger remnant, and following the arguments above, is expected to correlate with .
Fig. 2 shows, based on the suite of numerical relativity simulations of Radice:2018pdn (), that there indeed exists a relatively tight correlation between the accretion disk mass and . For our analysis, we will use
Connecting the NS radius to the chirp mass and tidal deformability, De:2018uhw (), we conclude that the disk mass (and thus the disk wind ejecta) is a function of the tidal deformability, total binary mass, and the maximum TOV mass, . Therefore, both information, those on the densest portion of the EOS, which controls , and those from lower densities, as encoded in or , play a role in controlling the disk mass and kilonova properties. The inclusion of these parameters and slight changes in the functional form of the phenomenological relation decrease the average fractional errors by more than a factor of 3 relative to previous disk mass estimates based on alone Radice:2018xqa (), thus reducing uncertainties and errors on the EOS constraints obtained from kilonova observations.
Another key ingredient in our analysis is the role of the dynamical ejecta as the first kilonova ejecta component. Based on a suite of numerical relativity simulations obtained by different groups and codes, Ref. Dietrich:2016fpt () derived the first phenomenological fit for the dynamical ejecta for BNS systems. This fit (in its original or upgraded version) has been employed in a number of studies, including the analysis of GW170817 Abbott:2017wuw (); Coughlin:2018miv (), and they have been updated in Coughlin:2018miv () and Radice:2018pdn (). Here, we present a further upgrade which incorporates the new numerical relativity dataset of Radice:2018pdn () and uses the fitting function of Coughlin:2018miv () (which fits instead of ). The extend dataset contains a total of 259 numerical relativity simulations. The final fitting function is
with , , , and and denoting the compactnesses of the individual stars, a more detailed discussed can be found in the Appendix.
A final ingredient in relating observational data to the binary parameters are phenomenological fits for the BH mass and spin. One finds that with an increasing total mass , the final black hole mass and angular momentum increases almost linearly. For unequal mass mergers, and decrease with . Considering the imprint of the EOS, we find that for larger values of , the final black hole mass decreases, which follows from the observation that the disk mass increases with . We finally obtain:
with and and
with , , and ; further details are given in the Appendix.
In addition to using these fits, we use the results of Margalit:2017dij (), who derive a upper limit on the mass of a non-rotating NS of based on energetic considerations from the GRB and kilonova which rule out a long-lived supramassive NS remnant, to place a prior on between 2–2.17 . Combining these phenomenological relations with the lightcurve data, our analysis strongly favors equal or nearly equal mass systems and (see appendix). We conclude that roughly of the first ejecta component is associated with dynamical ejecta, while about of the disk mass must be ejected in winds to account for the second ejecta component. The latter agrees with the results of long-term general relativistic magnetohydrodynamical simulations of the post-merger accretion disk (e.g. Siegel:2017nub ()). If we do not enforce constraints on , we obtain similar constraints in the binary parameters but with allowed values . This is broadly consistent with the results presented in Margalit:2017dij (); Rezzolla:2017aly (); Ruiz:2017due () and provides a new and independent measurement of the maximum TOV mass, which will become more accurate with future multi-messenger events.
Our third and final result uses Bayesian parameter estimation of GRB170817A directly (green shaded region in Fig. 1). We assume that the GRB jet is powered by the accretion of matter from the debris disk onto the BH Eichler:1989ve (); Paczynski:1991aq (); Meszaros:1992ps (); Narayan:1992iy () and that the jet energy is proportional to the disk mass. Accounting for the loss of disk mass to winds, we connect our estimates of the disk wind ejecta from the analysis of AT2017gfo to the following GRB parameter estimation analysis. In order to assess the robustness of our conclusions, and to evaluate potential systematic uncertainties, we show results for three different fits to the GRB afterglow: van Eerten et al vanEerten:2018vgj (), Wu and MacFadyen Wu:2018bxg (), and Wang et al Wang:2018nye (). While the analyses of Refs. vanEerten:2018vgj (); Wu:2018bxg () differ on the energy of the GRB, the use of either one further constrains the value of and the binary mass ratio, shifting both to slightly higher values than obtained through the analysis of AT2017gfo alone.
To obtain the final constraints on the EOS and binary properties, we combine the posteriors obtained from GW170817 and AT2017gfo+GRB170817A. The analysis of AT2017gfo and GRB170817A are highly correlated, as both use the same phenomenological description for the disk mass and the AT2017gfo posteriors are employed as priors for the GRB analysis. However, we assume the parameter estimations results from the GW and EM analysis for and are independent from one another. Thus, the final multi-messenger probability density function is given by:
We summarize our constraints on the binary parameters and EOS in Table 1. The final constraints on the tidal deformability and the mass ratio are shown at the bottom of Figure 1, where we use the GRB model of Ref. vanEerten:2018vgj () (similar constraints are obtained with the other GRB models). According to our analysis, the confidence interval for the tidal deformability is . The distribution has its percentile at . Relating the measured confidence interval to the NS radius De:2018uhw (), we obtain a constraint on the NS radius of (with a uncertainty of the quasi-universal relation De:2018uhw (); Radice:2018ozg () connecting and ). This result is in good agreement with that recently obtained by the multi-messenger analysis presented in Radice and Dai Radice:2018ozg (). Considering the constraint on the mass ratio, we find that at confidence. Combining this with the measured chirp mass, the total binary mass lies in the range . The radius constraint, together with the constraint on the maximum TOV-mass, can be used to rule out or favor a number of proposed NS EOSs, as illustrated in Fig. 3.
While many analyses of GW170817 and its electromagnetic signatures have been presented in the literature, that presented here is the first to combine information from all three channels: GW170817, GRB170817A, and AT2017gfo. While this article was in preparation, Ref. Radice:2018ozg () extended the work of Radice:2017lry () in performing a multi-messenger parameter estimation based on relative binning incorporating information from the disk mass derived in Radice:2018pdn (). Our work makes use of more available knowledge than employed in any previous multi-messenger analyses. In particular, our final posteriors describe the observed lightcurve data of AT2017gfo, which have been improved significantly compared to previous works, and explain the properties of GRB170817A. Future improvements on this analysis will require a larger suite of light curve models, which account for additional details such as the potentially aspherical geometry of the merger ejecta. Radiative transfer simulations which account for these asymmetries will enable additional constraints on the system viewing angle. Future simulation work should also explore the inevitable interplay between the ejecta components, which could complicate the simple dynamical ejecta/disk wind ejecta dichotomy adopted in our analysis.
Acknowledgements.We that Zoheyr Doctor, Michael Pürrer, and David Radice for helpful discussions and comments on the manuscript. MWC is supported by the David and Ellen Lee Postdoctoral Fellowship at the California Institute of Technology. TD acknowledges support by the European Union’s Horizon 2020 research and innovation program under grant agreement No 749145, BNSmergers. BDM acknowledges support from NASA (grant number NNX16AB30G).
- (1) Abbott B P et al. (Virgo, LIGO Scientific) 2017 Phys. Rev. Lett. 119 161101 (Preprint eprint 1710.05832)
- (2) Abbott B P et al. (Virgo, Fermi-GBM, INTEGRAL, LIGO Scientific) 2017 Astrophys. J. 848 L13 (Preprint eprint 1710.05834)
- (3) Abbott B P et al. (GROND, SALT Group, OzGrav, DFN, DES, INTEGRAL, Virgo, Insight-Hxmt, MAXI Team, Fermi-LAT, J-GEM, RATIR, IceCube, CAASTRO, LWA, ePESSTO, GRAWITA, RIMAS, SKA South Africa/MeerKAT, H.E.S.S., 1M2H Team, IKI-GW Follow-up, Fermi GBM, Pi of Sky, DWF (Deeper Wider Faster Program), MASTER, AstroSat Cadmium Zinc Telluride Imager Team, Swift, Pierre Auger, ASKAP, VINROUGE, JAGWAR, Chandra Team at McGill University, TTU-NRAO, GROWTH, AGILE Team, MWA, ATCA, AST3, TOROS, Pan-STARRS, NuSTAR, ATLAS Telescopes, BOOTES, CaltechNRAO, LIGO Scientific, High Time Resolution Universe Survey, Nordic Optical Telescope, Las Cumbres Observatory Group, TZAC Consortium, LOFAR, IPN, DLT40, Texas Tech University, HAWC, ANTARES, KU, Dark Energy Camera GW-EM, CALET, Euro VLBI Team, ALMA) 2017 Astrophys. J. 848 L12 (Preprint eprint 1710.05833)
- (4) Arcavi I et al. 2017 Nature 551 64 (Preprint eprint 1710.05843)
- (5) Coulter D A et al. 2017 Science [Science358,1556(2017)] (Preprint eprint 1710.05452)
- (6) Lipunov V M et al. 2017 Astrophys. J. 850 L1 (Preprint eprint 1710.05461)
- (7) Soares-Santos M et al. (DES) 2017 Astrophys. J. 848 L16 (Preprint eprint 1710.05459)
- (8) Tanvir N R et al. 2017 Astrophys. J. 848 L27 (Preprint eprint 1710.05455)
- (9) Valenti S, Sand D J, Yang S, Cappellaro E, Tartaglia L, Corsi A, Jha S W, Reichart D E, Haislip J and Kouprianov V 2017 Astrophys. J. 848 L24 (Preprint eprint 1710.05854)
- (10) De S, Finstad D, Lattimer J M, Brown D A, Berger E and Biwer C M 2018 Phys. Rev. Lett. 121 091102 (Preprint eprint 1804.08583)
- (11) Abbott B P et al. (Virgo, LIGO Scientific) 2018 arXiv: 1805.11581 (Preprint eprint 1805.11581)
- (12) Radice D, Perego A, Zappa F and Bernuzzi S 2018 Astrophys. J. 852 L29 (Preprint eprint 1711.03647)
- (13) Radice D and Dai L 2018 arXiv:1810.12917 (Preprint eprint 1810.12917)
- (14) Coughlin M W, Dietrich T, Doctor Z, Kasen D, Coughlin S, Jerkstrand A, Leloudas G, McBrien O, Metzger B D, OâShaughnessy R and Smartt S J 2018 Monthly Notices of the Royal Astronomical Society 480 3871–3878 (Preprint eprint 1805.09371)
- (15) Bauswein A, Just O, Janka H T and Stergioulas N 2017 Astrophys. J. 850 L34 (Preprint eprint 1710.06843)
- (16) Annala E, Gorda T, Kurkela A and Vuorinen A 2018 Phys. Rev. Lett. 120 172703 (Preprint eprint 1711.02644)
- (17) Most E R, Weih L R, Rezzolla L and Schaffner-Bielich J 2018 Phys. Rev. Lett. 120 261103 (Preprint eprint 1803.00549)
- (18) Ruiz M, Shapiro S L and Tsokaros A 2018 Phys. Rev. D97 021501 (Preprint eprint 1711.00473)
- (19) Margalit B and Metzger B D 2017 Astrophys. J. 850 L19 (Preprint eprint 1710.05938)
- (20) Rezzolla L, Most E R and Weih L R 2018 Astrophys. J. 852 L25 (Preprint eprint 1711.00314)
- (21) Shibata M, Fujibayashi S, Hotokezaka K, Kiuchi K, Kyutoku K, Sekiguchi Y and Tanaka M 2017 Phys. Rev. D96 123012 (Preprint eprint 1710.07579)
- (22) Abbott B P et al. (Virgo, LIGO Scientific) 2018 arXiv: 1805.11579 (Preprint eprint 1805.11579)
- (23) Abbott B P et al. (LIGO Scientific, Virgo) 2018 (Preprint eprint 1811.12907)
- (24) Finstad D, De S, Brown D A, Berger E and Biwer C M 2018 Astrophys. J. 860 L2 (Preprint eprint 1804.04179)
- (25) Dai L, Venumadhav T and Zackay B 2018 arXiv: 1806.08793
- (26) Dudi R, Pannarale F, Dietrich T, Hannam M, Bernuzzi S, Ohme F and Brügmann B 2018 Phys. Rev. D98 084061 (Preprint eprint 1808.09749)
- (27) Samajdar A and Dietrich T 2018 arXiv: 1810.03936
- (28) Radice D, Perego A, Bernuzzi S and Zhang B 2018 arXiv:1803.10865 (Preprint eprint 1803.10865)
- (29) Kasen D, Metzger B, Barnes J, Quataert E and Ramirez-Ruiz E 2017 Nature, 10.1038/nature24453 (Preprint eprint 1710.05463)
- (30) Doctor Z, Farr B, Holz D E and PÃ¼rrer M 2017 Phys. Rev. D96 123011 (Preprint eprint 1706.05408)
- (31) Pürrer M 2014 Class. Quant. Grav. 31 195010 (Preprint eprint 1402.4146)
- (32) van Eerten E T H, Ryan G, Ricci R, Burgess J M, Wieringa M, Piro L, Cenko S B and Sakamoto T 2018 arXiv:1808.06617 (Preprint eprint 1808.06617)
- (33) Smartt S J et al. 2017 Nature, 10.1038/nature24303 (Preprint eprint 1710.05841)
- (34) Metzger B D and FernÃ¡ndez R 2014 Mon.Not.Roy.Astron.Soc. 441 3444 (Preprint eprint 1402.4803)
- (35) Kohri K, Narayan R and Piran T 2005 Astrophys. J. 629 341–361 (Preprint eprint astro-ph/0502470)
- (36) Surman R, McLaughlin G C and Hix W R 2006 Astrophys. J. 643 1057–1064 (Preprint eprint astro-ph/0509365)
- (37) Metzger B, Piro A and Quataert E 2008 Mon.Not.Roy.Astron.Soc. 390 781 (Preprint eprint 0805.4415)
- (38) Dessart L, Ott C, Burrows A, Rosswog S and Livne E 2009 Astrophys.J. 690 1681 (Preprint eprint 0806.4380)
- (39) FernÃ¡ndez R and Metzger B D 2013 Mon. Not. Roy. Astron. Soc. 435 502 (Preprint eprint 1304.6720)
- (40) Perego A, Rosswog S, Cabezon R, Korobkin O, Kaeppeli R et al. 2014 Mon.Not.Roy.Astron.Soc. 443 3134 (Preprint eprint 1405.6730)
- (41) Siegel D M, Ciolfi R and Rezzolla L 2014 Astrophys. J. 785 L6 (Preprint eprint 1401.4544)
- (42) Just O, Bauswein A, Pulpillo R A, Goriely S and Janka H T 2015 Mon. Not. Roy. Astron. Soc. 448 541–567 (Preprint eprint 1406.2687)
- (43) Rezzolla L and Kumar P 2015 Astrophys. J. 802 95 (Preprint eprint 1410.8560)
- (44) Ciolfi R and Siegel D M 2015 Astrophys.J. 798 L36 (Preprint eprint 1411.2015)
- (45) Siegel D M and Metzger B D 2017 Phys. Rev. Lett. 119 231102 (Preprint eprint 1705.05473)
- (46) Radice D, Perego A, Hotokezaka K, Fromm S A, Bernuzzi S and Roberts L F 2018 arXiv: 1809.11161
- (47) Bauswein A, Baumgarte T and Janka H T 2013 Phys.Rev.Lett. 111 131101 (Preprint eprint 1307.5191)
- (48) Dietrich T and Ujevic M 2017 Class. Quant. Grav. 34 105014 (Preprint eprint 1612.03665)
- (49) Abbott B P et al. (Virgo, LIGO Scientific) 2017 Astrophys. J. 850 L39 (Preprint eprint 1710.05836)
- (50) Eichler D, Livio M, Piran T and Schramm D N 1989 Nature 340 126–128
- (51) Paczynski B 1991 Acta Astron. 41 257–267
- (52) Meszaros P and Rees M J 1992 Astrophys. J. 397 570–575
- (53) Narayan R, Paczynski B and Piran T 1992 Astrophys. J. 395 L83–L86 (Preprint eprint astro-ph/9204001)
- (54) Wu Y and MacFadyen A 2018 arXiv:1809.06843 (Preprint eprint 1809.06843)
- (55) Wang Y Z, Shao D S, Jiang J L, Tang S P, Ren X X, Jin Z P, Fan Y Z and Wei D M 2018 arXiv:1811.02558 (Preprint eprint 1811.02558)
- (56) Stovall K et al. 2018 Astrophys. J. 854 L22 (Preprint eprint 1802.01707)
- (57) Dietrich T et al. 2018 arXiv: 1804.02235 (Preprint eprint 1804.02235)
- (58) Yagi K and Yunes N 2016 Class. Quant. Grav. 33 13LT01 (Preprint eprint 1512.02639)
- (59) Chatziioannou K, Haster C J and Zimmerman A 2018 Phys. Rev. D97 104036 (Preprint eprint 1804.03221)
- (60) Mooley K P et al. 2018 Nature 554 207 (Preprint eprint 1711.11573)
- (61) Rosswog S, Feindt U, Korobkin O, Wu M R, Sollerman J, Goobar A and Martinez-Pinedo G 2017 Class. Quant. Grav. 34 104001 (Preprint eprint 1611.09822)
- (62) Pian E et al. 2017 Nature (Preprint eprint 1710.05858)
- (63) Maselli A, Cardoso V, Ferrari V, Gualtieri L and Pani P 2013 Phys. Rev. D88 023007 (Preprint eprint 1304.2052)
- (64) Yagi K and Yunes N 2017 Phys. Rept. 681 1–72 (Preprint eprint 1608.02582)
- (65) Metzger B D, Piro A L and Quataert E 2009 Mon. Not. Roy. Astron. Soc. 396 304 (Preprint eprint 0810.2535)
- (66) Lee W H, Ramirez-Ruiz E and Diego-Lopez-Camara 2009 Astrophys.J. 699 L93–L96 (Preprint eprint 0904.3752)
- (67) Martin D, Perego A, Arcones A, Thielemann F K, Korobkin O and Rosswog S 2015 Astrophys. J. 813 2 (Preprint eprint 1506.05048)
- (68) Wu M R, FernÃ¡ndez R, MartÃnez-Pinedo G and Metzger B D 2016 Mon. Not. Roy. Astron. Soc. 463 2323–2334 (Preprint eprint 1607.05290)
- (69) Lippuner J, FernÃ¡ndez R, Roberts L F, Foucart F, Kasen D, Metzger B D and Ott C D 2017 Mon. Not. Roy. Astron. Soc. 472 904–918 (Preprint eprint 1703.06216)
- (70) Fujibayashi S, Sekiguchi Y, Kiuchi K and Shibata M 2017 Astrophys. J. 846 114 (Preprint eprint 1703.10191)
- (71) Fujibayashi S, Kiuchi K, Nishimura N, Sekiguchi Y and Shibata M 2018 Astrophys. J. 860 64 (Preprint eprint 1711.02093)
- (72) Siegel D M and Metzger B D 2018 Astrophys. J. 858 52 (Preprint eprint 1711.00868)
- (73) Metzger B D, Thompson T A and Quataert E 2018 Astrophys. J. 856 101 (Preprint eprint 1801.04286)
- (74) Meszaros P 2006 Rept. Prog. Phys. 69 2259–2322 (Preprint eprint astro-ph/0605208)
- (75) Lee W H and Ramirez-Ruiz E 2007 New J. Phys. 9 17 (Preprint eprint astro-ph/0701874)
- (76) Giacomazzo B, Perna R, Rezzolla L, Troja E and Lazzati D 2013 Astrophys. J. 762 L18 (Preprint eprint 1210.8152)
- (77) Ascenzi S, De Lillo N, Haster C J, Ohme F and Pannarale F 2018 arXiv: 1808.06848
- (78) Bernuzzi S, Radice D, Ott C D, Roberts L F, Moesta P and Galeazzi F 2016 Phys. Rev. D94 024023 (Preprint eprint 1512.06397)
- (79) Bernuzzi S, Nagar A, Dietrich T and Damour T 2015 Phys.Rev.Lett. 114 161103 (Preprint eprint 1412.4553)
- (80) Dietrich T, Bernuzzi S, Ujevic M and Brügmann B 2015 Phys. Rev. D91 124041 (Preprint eprint 1504.01266)
- (81) Dietrich T, Ujevic M, Tichy W, Bernuzzi S and Brügmann B 2017 Phys. Rev. D95 024029 (Preprint eprint 1607.06636)
- (82) Dietrich T, Radice D, Bernuzzi S, Zappa F, Perego A, BrÃ¼gmann B, Chaurasia S V, Dudi R, Tichy W and Ujevic M 2018 arXiv: 1806.01625 (Preprint eprint 1806.01625)
- (83) Dietrich T and Hinderer T 2017 Phys. Rev. D95 124006 (Preprint eprint 1702.02053)
- (84) Healy J, Lousto C O and Zlochower Y 2017 Phys. Rev. D96 024031 (Preprint eprint 1705.07034)
Appendix A GW170817 analysis
We build our analysis of GW170817 on the publicly available posteriors released by the LVC. We proceed in three steps: (i) we review the original samples, (ii) we incorporate relations to ensure that both stars employ the same EOS, (iii) we restrict the viewing angle based on state-of-the-art GRB and afterglow models.
I. Original LIGO posteriors: For GW170817, we use the published results of the LVC TheLIGOScientific:2017qsa (); Abbott:2018wiz (); Abbott:2018exr (). In particular, the publicly available posterior samples of Abbott:2018wiz () provide the starting point for our GW interpretation. We decide to employ the “low spin” assumption, which restricts the rotational frequency of individual NSs so that the individual dimensionless spins are restricted to . This restriction is motivated by the observed BNS in our galaxy where the fastest spinning NS in a BNS system (PSR J1946+2052 Stovall:2018ouw ()) will have a dimensionless spin of at merger. Thus, we use the samples low_spin_PhenomPNRT_posterior_samples.dat.gz of https://dcc.ligo.org/LIGO-P1800061. The analysis of the follow up LVC results Abbott:2018wiz () improves over the initial results of the initial LVC results TheLIGOScientific:2017qsa (), as a broader frequency band of Hz, recalibrated Virgo data, more sophisticated waveform models Dietrich:2018uni (), and the known source location from EM observations have been used Monitor:2017mdv ().
II. Quasi-universal relations: While Abbott:2018wiz () did not make any assumption of the nature of the two merging compact objects, providing a general analysis, it seems natural to assume that the merging objects are two NSs described by the same EOS, an assumption similar to De:2018uhw (); Abbott:2018exr (); Radice:2018ozg ()111We note the caveat that the employed quasi-universal relations do not incorporate first-order phase transitions.. We follow Abbott:2018exr () by assuming GW170817 was caused by the merger of two NSs described by the same EOS. However, instead of using the relations presented in Yagi and Yunes Yagi:2015pkc () as employed in Abbott:2018exr () and Chatziioannou:2018vzf (), we use the the relation from De:2018uhw (). To do so, we discard all samples for which the quasi-universal relation is violated by more than . Note that we do not employ directly the posterior samples provided together with Abbott:2018exr (), as the publicly available samples do not contain the viewing angle. Furthermore, although an overall uncertainty of the quasi-universal relation of was estimated in De:2018uhw () we increase this value by a factor to be more conservative in our analysis.
As discussed in the literature Chatziioannou:2018vzf (); De:2018uhw (); Abbott:2018wiz (), enforcing the two objects to be NSs described by the same EOS leads to slightly tighter constrains on the tidal deformability (see the top panel of Fig. 4). We show as a dashed line the original posterior (without assuming the quasi-universal relations) and the solid line marks the posterior incorporating the quasi-universal relation. Based on this result, we see that the GW data do not support . This motivates that in the following analysis of the EM counterparts we will restrict the tidal deformabilities to .
III. The viewing angle: As a last step to restrict the GW posterior samples, we incorporate a viewing angle constraint based on the work of vanEerten:2018vgj () (note that also other works analyzing GRB170817A and its afterglow can be employed and lead to similar results Finstad:2018wid (); Mooley:2017enz ()). Ref. vanEerten:2018vgj () finds that the viewing angle of the GRB170817A was degrees. For a conservative estimate, we assume degrees to allow for additional uncertainties. All posterior samples with viewing angles that do not fall inside this interval are discarded. The final result is shown as the shaded region in Fig. 4.
Appendix B Analyzing AT2017gfo
Similar to the analysis of GW170817, we proceed in multiple steps to obtain a posterior distribution of the binary properties of AT2017gfo. For this purpose, we follow and extend the discussion in Coughlin:2018miv () to which we refer for further details and extensive checks of the underlying algorithms.
I. Modeling AT2017gfo: The observational data Coughlin:2018miv (); Smartt:2017fuw (); GBM:2017lvd () (see Fig. 5) are fit with the radiative transfer model of Kasen:2017sxr (). The model employs a multi-dimensional Monte Carlo code to solve the multi-wavelength radiation transport equation for an expanding medium. The model allows for the usage of a two component ejecta description. Each component depends parametrically on the ejecta mass , the mass fraction of lanthanides , and the ejecta velocity . These individual parameters will depend on the merger process and the binary parameters. The usage of at least two components is motivated by the presence of different ejecta types. Among the biggest drawbacks of our analysis is the assumption that both components are treated spherically symmetric with a uniform composition. Neglecting mixing of different ejecta types Rosswog:2016dhy (), we add the two separate model components. We have tested the recoveries of the non-spherical models presented in Kasen:2017sxr () using the spherical model. Based on these data, there will be a viewing angle bias in the parameter recoveries here depending on the degree of asymmetry in the ejecta.
We employ a grid with ejecta masses = 0.001, 0.0025, 0.005, 0.0075, 0.01, 0.25, 0.05, and 0.1, ejecta velocities = 0.03, 0.05, 0.1, 0.2, and 0.3, and mass fraction of lanthanides = 0, , , , , and . In order to draw inferences about generic sources not corresponding to one of these gridpoints, we adapt the approach outlined in Doctor:2017csx (); Purrer:2014fza (), where GPR is employed to interpolate principal components of gravitational waveforms. The reliability of the method has been tested in Coughlin:2018miv ().
|AT2017gfo||GRB170817A- vanEerten:2018vgj ()||GRB170817A- Wu:2018bxg ()||GRB170817A- Wang:2018nye ()|
|[0,1100]||KN analysis||KN analysis||KN analysis|
|[1,2]||KN analysis||KN analysis||KN analysis|
|[2.0,2.17]||KN analysis||KN analysis||KN analysis|
|50.30 0.84||0.81 0.39|
|KN analysis||KN analysis|
For completeness, we present the lightcurves together with the observational data in Fig. 5. The posterior for the ejecta properties is shown on the left of Fig. 6. The priors for the analysis are given in Table 2. We find that we are able to fit the observational data within the assumed magnitude uncertainty Coughlin:2018miv (). We want to point out that although we are able to fit and describe the general X-shooter spectra Smartt:2017fuw (); Pian:2017gtc (), the current model is unable to accurately represent observed wavelength specific features (see the discussion in Coughlin:2018miv ()).
II. Relating ejecta properties to the binary parameters:
To connect the individual ejecta components to the different ejecta mechanisms, we assume that the first ejecta component is proportional to dynamical ejecta, i.e., it gets released during the merger process and is proportional to . The second ejecta component arises from disk winds. We find that constraints on the mass ratio mostly follow from this first assumption, and constraints on the tidal deformability arise mainly from the second ejecta component. While the analysis has velocity and lanthanide fraction priors to make these components physically motivated, in the case where these assumptions are incorrect, the analysis will break down.
With the uncertainty in the modeling of ejecta in numerical relativity simulations and the potential systematic biases due to the missing input physics, we only assume that the dynamical ejecta describes a fraction of the total first component:
To allow for a direct comparison with the GW analysis, we express in terms of . This can be achieved by writing the compactnesses of the individual stars as Maselli:2013mva (); Yagi:2016bkt (), employing again , and using the definition of the tidal deformability
The second ejecta component is related to ejecta arising from disk winds. Long-term simulations find that about Dessart:2008zd (); Metzger:2008av (); Metzger:2008jt (); Lee:2009uc (); Fernandez:2013tya (); Siegel:2014ita (); Just:2014fka (); Metzger:2014ila (); Perego:2014fma (); Martin:2015hxa (); Wu:2016pnw (); Siegel:2017nub (); Lippuner:2017bfm (); Fujibayashi:2017xsz (); Fujibayashi:2017puw (); Siegel:2017jug (); Metzger:2018uni (); Radice:2018xqa () of the overall disk mass can be ejected. Thus, it seems plausible to assume
i.e., the disk wind ejecta are overall proportional to the mass of the debris disk surrounding the remnant BH. Knowing that a large fraction of the disk falls into the BH directly after BH formation and that not all matter gets ejected, we restrict to lie within .
The right part of figure 6 shows the findings of our AT2017gfo analysis,
which we shortly summarize below:
(i) our study favors equal or nearly equal mass systems,
where the constraint on the mass ratio mainly arises from the
correlation between the first component ejecta and the dynamical ejecta.
(ii) shows a clear jump at .
This constraint arises mainly from the second component ejecta and is related
to the disk mass increase for larger values of .
(iii) Only about of the first ejecta component is associated to
(iv) About of the disk mass has to be ejected to account for the
second ejecta component, which agrees with the disk wind ejecta found
in long term numerical relativity simulations.
(v) The analysis shows no strong constraint
on the maximum allowed TOV mass.
Appendix C Analyzing GRB170817A
In addition to including information about the viewing angle to restrict the GW posteriors (Fig. 4), we will also present a Bayesian parameter estimation for GRB170817A directly.
To relate the GRB properties to the properties of the binary, we employ the typical assumption that the GRB is driven by the accretion of matter from the debris disk onto the BH Eichler:1989ve (); Paczynski:1991aq (); Meszaros:1992ps (); Narayan:1992iy (); Meszaros:2006rc (); Lee:2007js (); Giacomazzo:2012zt (); Ascenzi:2018mwp () and therefore the energy is proportional to the disk rest mass, i.e.,
We note that based on our previous discussion, a fraction of the disk is ejected by disk winds. This part of the original disk cannot drive the GRB, so we set
To connect the GRB and kilonova analysis, we reuse the posterior obtained in the previous subsection. Similarly, we also employ the posterior distributions of as priors for our future GRB parameter estimation analysis.
I. The structured GRB-jet model of van Eerten et al vanEerten:2018vgj ():
van Eerten et al vanEerten:2018vgj () show that the latest observations of the GRB170817A afterglow is consistent with the emergence of a relativistic structured jet (with a jet width of degrees) seen at an angle of degree from the jet axis. This structured jet model fits well within the range of properties of cosmological sGRBs. Incorporating the Gaussian structured form of the jet as proposed in vanEerten:2018vgj (), the final GRB energy arising from the measured isotropic energy is given by
According to the analysis of vanEerten:2018vgj (), one obtains . We use this result as an input for Eq. (10) and sample over the final value employing a Gaussian distribution with a width identical to the stated uncertainty in vanEerten:2018vgj ().
The left panel of Fig. 7 shows our results. We find that , i.e., of the disk rest-mass is converted into GRB energy. This generally agrees with existing theoretical studies Lee:2007js (); Giacomazzo:2012zt () and increases the confidence in our GRB analysis. Furthermore, we find that the estimate and the constraint on the mass ratio shifts to slightly larger values than studying purely AT2017gfo.
II. The GRB model of Wu and MacFadyen Wu:2018bxg (): An alternative description of GRB170817A is presented in Wu and MacFadyen Wu:2018bxg (). They employ the analytic two-parameter “boosted fireball” model. The model consists of a variety of outflows varying smoothly between a highly collimated ultra-relativistic jet and an isotropic fireball. Developing a synthetic light curve generator, they fit the observational data by performing a Markov-Chain Monte Carlo (MCMC) analysis. Similar to vanEerten:2018vgj (), Wu and MacFadyen Wu:2018bxg () favor a relativistic structured jet. The jet opening angle is degrees seen from a viewing angle of degree.
The middle panel of Fig. 7 shows our results. The estimated GRB energy is , i.e., more than an order of magnitude below the estimated GRB energy of van Eerten et al vanEerten:2018vgj (). Consequently, the estimated value of is smaller. Nevertheless, the constraints on the binary properties and the EOS constraints are in agreement between van Eerten et al vanEerten:2018vgj () and Wu and MacFadyen Wu:2018bxg (), i.e., the constraints are robust to the systematic difference in energy estimates.
III. GRB due to the Blanford-Znajek mechanism Wang:2018nye ():
As a final way to interpret the observed GRB, we follow Wang:2018nye (). In this model, the energy to launch the GRB, assuming neutrino annihilation as the central engine, requires massive disks masses of the order of . Such massive disks are in tension with state-of-the-art numerical relativity simulations and disfavor neutrino annihilation as the mechanism responsible for the jet-launch. On the other hand, magnetic energy extraction requires disk masses about one order of magnitude smaller and therefore could act as the central engine for GRB170817A. Following the discussion of Wang:2018nye (), the disk mass necessary to explain the observation of GRB170817A based on the Blanford-Znajek (BZ) mechanism is given by
In contrast to Wang:2018nye (), we substitute by a single parameter for which we assume a flat prior with . Furthermore, we extend the analysis of Wang:2018nye (), who used a flat distribution for the BH spin within , by employing Eq. (19) to estimate the BH spin, and we substitute the disk mass fits of Radice:2018pdn () by Eq. (14). As in the previous discussions, we also incorporate the disk wind ejecta via Eq. (10). The right panel of Fig. 7 shows our results. The final results on the tidal deformability, mass ratio, and are similar to the previous results.
Very recently, Ref. Wang:2018nye () provided constraints on the EOS obtained from a new interpretation of the GRB and its afterglow phase, quoting a constraint of . Our tests show that this constraint is highly dependent on the particular choice of made in Wang:2018nye (). Assuming flat priors on all unknown parameters in Eq. (12) creates a prior peaking at . This prior choice is the driving mechanism for the very tight constraint presented in Wang:2018nye () and seems in our opinion to be an artifact of the sampling rather than a physical observation.
Appendix D Fits to numerical relativity
We now present the fits to numerical relativity we performed which are required for our analyses.
I. Disk mass A crucial ingredient for the analysis in this work and also the recent works of Radice:2018ozg (); Wang:2018nye () is the estimate of the debris disk mass . Here, we revisit the derivation of the phenomenological fit presented in Radice:2018pdn () and employ their fiducial dataset to derive an updated version of the fit. Figure 15 of Radice:2018pdn () shows a clear correlation between and the tidal deformability of the binary. We suggest that the reason for this clear and prominent correlation is related to the limited sample of only four EOSs in comparison to the wide range of sampled masses, and the fact that the tidal deformability depends strongly on the total mass of the binary, De:2018uhw (). Already from Fig. 15 of Radice:2018pdn (), one sees that setups with the same EOS (and hence roughly same radii) but different NS masses, lead to different disk masses.
Here we propose an alternative explanation which naturally accounts for the observed phenomenology and scaling with . Merger simulations suggest that the disk mass is accumulated primarily through redistribution of matter in the post-merger remnant. Thus, the remnant lifetime prior to collapse is found to strongly correlate with the amount of disk mass Radice:2018xqa (). Here we suggest that this lifetime is governed to a large degree by , where is the threshold mass above which the merger undergoes a prompt-collapse (on dynamical timescales). Thus is a measure of the stability of the post-merger remnant, and following the arguments above, should correlate with .
We show in Fig. 2 in the main text the correlation between the disk mass and the threshold mass for prompt BH formation, where we estimate the prompt collapse threshold as Bauswein:2013jpa ():
denotes the maximum mass of a non-rotating (TOV) NS for a given EOS and is the radius of a star. With the help of Fig. 2 in the main text, it becomes clear that the reduction of the disk mass relates to the stability of the merger remnant and consequently, the disk mass drops abruptly when and the remnant undergoes a prompt-collapse. This naturally explains the location of the turnover in .
Based on these observations and the fact that the NS radius can be related to the tidal deformability by (with the chirp mass ) De:2018uhw (), we conclude that the disk mass is a function of the tidal deformability, the total mass of the system, and the maximum TOV mass . We do not find a strong dependence on the mass ratio and neglect mass ratio effects since the dataset does not resolve those effects well. Therefore, information about the densest part of the EOS, encoded in , and the information at smaller densities, encoded in or , are essential for a reliable description of the disk mass.
To include the dependence of , we also find that fitting instead of leads to a significant reduction of the fractional error . We choose the following functional form
The mean absolute error of with respect to the original numerical relativity data is , and we obtain a fractional error of in ; for comparison, the original fit presented in Radice:2018pdn () has absolute errors of and average fractional errors of . The large fractional error is caused by a small number of setups with very small disk masses. Fitting the logarithm of the disk mass improves the fit in this region of the parameter space, as already discussed in Coughlin:2018miv () for the dynamical ejecta. We present the fit as a function of the in Fig. 2 in the main text and also present the absolute (middle panel) and fractional errors (bottom panel) for Eq. (14) in comparison with the results of Radice:2018pdn (). We point out that in a region around the turning point into prompt collapse scenarios, the absolute and fractional errors of the new fit are noticeable smaller than in the original version.
Finally, since we want to relate information extracted from the disk mass estimates with the GW measurement, we propose to relate the NS radius to the tidal deformability via De:2018uhw ()
While informing this relation adds an additional uncertainty, we find that the fitting residuals increase simply to and .
II. Dynamical Ejecta We approximate the mass of the dynamical ejecta by
with , , , and and denoting the compactnesses of the individual stars. The absolute uncertainty of the fit, i.e., is . Furthermore, we note that while the fractional error of is only , the fractional error with respect to is caused by datapoints with very small ejecta masses ().
The velocity of the dynamical ejecta is given by
where , , and . The average absolute error of the fit is and the fractional error is .
III. BH properties The large set of numerical relativity data publicly released in the CoRe catalog Dietrich:2018phi () together with results published in Bernuzzi:2015opx () allows us to derive phenomenological fits for the BH mass and spin. A detailed list of the employed simulations and the BH properties is presented in Tab. 3. We restrict our consideration to non-spinning NSs, but plan to extend the presented results in the future once a larger set of spinning BNS configurations is available. Furthermore, we consider only cases for which an almost stationary state is reached after BH formation, so that remnant properties can be extracted reliably. Thus, we do not consider setups for which the BH mass increases significantly due to accretion or for which the black hole mass decreases due to insufficient resolution.
Trivially, we find that with an increasing total mass, the final black hole mass and angular momentum increases almost linearly. For unequal mass mergers, and decrease. Based on this observation, we propose a functional dependence of (where refers to the symmetric mass ratio). The coefficient is chosen to be two, which is motivated by predictions for BBH systems Healy:2017mvh (). Considering the imprint of the EOS, we find that for larger values of , the final black hole mass decreases, which follows from the observation that the disk mass increases with . As a simple ansatz, we choose:
with and . The mean average absolute error of the fit is and the fractional error is .
We find in our dataset that the BH mass and spin are strongly correlated. This motivates the use of a similar functional behavior for the BH spin as for the BH mass. However, we extend Eq. (18) by (i) adding an additional constant, and (ii) incorporating the fact that the dimensionless spin is restricted to be . The final fitting function is
with , , and . The mean average absolute error of the fit is and the fractional error is . The remnant property dataset and the corresponding fit results are presented in Fig. 8.