The Challenges Ahead for Multimessenger Analyses of Gravitational Waves and Kilonova: a Case Study on GW190425
Geert Raaijmakers, Samaya Nissanke, Francois Foucart, Mansi M. Kasliwal, Mattia Bulla, Rodrigo Fernandez, Amelia Henkel, Tanja Hinderer, Kenta Hotokezaka, Kamil? Lukošiūt?, Tejaswi Venumadhav, Sarah Antier, Michael W. Coughlin, Tim Dietrich, Thomas D. P. Edwards
DDraft version February 24, 2021
Typeset using L A TEX twocolumn style in AASTeX62
The Challenges Ahead for Multimessenger Analyses of Gravitational Waves and Kilonova: a Case Study onGW190425
Geert Raaijmakers, Samaya Nissanke,
1, 2
Francois Foucart, Mansi M. Kasliwal, Mattia Bulla, Rodrigo Fern´andez, Amelia Henkel, Tanja Hinderer,
1, 7
Kenta Hotokezaka, Kamil ˙e Lukoˇsi¯ut ˙e, Tejaswi Venumadhav,
9, 10, 11
Sarah Antier, Michael W. Coughlin, Tim Dietrich,
14, 15 andThomas D. P. Edwards
16, 1 GRAPPA, Anton Pannekoek Institute for Astronomy and Institute of High-Energy Physics, University of Amsterdam, Science Park 904,1098 XH Amsterdam, The Netherlands Nikhef, Science Park 105, 1098 XG Amsterdam, The Netherlands Department of Physics, University of New Hampshire, 9 Library Way, Durham NH 03824, USA Division of Physics, Mathematics, and Astronomy, California Institute of Technology, Pasadena, CA 91125, USA The Oskar Klein Centre, Department of Astronomy, Stockholm University, AlbaNova, SE-106 91 Stockholm, Sweden Department of Physics, University of Alberta, Edmonton, AB T6G 2E1, Canada Institute for Theoretical Physics, Utrecht University, Princetonplein 5, 3584 CC Utrecht, The Netherlands Research Center for the Early Universe, Graduate School of Science, University of Tokyo, Bunkyo-ku, Tokyo 113-0033, Japan Department of Physics, University of California at Santa Barbara, Santa Barbara, California 93106, USA School of Natural Sciences, Institute for Advanced Study, 1 Einstein Drive, Princeton, New Jersey 08540, USA International Centre for Theoretical Sciences, Tata Institute of Fundamental Research, Bangalore 560089, India Universit´e de Paris, CNRS, Astroparticule et Cosmologie, F-75013 Paris, France School of Physics and Astronomy, University of Minnesota, Minneapolis, Minnesota 55455, USA Institute of Physics and Astronomy, University of Potsdam, Karl-Liebknecht-Str. 24/25, 14476, Potsdam, Germany Max Planck Institute for Gravitational Physics (Albert Einstein Institute), Am M¨uhlenberg 1, Potsdam 14476, Germany The Oskar Klein Centre, Department of Physics, Stockholm University, AlbaNova, SE-106 91 Stockholm, Sweden (Received XXXXXX; Revised XXXXXX; Accepted February 24, 2021)
ABSTRACTIn recent years, there have been significant advances in multi-messenger astronomy due to the discov-ery of the first, and so far only confirmed, gravitational wave event with a simultaneous electromagnetic(EM) counterpart, as well as improvements in numerical simulations, gravitational wave (GW) detec-tors, and transient astronomy. This has led to the exciting possibility of performing joint analysesof the GW and EM data, providing additional constraints on fundamental properties of the binaryprogenitor and merger remnant. Here, we present a new Bayesian framework that allows inferenceof these properties, while taking into account the systematic modeling uncertainties that arise whenmapping from GW binary progenitor properties to photometric light curves. We extend the relativebinning method presented in Zackay et al. (2018) to include extrinsic GW parameters for fast analysisof the GW signal. The focus of our EM framework is on light curves arising from r-process nucle-osynthesis in the ejected material during and after merger, the so called kilonova, and particularly onblack hole - neutron star systems. As a case study, we examine the recent detection of GW190425,where the primary object is consistent with being either a black hole (BH) or a neutron star (NS).We show quantitatively how improved mapping between binary progenitor and outflow properties,and/or an increase in EM data quantity and quality are required in order to break degeneracies in thefundamental source parameters. INTRODUCTION
Corresponding author: G. [email protected]
The first gravitational-wave detection of two mergingcompact objects, during the first observing run of Ad-vanced LIGO, gave us a new direct probe into the fun-damental physics governing such objects (Abbott et al.2016). In the subsequent observing runs of LIGO andVirgo, tens of binary black hole (BBH) mergers, and sev- a r X i v : . [ a s t r o - ph . H E ] F e b Raaijmakers et al. eral binary neutron star (BNS) and black hole - neutronstar (BHNS) mergers have been reported, allowing forpopulation level studies (Abbott et al. 2019, 2020c). Sofar, however, the BNS merger GW170817 (Abbott et al.2017a) (combined signal-to-noise ratio (SNR) of 32.4)has been the only event for which electromagnetic (EM)signatures were detected across the frequency spectrum(e.g., Abbott et al. 2017b,c; Coulter et al. 2017; Kasliwalet al. 2017; Hallinan et al. 2017). The extensive follow-up campaign of the counterpart resulted in additionalinformation on many aspects of the merger, rangingfrom the dense matter equation of state (EOS) of NSs(e.g., Radice et al. 2018c; Raaijmakers et al. 2020; Tewset al. 2020) to measurements of the Hubble constant(e.g., Abbott et al. 2017d; The LIGO Scientific Collab-oration et al. 2019; Hotokezaka et al. 2019; Mukherjeeet al. 2019; Howlett & Davis 2020; Nicolaou et al. 2020;Dietrich et al. 2020; Dhawan et al. 2020).Recent GW events potentially involving NSs —GW190425 (Abbott et al. 2020a), GW190426 (Abbottet al. 2020c) and GW190814 (Abbott et al. 2020d) —have all led to similar large-scale searches for an EMcounterpart (Coughlin et al. 2019b; Goldstein et al.2019; Andreoni et al. 2020; Antier et al. 2020; Pageet al. 2020; Gompertz et al. 2020). These searches havenot resulted in any successful detections so far, thoughthe absence of a counterpart can add weak additionalconstraints on the binary parameters of the system(Coughlin et al. 2020a,b; Anand et al. 2020; Ackleyet al. 2020). The absence of a detection could be due tothe low SNR of some of these events, e.g. GW190425with a combined SNR of 12.9, producing much largerskymaps compared to GW170817.One of the counterparts to a BHNS or BNS merger isthe so-called kilonova or macronova (Li & Paczy´nski1998; Kulkarni 2005), a thermal ultraviolet-optical-infrared transient (UVOIR) arising from r-process nu-cleosynthesis (Burbidge et al. 1957; Cameron 1957) inthe neutron rich material that is ejected during andafter the merger (Lattimer & Schramm 1974; Rosswoget al. 1999; Metzger et al. 2010). If a kilonova is de-tected, critical questions remain over whether one couldimprove upon the constraints of the binary parametersand possibly uncover the ambiguous nature of one ofthe binary components in e.g., GW170817 (see, e.g.,Hinderer et al. 2019; Coughlin & Dietrich 2019) andGW190425 (see, e.g., Barbieri et al. 2020a; Most et al.2020; Kyutoku et al. 2020). Leaving out the unconfirmed detection of a candidate elec-tromagnetic counterpart to a binary black hole merger (Grahamet al. 2020).
There have been a number of studies performing amultimessenger analysis of GW170817 and the associ-ated kilonova, most of which focus on the dense mat-ter EOS, using either a limit on the minimum amountof ejected mass derived from the UVOIR light curves(e.g., Radice et al. 2018c; Radice & Dai 2019; Hindereret al. 2019; Raaijmakers et al. 2020; Capano et al. 2020)or a joint Bayesian modeling framework (e.g., Cough-lin et al. 2018; Coughlin et al. 2019a; Dietrich et al.2020; Breschi et al. 2021b). The methods in Coughlinet al. (2019a) and Dietrich et al. (2020) allow for addi-tional uncertainty in the analysis, attributed to the un-certainty in the modeling of ejected material, by addingan unknown component to the kilonova. There havealso been more general studies developing a Bayesianframework for multimessenger analysis, focusing on sev-eral parameters describing the binary progenitor system(e.g., Barbieri et al. 2019; Breschi et al. 2021a; Nichollet al. 2021). Coughlin et al. (2019a), Dietrich et al.(2020) and Barbieri et al. (2019) also consider the shortgamma-ray burst (GRB) signal (e.g., Goldstein et al.2017) associated with compact object mergers, whichcan help constrain the inclination angle of the system(e.g., Ryan et al. 2020). Barbieri et al. (2019) further-more consider the radio remnant due to the kilonova andGRB afterglow (see e.g. Hotokezaka et al. 2016). How-ever, Barbieri et al. (2019) do not include uncertaintiesin their ejecta modeling or light curve computation whenperforming a joint analysis of the GW and EM data.More recently, there have been two new Bayesian anal-yses frameworks developed by Breschi et al. (2021b) andNicholl et al. (2021). In Breschi et al. (2021b), the au-thors use a three-component model with angular depen-dence, first presented in Perego et al. (2017), and per-form a joint analysis of GW170817 and the associatedkilonova. In Nicholl et al. (2021), the authors use thepublic light curve generation code in
MOSFiT (Guillo-chon et al. 2018) and divide the ejecta into two compo-nents. Additionally they allow for a third componentarising from shock-heated material by a GRB jet. Theythen perform a similar joint analysis of GW170817 andthe associated kilonova, focusing on extracting the densematter EOS.In this work, we will present a new multimessengerframework to analyse gravitational waves and the associ-ated UVOIR light curves jointly. We leave the inclusionof the GRB and radio counterpart to future work, andfurthermore will not include any interaction between thekilonova material and the jet (see e.g. Kasliwal et al.2017; Piro & Kollmeier 2018; Nativi et al. 2021; Klionet al. 2021). We will use the analytical formulae, fit-ted to numerical simulation of BNS and BHNS, derivedin Kr¨uger & Foucart (2020) to connect outflow prop-erties to binary progenitor properties. We incorporatemodeling uncertainties by estimating the errors on theseformulae and by using a new relation that estimates theamount of material that is ejected from the disk sur-rounding the merger remnant, which is assumed to be aBH. Because of this assumption, our framework is onlyappropriate for BNS systems with high total masses andfor BHNS systems; BNS systems with lower total massesmight lead to a NS or hypermassive NS remnant whichwill significantly alter the light curve properties (see,e.g., Kawaguchi et al. 2020). The paper is structured asfollows; in Section 2, we will discuss how the propertiesof the binary progenitor are connected to the propertiesof the outflows during and after the merger; in Section3, we will discuss the kilonova model that is used tocompute UVOIR light curves from these outflow prop-erties; in Section 4, we will test our analysis frameworkon a low SNR GW190425-like merger with simulatedGW and EM data, and discuss the implications in Sec-tion 5. We note that the model presented in this workincludes a more careful consideration of the errors forthe BHNS scenario than the BNS scenario, which weleave to future work. CONNECTING BINARY PROGENITOR TOOUTFLOW PROPERTIESIn order to build a multimessenger analysis frame-work, we start by connecting the fundamental sourceproperties of the binary progenitor to parameters thatdescribe average properties of the outflows that occurduring and after the merger (see also Figure 1). Boththe connecting relations and the source properties ofinterest differ between a BNS and a BHNS system.For a BNS system, the parameters of interest are themasses of the two neutron stars and their tidal de-formabilities (see, e.g., Hinderer et al. 2010). In thiswork, the spin of the neutron star is assumed to besmall when computing outflow properties and thus willnot impact the kilonova. This assumption is based onthe spin distribution of Galactic BNSs, which indicatesthat the magnitude of the dimensionless spin parameter χ = cIω/ ( GM ) ≤ .
05 (Zhu et al. 2018). In general,however, the spin of the neutron star may have a signif-icant effect on the outflows (see, e.g. Most et al. 2019;East et al. 2019; Chaurasia et al. 2020). For a BHNSsystem, the parameters of interest for the neutron starare again the mass and the tidal deformability. Blackholes, however, are thought to have a tidal deformabilityof zero in general relativity (see e.g. Binnington & Pois-son 2009; Chia 2020; Chirenti et al. 2020). Contrary toneutron stars, however, the spins of black holes in bina- ries can be as high as | χ BH | (cid:46) . C NS = GM NS / ( c R NS ) (where G is the Gravitational constant, c is the speed of light, M NS and R NS the mass andradius of the NS respectively), as a function of its tidaldeformability, Λ NS : C NS = 0 . − . NS ) + 0 . NS ) . (1)While the accuracy of this relation is on the order of10% for nucleonic EOS, we will assume it to hold trueas it evades the need to directly integrate an EOS toobtain a radius. 2.1. Dynamical ejecta
In the case of a BNS system, the mass of the dynamicalejecta is estimated using the formula derived in Kr¨uger& Foucart (2020): M BNSdyn − M (cid:12) = (cid:18) aC + b M n M n + cC (cid:19) M + (1 ↔ , (2)where the coefficients are a = − . , b = 114 . , c = − .
56, and n = 1 . Raaijmakers et al.
Light curvesOutflow propertiesBinary properties { Fit formulae to numerical simulations { Semi-analytical light curve model { Figure 1.
Schematic drawing of how binary progenitor properties connect to outflow properties and UVOIR light curves. Ifwe consider a sample of binary properties (cid:126)x , obtained from, e.g., parameter estimation on the GW strain data, then each pointin this sample ( (cid:126)x , (cid:126)x . . . ) maps to outflow properties (cid:126)y via analytical formulae calibrated to numerical simulations. Becauseof the uncertainties in both the numerical simulations and in these formulae, one sample in (cid:126)x gives a range of possible outflowproperties. From each of these points in (cid:126)y , light curves can be computed through semi-analytical models, where again theuncertainty in these models leads to multiple light curves for a given point in (cid:126)y . that M > M . We assume the maximum error to bethe 2 σ standard deviation obtained from the residualsbetween the fitting formula and the numerical simula-tions, which corresponds to an error of ± .
008 M (cid:12) .When we assume that the progenitor is a BHNS sys-tem (with a mass ratio of Q = 1 /q = M BH /M NS ), we usethe following equation (Kr¨uger & Foucart (2020), basedon models by Kawaguchi et al. (2016)) to estimate thedynamical ejecta: M BHNSdyn M b NS = a Q n − C NS C NS − a Q n R ISCO M BH + a . (3)Here, M b NS is the baryonic mass of the NS and R ISCO isthe radius of the innermost stable circular orbit arounda black hole of mass M BH (Bardeen et al. 1972). Thefirst can be approximated, following Lattimer & Prakash(2001), by the simple formula, M b NS = M NS (cid:18) . C NS − . C NS (cid:19) , (4)while the latter is given by, R ISCO = 3 + Z − sgn( χ BH ) (cid:112) (3 − Z )(3 + Z + 2Z ) , (5) Z = 1 + (1 − χ ) / (cid:104) (1 + χ BH ) / + (1 − χ BH ) / (cid:105) , (6) Z = (cid:113) χ + Z . (7) In this scenario, we again assume the error in the massof the dynamical ejecta to be the 2 σ limit of the numer-ical error of the fit, i.e. ∆M BHNSdyn = 0 . (cid:12) . Wenote, however, that for both the BNS and BHNS sce-nario, this is an underestimate of the real uncertaintydue to the uncertainties in numerical modeling.The velocity of the dynamical ejecta and the error inthe velocity are assumed to be equal in the BHNS andBNS systems, which we justify by assuming that for theBNS system the dynamical ejecta is dominated by thetidal tail (see next paragraph). The velocity and theerror are approximated using the formulae derived inFoucart et al. (2017): v dyn = 0 . Q + 0 . , (8)∆ v dyn = 0 . v dyn . (9)Finally, the composition of the dynamical ejecta is ex-tremely neutron rich, which leads to a large lanthanide-fraction through rapid neutron capture. Since lan-thanides have many absorption lines, this will cause theejecta to be very opaque, and so we will use an effec-tive grey opacity for this component between 1 and 10cm g − for both BHNS and BNS systems. However,for (especially symmetric) BNS mergers, the dynamicalejecta includes shock heated ejecta on top of the tidaltail which are predicted to have a lower lanthanide frac- Independent of frequency. Q . . . . . . ξ = M e j / M d i s k Thermally-driven outflowsMagnetically-driven outflows M disk = 0 . disk = 0 . disk = 0 . disk = 0 . disk = 0 .
32 4 6 8 10 12 14 16 M BH (M (cid:12) ) Figure 2.
Uncertainty in the fraction of the initial diskM disk that is eventually ejected as a function of the mass ra-tio of the progenitor binary Q . The points show the results ofthe viscous hydrodynamic simulations of BH accretion disksfrom Fern´andez et al. (2020), which capture the thermally-driven outflow due to angular momentum transport (i.e., dis-sipation of MRI turbulence) and nuclear recombination. Thenumbers in the plot indicate the disk compactness C d andthe top axis is the black hole mass used for the initial setup.As an estimation, we connect to the mass ratio by dividingthe black hole mass by a fiducial neutron star mass of 1 . (cid:12) . The light-shaded band indicates the uncertainty cap-tured by Equation (12), while the darker-shaded band showsthe possible enhancement in the fraction of the disk massejected through magnetic driving, ranging between 0% and20% depending on the initial field geometry of the disk (c.f.Christie et al. 2019) tion, and thus a lower opacity (see, e.g., Sekiguchi et al.2016; Dietrich & Ujevic 2017; Tanaka et al. 2020; Nedoraet al. 2021). For simplicity, we assume that this compo-nent does not contribute significantly compared to thetidal tail. For the case study considered in Section 4, thisassumption can be justified by the asymmetric massesof the simulated signal.2.2. Disk wind ejecta
The second type of ejecta comes from matter outflowsin the accretion disk surrounding the merger remnant. Itis important to note that the following discussion onlyapplies when the merger remnant is a BH, which forBNS systems is only the case when the total mass is highenough that the merger remnant undergoes collapse to aBH. The first step to estimate the mass of these outflowsis to compute the mass of the disk, i.e. the material thatis still bound to the merger remnant and supported byrotation against collapse. In the case of a BNS systemwe use the formula derived in Kr¨uger & Foucart (2020) to compute: M BNSdisk = M max (cid:8) × − , ( aC + c ) d (cid:9) . (10)For the BHNS scenario, we first compute the remnantmass outside of the black hole after the merger using thefitting formula of Foucart et al. (2018):ˆ M rem = (cid:20) Max (cid:18) α − NS η / − β C NS η , (cid:19)(cid:21) δ , (11)where η is the symmetric mass ratio, i.e. η = Q/ (1+ Q ) .From the remnant mass, we can compute the mass ofthe disk surrounding the black hole by simply subtract-ing the dynamical ejecta, i.e. M BHNSdisk = ˆ M rem − M BHNSdyn .A certain fraction of the disk will not be accreted ontothe black hole but will instead become gravitationallyunbound in what we will call a disk wind. Numericalsimulations of a disk surrounding a BH merger rem-nant indicate that the ejected material from the diskcan roughly be divided into two further sub-componentsbased on their velocity. The first thermally-driven com-ponent arises from energy dissipation due to angularmomentum transport (i.e. MRI-driven turbulence) andnuclear recombination on timescales longer than the dy-namical time, which can be modeled well by viscoushydrodynamic simulations (see, e.g., Fern´andez et al.2020). A faster magnetically-driven component appearson shorter timescales in general relativistic magnetohy-drodynamic (GRMHD) simulations (see, e.g., Siegel &Metzger 2017; Fern´andez et al. 2019; De & Siegel 2020).For the thermally-driven outflow component generatedaround a BH remnant, Fern´andez et al. (2020) foundan inverse linear dependency of the ejected mass onthe compactness of the disk using viscous hydrodynamicsimulations. A disk formed closer to the BH (more com-pact) ejects a smaller fraction of the initial disk massthan an initially more extended disk because a largerfraction of the disk is accreted before the thermal out-flow can begin, given its proximity to R ISCO . We assumethat the compactness of the disk is roughly proportionalto the mass ratio of the progenitor binary so that we canestimate the fraction of the disk that is ejected thermallyas a function of the initial binary properties. In Figure2, we show the fraction ξ of ejected material M ej overthe initial disk mass M disk found in Fern´andez et al.(2020). On the top axis, we show the value of the blackhole mass used in the simulations, while the numbersin the plot indicate the disk compactness C d = ( M BH / (cid:12) )(50km /R d ), where R d is the radius of the disk. Theinitial conditions used in Fern´andez et al. (2020) are in-spired by numerical relativity simulations using neutronstars of mass M NS ∼ . − . M (cid:12) . Thus, we estimate the Raaijmakers et al. mass ratio by simply dividing the mass of the black holeby a fiducial neutron star with mass 1 .
35 M (cid:12) . To cap-ture the uncertainty in the viscous outflows as a functionof the mass ratio, we use the ansatz: ξ = M ej M disk = ξ + ξ − ξ e . Q − , (12)where for the lower bound ξ = 0 .
04 and ξ = 0 . ξ = 0 .
32 and ξ = 0 . ∼ . ≤ . c . In Fern´andez et al. (2019), the au-thors find that the material ejected due to the magneticfield has velocities ≥ . c , which is expected consideringtheir poloidal magnetic field configuration. The velocityof the thermally-driven outflows is much lower ( ≤ . c ),with the bulk of the material around 0 . c − . c (seealso Fern´andez et al. 2020). We take the average veloc-ity to be between 0 . c and 0 . c , with a minimum andmaximum velocity distribution cutoff above 0 . v wind andbelow 2 v wind respectively (see Table 1).The composition of the disk outflows is less accuratelyknown than for the dynamical ejecta. The thermal com-ponent is ejected on timescales long enough that it ismostly reprocessed by neutrino emission and absorp-tion. The composition is primarily light r -process ele-ments, with variable amounts of lanthanide-rich matterdepending on detailed properties of the disk and cen-tral object (e.g., Just et al. 2015; Martin et al. 2015;Wu et al. 2016; Lippuner et al. 2017). The composi-tion of the magnetic component is less well understoodgiven the limited number of long-term GRMHD simula-tions. The faster overall velocities imply shorter expan-sion times and less neutrino reprocessing, with a larger . . . M (M (cid:12) ) . . . . . . . . M ( M (cid:12) ) . . . . . . P ( B HN S | G W ) Figure 3.
Posterior samples of the two source masses inGW190425 for the high-spin priors using the
IMRPhenomDNRT waveform model. The samples are color coded by the prob-ability of GW190425 being a BHNS or a BNS system, giventhe constraints on the EOS by GW170817 (Abbott et al.2018). expected proportion of lanthanide-rich matter than thethermal component (Siegel & Metzger 2017; Fern´andezet al. 2019). Neutrino absorption has been included inone GRMHD study only, with results of the early diskevolution indicating that it is indeed an important pro-cess in regulating the composition (Miller et al. 2019).To model the opacity we choose an effective grey opacitybetween 0 . − . g − for the thermal componentand between 1 −
10 cm g − for the magnetic component.An important note to make here is that in the case ofthe merger remnant being a hypermassive neutron star(HMNS), both the fraction of the disk that is ejectedand the composition change drastically. Due to absorp-tion of neutrinos emitted by the HMNS and absence ofa BH mass sink, the disk could be entirely ejected, de-pending on the lifetime of the HMNS (see, e.g., Metzger& Fern´andez 2014; Fujibayashi et al. 2018; Fahlman &Fern´andez 2018; Fujibayashi et al. 2020; M¨osta et al.2020). We do not consider this possibility here, but willleave it for future work.2.3. Ejecta mass distributions GW190425
To better understand the formulae discussed in Sec-tion 2.2, we apply them to the gravitational wave eventGW190425 (Abbott et al. 2020a). GW190425 was de-tected in the LIGO Livingston and the Virgo detectors.The two component masses, also shown in Figure 3,are consistent with the system being either a BNS ora BHNS, depending on the maximum mass of NSs asdetermined by the dense matter EOS. − − log ( M dyn ) − − − − − l og ( M w i nd ) − − log ( M wind ) BNS
This workFormulae in Barbieri et al. (2019)Formulae in Dietrich et al. (2020) − − − − log ( M dyn ) − . − . − . − . l og ( M w i nd ) − . − . − . − . log ( M wind ) BHNS
This workFormulae in Barbieri et al. (2019)
Figure 4.
Comparison of ejecta masses for the dynamical component and the disk wind component for the BNS scenario (leftpanel) and the BHNS scenario (right panel). The masses are computed by using different sets of fitting formulae (Coughlinet al. 2019a; Barbieri et al. 2020b) on the high-spin posterior distribution of GW190425 (Abbott et al. 2020a). The shadedregions contain respectively 68% and 95% of all ejecta masses after discarding the samples where M dyn or M wind is zero. Thefractions of discarded samples are 58%, 78%, and 30% for the BNS formulae in this work, in the work of Barbieri et al. (2019),and in the work of Dietrich et al. (2020), respectively. For the BHNS scenario, the fractions of discarded samples are 40% and25% for the formulae in this work and in the work of Barbieri et al. (2019), respectively.
In order to apply the formulae to GW190425, werandomly draw O (10 ) samples from the public poste-rior probability density function (PDF) for GW190425.For each draw, we simultaneously take a random EOSsample from the public posterior EOS samples ofGW170817 (Abbott et al. 2018). Depending on whetherthe primary mass of our draw is above or below the max-imum mass of the chosen EOS, we can thus assume thenature of the primary component of GW190425 is aBH or NS respectively. The dimensionless BH spin, χ BH , needed in case of a BHNS system is part of thesample that is drawn from GW190425. However, thetidal deformability of the two objects in GW190425 wasunconstrained, so instead we set the tidal deformabilityof the NS(s) by evaluating the drawn EOS at the valuesdrawn for the mass(es). We then apply Equations 2 -11 to compute the corresponding ejecta masses where,for the disk wind ejecta, we randomly pick a value for ξ between the upper and lower bound of Equation (12)to determine how much of the disk is unbound. https://dcc.ligo.org/LIGO-P2000223/public, high-spin priorwith waveform model IMRPhenomDNRT https://dcc.ligo.org/LIGO-P1800115/public, low-spin priorwith waveform model IMRPhenomPNRT
The resulting distributions of ejecta masses are shownin Figure 4 in pink, where the different shades illustratethe regions containing 68% and 95% of all values. Forcomparison, we also show the distributions one wouldcompute when using different sets of fitting formulae. Inbrown, we show the distributions following the methodfrom Barbieri et al. (2019), where the authors use fitsfor the dynamical ejecta and disk mass from Kawaguchiet al. (2016) and Foucart et al. (2018) respectively in aBHNS merger and from Radice et al. (2018b) and Barbi-eri et al. (2019) respectively in the case of a BNS merger.For the BHNS scenario, there is reasonable agreementbetween the distribution using the formulae in Barbi-eri et al. (2019) and the distribution computed in thiswork. This is due to the fact that the same formula forthe disk mass is used, and the formula for the dynamicalejecta in Kr¨uger & Foucart (2020) is an adjusted versionof the formula in Kawaguchi et al. (2016). For the BNSscenario, the formula used in Barbieri et al. (2019) pre-dicts a lower dynamical ejecta mass. The predictionsfor disk wind mass is consistent between Barbieri et al.(2019) and this work. In green we also show the dis-tributions obtained from using the same method as inDietrich et al. (2020) for BNS systems. The authors useformulae for the dynamical ejecta mass and disk wind
Raaijmakers et al. mass derived in Coughlin et al. (2019a) and Dietrichet al. (2020), respectively. Although the predictions fordynamical ejecta mass are slightly higher than in thiswork, there is general agreement between the two distri-butions. LIGHT CURVE MODELINGFrom the inferred ejecta parameters discussed in theprevious section, we can generate the correspondinglight curves. To do so, we use the publicly availablecode by Hotokezaka & Nakar (2020), a semi-analyticmodel based on the work by Li & Paczy´nski (1998). Incontrast to other semi-analytic light curve models (see,e.g. Waxman et al. 2019), both the effect of the time de-lay between photon production and photon emergenceand the gradient in the velocity are incorporated. Wenote, however, that the code does not take into accountthe asymmetric geometry of the outflows (see, e.g., Bar-bieri et al. 2019; Bulla 2019; Kawaguchi et al. 2020; Zhuet al. 2020; Heinzel et al. 2021, for viewing angle de-pendent light curve models) which can have a signifi-cant impact on the luminosity and color of the observedlight curve (Bulla 2019; Darbha & Kasen 2020; Korobkinet al. 2020). We also fix, for an increase in computa-tional efficiency, the nuclear heating rate to the modelby Korobkin et al. (2012) (with thermalization efficiency (cid:15) = 0 . dyn , v dyn , and M wind computedfrom the posterior PDFs of GW190425, we estimate thecorresponding light curve. Since the software by Ho-tokezaka & Nakar (2020) does not only take an averagevelocity as input but also requires a velocity distribu-tion, which we have chosen to be a power law distribu-tion, we have to consider a few extra parameters. Theseconsist of a power law index n , varied between 3 . − . . dyn and v dyn , and the maximum ve-locity cutoff between 1 . dyn and 2 . dyn . For the post-merger disk winds, the average velocity of the ejectais not connected to the properties of the binary but istaken to be between 0 . c and 0 . c . The minimum veloc-ity cutoff is the same as for the dynamical component, https://github.com/hotokezaka/HeatingRate The code by Hotokezaka & Nakar (2020) does include a cal-culation of the thermalization based on the ejecta properties. while the maximum velocity cutoff is slightly lower, i.e.,between 1 . wind and 2 . wind .The final parameters that we need to include are theeffective opacities of the ejected material. The dynam-ically ejected material is very neutron-rich and will,therefore, produce lanthanides through r-process nucle-osynthesis, leading to a high opacity. We allow theeffective grey opacity for this component to be withinthe range of 1 −
10 cm g − . Since the post-mergerdisk winds consist of ejecta with varying composition,we assign two effective grey opacities, one in the range0 . − . g − for v ≤ v κ and one in the range 1 − g − for v > v κ . The transition velocity v κ is ran-domly drawn between the minimum and maximum ofthe velocity distribution. Table 1 shows a summary ofall parameters and their allowed ranges.The final light curves are then obtained by adding thetwo components, where estimates for different photo-metric bands rely on the assumption of blackbody emis-sion. The uncertainty in the distance is also includedby randomly drawing a distance measurement from theposterior PDF of GW190425. In Figures 5 and 6, weshow the light curves corresponding to posterior PDFdistributions of GW190425 in the g-band and r-band,for both the BNS scenario (left panels) and the BHNS(scenario). The dark and light shaded regions indicatethe regions where 68% and 95% of the light curves arecontained respectively. The black horizontal line showsthe upper limits obtained by the GROWTH collabora-tion using ZTF (Coughlin et al. 2019c).As a comparison, and to better assess the currentmodeling uncertainties, we also include light curves ob-tained with our code when the input parameters are de-rived from other methods (upper panels). First, we com-pute light curves from the distribution of ejecta massesand velocities following the methods outlined in Cough-lin et al. (2019a) with updated fitting formulae fromDietrich et al. (2020), which are tailored to BNS sys-tems. In Coughlin et al. (2019a) the authors also usea two-component model, divided into dynamical anddisk wind ejecta. They relate the dynamical ejectamass to the total ejected mass in the first componentby asserting the proportionality M ej , = α M dyn , where0 . ≤ α ≤
1. For the second component, the mass is as-sumed to be a fraction of the disk mass M ej , = ξ M disk , Again with the caveat that more symmetric NS mergers willalso produce shock heated material with lower opacities (see Sec-tion 2.1). The temperature is estimated through the bolometric lu-minosity L bol and the radius of the photosphere r ph as T = L bol / (4 πσr ). Table 1.
All parameters used in the model described in Section 2 and their prior support in the full analysis of the GW signaland it’s EM counterpart. The notation U ( a, b ) here means uniformly drawn between boundaries a and b . We also show theadjusted priors for the Case 4 run discussed in Section 4.2Parameters Description Prior density and support (Case 4) Binary properties M [M (cid:12) ] Mass of the primary object ∈ P ( x | ˜d GW190425 )M [M (cid:12) ] Mass of the secondary object ∈ P ( x | ˜d GW190425 ) χ Spin parameter of the primary object (BHNS) ∈ P ( x | ˜d GW190425 )Λ Tidal deformability of the primary object Λ(M ; EOS GW )Λ Tidal deformability of the secondary object (BNS) Λ(M ; EOS GW ) Ejecta and light curve properties M dyn [M (cid:12) ] Mass of the dynamical ejecta Eq. (2) (BNS) or (3) (BHNS)v dyn [c] Velocity of the dynamical ejecta Eq. (9) v min , dyn [c] Minimum velocity of the dynamical ejecta ∼ U(0.1, 1.0) v dyn ∼ U(0.7, 0.9) v dyn v max , dyn [c] Maximum velocity of the dynamical ejecta ∼ U(1.5, 2.5) v dyn ∼ U(1.5, 1.7) v dyn n dyn Power law index of density distribution ∼ U(3.5, 4.5) κ dyn [cm g − ] Effective grey opacity of the dynamical ejecta ∼ U(1.0, 10.0) ∼ U(5.0, 8.0)M wind [M (cid:12) ] Mass of the disk wind ejecta Eq. (10) (BNS) or (11) (BHNS) and Eq. 12v wind [c] Velocity of the disk wind ejecta ∼ U(0.1, 0.3) ∼ U(0.12, 0.18) v min , wind [c] Minimum velocity of the disk wind ejecta ∼ U(0.1, 1.0) v wind ∼ U(0.3, 0.5) v wind v max , wind [c] Maximum velocity of the disk wind ejecta ∼ U(1.5, 2.0) v wind ∼ U(1.5, 1.7) v wind n wind Power law index of density distribution ∼ U(3.5, 4.5)v κ [c] Transition velocity between low and high κ ∼ U(v min , wind , v max , wind ) κ low [cm g − ] Effective grey opacity for v ≤ v κ ∼ U(0.1, 1.0) ∼ U(0.1, 0.3) κ high [cm g − ] Effective grey opacity for v > v κ ∼ U(1.0, 10.0) ∼ U(6.0, 10.0) where 0 ≤ ξ ≤ .
5. The velocity of the second com-ponent and the lanthanide fraction of both componentsare not related to any binary parameters, but consid-ered as free parameters in the ranges 0 . ≤ v disk ≤ . − ≤ log ( X lan ) ≤ −
1. The black dashed (pinkdashed) lines in Figure 5 (5) illustrate the region where95% of the light curves are contained. Since the parame-ters α and ξ allow for a larger range of ejecta masses, theresulting kilonova light curves are relatively broad, butconsistent with the light curves presented in this paper.Second, we show the light curves computed with ourcode using the ejecta properties obtained by followingthe methods in Barbieri et al. (2019). The authors usea similar approach to the methods presented in this pa-per, though with a three-component model; a dynamicalcomponent, a neutrino-driven disk wind component, anda viscous (secular, corresponding to ‘thermally-driven’in our framework, Equation 12) component. From thedisk mass, they estimate that for a BHNS system, 1%is ejected through neutrino-driven winds while 20% isejected through viscous processes, both with a veloc-ity of 0 . . .
04c respectively.An effective grey opacity is used, which is set to 15 cm g − for the dynamical ejecta and to 1 and 5 cm g − for respectively the wind and viscous ejecta in a BHNSscenario. For BNSs, only the opacity of the disk windis lowered to 0 . g − . The authors do include anangle dependence in their light curve computations butthis is not included here. The effect of including an-gle dependencies is larger for the BHNS scenario, wherethe geometry of the dynamical ejecta is more likely to beconfined to the equatorial plane and not in all azimuthaldirections. As a result the kilonova appears brighterwhen we observe in the polar direction and fainter whenwe observe in the equatorial plane. The yellow (black)dot-dashed lines in Figure 5 (6) again indicate the re-gion where 95% of all light curves are contained. Forthe BHNS scenario, these regions are comparable to ourresults, mostly due to the fact that the same underly-ing fitting formula for the disk mass is used (see Fou-cart et al. 2017). Their method does, however, predictfainter light curves for the BNS scenario. Although Fig-ure 4 suggests that the disk wind mass is comparableto our results, there is a large fraction of samples not0 Raaijmakers et al. A B m ag n i t ud e BNS σ , based on Barbieri et al. (2020)2 σ , based on Dietrich et al. (2020)This workupper limits ZTF NSBH σ , based on Barbieri et al. (2020) . . . . . time (days) Bulla (2019), equatorial viewBulla (2019), polar view . . . . . time (days) Kyutoku et al. (2020)Bulla (2019), equatorial viewBulla (2019), polar view
Figure 5.
The light curves in the g-band from assuming that GW190425 is a BNS system (left panels) or a BHNS system (rightpanels), given the equation of state constraints from GW170817. The different shaded bands indicate the regions that contain68% and 95% of all light curves respectively. Also shown in the upper panels are the 95% regions for light curves computedusing the code by Hotokezaka & Nakar (2020), but using the methods in Barbieri et al. (2020b) (orange dashed-dotted lines)and Dietrich et al. (2020) (black dashed line). Note that, contrary to the upper bounds, the lower bounds of the 95% credibleregions should not be taken definitively as part of the parameter space consistent with GW190425 results in (almost) no EMradiation. The brown solid lines in the lower right panel are two models from Kyutoku et al. (2020); one for a lanthanide-richdisk outflow and dynamical outflow at 140 Mpc and one for a lanthanide-poor disk outflow and dynamical outflow at 130 Mpc.The blue (grey) dashed lines in the lower panels correspond to an equatorial view (polar view) of light curves computed with
POSSIS (Bulla 2019), with ejecta masses set by 0 . ≤ M dyn , M wind ≤ .
03 M (cid:12) for BHNS, 0 . ≤ M dyn ≤ .
02 M (cid:12) and 0 . ≤ M wind ≤ .
03 M (cid:12) for BNS, and distances in the range 120 −
200 Mpc. The black horizontal lines correspond to the upper limitsfound with ZTF by the GROWTH collaboration, covering 21% integrated probability of the skymap for this event. (Coughlinet al. 2019c). shown, where the dynamical ejecta mass is zero and thedisk wind mass is very low, that dominate the resultingregion containing 95% of all light curves. Another factorthat contributes to the differences in the resulting lightcurves are the higher opacities chosen for dynamical anddisk wind ejecta.We also show light curves that are computed using dif-ferent approaches to generate light curves. In the caseof a BHNS system, we show the light curves computedby Kyutoku et al. (2020) (shown in their Figure 1). Theauthors use a Monte-Carlo radiation-transfer code (seeTanaka & Hotokezaka 2013; Kawaguchi et al. 2020) to estimate light curves from a representative disk massand dynamical mass outflow. Instead of using an effec-tive grey opacity, as in this work, the opacity is time-dependent and derived with updated atomic structurecalculations of r-process elements (Tanaka et al. 2020).The variation in their light curves is a combination of dif-ferent viewing angles, composition, and distances. Thelower limit shown here corresponds to their equatorialview of a lanthanide-rich disk outflow and dynamicalejecta at 140 Mpc, while the upper limit correspondsto a polar view of a lanthanide-poor disk outflow withdynamical ejecta at 130 Mpc.1 A B m ag n i t ud e NSNS σ , based on Barbieri et al. (2020)2 σ , based on Dietrich et al. (2020)This workupper limits ZTF NSBH σ , based on Barbieri et al. (2020) . . . . . time (days) Bulla (2019), equatorial viewBulla (2019), polar view . . . . . time (days) Kyutoku et al. (2020)Bulla (2019), equatorial viewBulla (2019), polar view
Figure 6.
Same as Figure 5, but for the r-band.
Lastly, we show light curves computed with
POSSIS (Bulla 2019), a three dimensional Monte-Carlo codethat models the radiation transport in kilonovae, tak-ing into account the geometry of the merger outflowsand the viewing angle. For faster light curve genera-tion the code does not solve the full radiative transferequation but rather uses time- and wavelength depen-dent opacities as direct input. For the BHNS case, weuse the grid presented in Anand et al. (2020) and re-strict to ejecta masses that are consistent with the dis-tribution of ejecta masses computed in Section 2, i.e.M dyn = 0 . , . , .
03 M (cid:12) and M wind = 0 . , . , . (cid:12) . The distance is varied between 120 −
200 Mpc, cor-responding to the 1 σ range of the luminosity distancemeasure inferred in Abbott et al. (2020a). We show theresulting bands for two viewing angles, one in the equa-torial and one in the polar plane. Similarly, for the BNScase, we use the grid computed in Dietrich et al. (2020)and restrict to M dyn = 0 . , . , . , .
02 M (cid:12) andM wind = 0 . , .
03 M (cid:12) .We conclude that although the methods to computethe light curves in Figure 5 and 6 are quite different,there is broad agreement on the photometric magnitude range that is predicted for GW190425. This is due tothe large statistical uncertainties in the GW posteriorPDF, and specifically the large spread in distance thatdirectly affects the magnitude of the light curve. Thereis, however, some discrepancy between the light curvescomputed from the ejecta properties based on the fittingformulae in Barbieri et al. (2019) and the formulae in Di-etrich et al. (2020) and this work. This is again a resultof the lower predicted M dyn and M wind in this particu-lar part of parameter space using the formulae in Radiceet al. (2018b) and Barbieri et al. (2019). When compar-ing the light curve ranges between BHNS and BNS it isnot clear that an EM detection would be able to distin-guish the nature of the system, but this warrants a moredetailed study which we will leave for future work. MULTIMESSENGER PARAMETERESTIMATION4.1.
Gravitation wave analysis
For fast analysis of the GW signal, we use the rela-tive binning method introduced by Zackay et al. (2018),which exploits the fact that the difference between grav-2
Raaijmakers et al. .
487 1 . M c D L . . . χ e ff ˜ Λ . . . q . . . . M . . . M M . . . . M . . . q ˜Λ . . . χ eff
100 200 D L BHNS priorBNS prior, high spinBNS prior, low spinGW190425
Figure 7.
Posterior distributions of binary properties inferred from a gravitational wave signal. The grey contours correspondto the posteriors obtained in Abbott et al. (2020a) with a high spin prior. The colored contours are the posteriors obtainedwith relative binning from our simulated signal. The red, blue and green colors indicate a BHNS prior, a BNS prior with highspins and a BNS prior with low spins respectively. The dashed grey lines indicate the injected values. itational waveforms with nonzero posterior probabilitydensity in the frequency domain can be described by asmoothly varying perturbation. By using the ratio ofa gravitational waveform with some fiducial waveformclose to where the likelihood peaks, the number of fre-quency points where the waveform is evaluated reducesto O (10 ), compared to O (10 ) for traditional GW pa-rameter inference. However, Zackay et al. (2018) did notinclude any extrinsic geometric parameters (see, e.g., Abbott et al. 2020b) in their method so here we haveextended the relative binning method used in our codeto include parameters such as distance, sky location, in-clination, and polarization (see also Finstad & Brown2020, for a similar approach).We simulate a waveform similar to GW190425(see Table 2 for the source parameters), using the IMRPhenomD NRTidalv2 model (Dietrich et al. 2019),and inject this into the LIGO Livingston and Virgo de-3 . . . . . time (days) A B m ag n i t ud e Posterior drawsBHNS priorBNS priorSimulated data
Figure 8.
The green points indicate the EM data pointsin the g-band used for the analysis in Section 4, while theshaded bands indicate the 95% credible regions for the priors,assuming either a BHNS or a BNS system. In black we showa selection of posterior draws from the posterior distributionobtained in Case 2 (see main text Section 4.2).
Table 2.
Source properties of the simulated gravitationalwave signal.Parameter Injected value M c (source frame) 1.441 M (cid:12) M c (detector frame) 1.4871 M (cid:12) M (source frame) 2.25 M (cid:12) M (source frame) 1.24 M (cid:12) q 0.55Λ χ ,z χ ,z -0.02Distance 145 Mpc tectors with an SNR of 11.6 and 4.7 respectively. Wenote that the mass of the primary object, M , is withinthe range of possible NS masses (see e.g. Raaijmakerset al. 2020), but by setting the tidal deformability to beΛ = 0, we set this object to be a BH. Furthermore,it is important to mention that the spin of the BH isoutside of the range of the low-spin prior, which waschosen as a test case to see how the typically employedprior choice of LVC analyses might biases the parameterestimation if the type of the progenitor system is not We note that there exist now publicly available waveformmodels dedicated to BHNS signals that were not available whenthis work started (see, e.g., Matas et al. 2020; Thompson et al.2020). certainly determined. To generate the noise, we use arealistic Power Spectral Density (PSD) from O3a andassume it to be stationary and Gaussian. We then com-pute the summary data necessary for relative binningand set up our likelihood function. If we assume thenoise in the detectors to be Gaussian and defined by apower spectrum S ( f ), we can write the likelihood as, L ( d GW | θ GW ) ∝ (13)exp (cid:20) − (cid:104) d GW − h ( θ GW ) | d GW − h ( θ GW ) (cid:105) (cid:21) , where d GW is the GW strain data, h ( θ GW ) a grav-itational waveform defined by the vector set of modelparameters θ GW , and (cid:104) a | b (cid:105) = 4 (cid:60) (cid:90) f max f min ˜ a ∗ ( f )˜ b ( f ) S ( f ) df. (14)We use Bayes’ theorem to express, P ( θ GW | d GW ) ∝ π ( θ GW ) L ( d GW | θ GW ) , (15)where P ( θ GW | d GW ) is the posterior distributionon the gravitational wave parameters given the data,and π ( θ GW ) is the prior distribution on these param-eters. To sample from the posterior PDF distribution,we use the nested sampling algorithm MultiNest (Ferozet al. 2009, 2013; Buchner et al. 2014). Instead of sam-pling in component masses, we sample in the parame-ters chirp mass and mass ratio to speed up the conver-gence of our posterior distribution. We choose a uniformprior in both, within the range M c ∈ [1 . , . q = 1 /Q ∈ [0 . , , ∈ [0 , = 0. Theprior for the two dimensionless spins is uniform in mag-nitude, while the orientation is isotropically distributed,although here we only consider spins that are alignedwith the orbital angular momentum to speed up the pa-rameter inference. When assuming a BNS system, weeither choose a low spin, | χ | ∈ [0 , . | χ | ∈ [0 , . | χ | ∈ [0 , . | χ | ∈ [0 , . π ( D L ) ∝ D L , with bounds cho-sen to be far from where the posterior probability isnon-negligible. The priors for the right ascension and https://dcc.ligo.org/LIGO-T2000012/public Raaijmakers et al. declination are taken uniformly across the sky, and thepriors for the inclination and polarization are uniformlydistributed in angle from 0 to 2 π or π , respectively. Weshow the posterior PDF distributions of all three priorchoices in Figure 7 and include the posterior distribu-tion of GW190425 for the high spin prior case in Abbottet al. (2020a) for comparison, which is the most agnosticcase when the nature of the binary is not known.Figure 7 shows the effects of changing the priors on thespins as seen by comparing the green and blue curves.When comparing the red to the other colored curves, italso shows the effects of the assumption on the natureof the binary, which involves different priors on tidal de-formabilities (for the BH it is set to zero) as well as onthe spins (high spin for the BH, low spin for the NS).As expected, the luminosity distance measured from theGW amplitude is largely unaffected by changes to thepriors on parameters that mainly impact the phasing.We note that even with the BHNS priors correspondingto the injected signal, the true values of the parametersare not correctly recovered due to the small SNR. Wealso note that the BHNS assumption actually does notyield results close to the injected values for any of the pa-rameters except for the tidal deformability. Instead, forthe masses and spin parameter the high-spin BNS priorslead to results closest to the true values. As expected,the mass and spin measurements are most impacted bythe spin priors, as seen by comparing the blue and greencurves. This is a known degeneracy due to the fact thateffects of the mass ratio and spins first enter the phas-ing only separated by half of a post-Newtonian orderand thus accumulate similarly with the GW frequency.Our results also show that the PDF for tidal deformabil-ity is most strongly dependent on the prior assumptionson the type of binary system. For reference, the resultsfrom Abbott et al. (2020a) with a high spin prior are alsoshown as the grey contours, which are in good agreementwith our corresponding simulated results depicted by theblue curves. Note that the GW signal considered herehas a low SNR, similar to GW190425. Therefore the fol-lowing discussion on additional constraining power fromEM observations will only apply to such low-SNR sig-nals. High-SNR signals, and signals detected across abroader frequency range by third generation GW detec-tors, already provide tighter constraints on binary pa-rameters, which affect the additional constraining powerof EM observations.4.2. Multimessenger analysis
In order to test how an EM counterpart to a GW eventcan constrain the properties of the merging binary, wecombine our EM model defined in Section 2 and 3 with the GW analysis mentioned above in Section 4.1. Takingthe values of our simulated GW event, we compute thecorresponding light curve by using Equations (2)-(11),and additionally setting v min , dyn = 0 . v dyn , v max , dyn =1 . v dyn , κ dyn = 5 cm g − for the dynamically ejectedcomponent and v wind = 0 . v min , wind = 0 . v wind , v max , wind = 1 . v wind , κ low = 0 . g − , κ high = 8 cm g − , and v κ = 0 .
18c for the disk wind ejecta. For the def-inition of these parameters, we refer back to Section 3 orTable 1. We evaluate this light curve at t = [0 . , , , P ( θ | d ) ∝ π ( θ ) L ( d | θ ) (16) ∝ π ( θ ) L ( d GW | θ GW , EM ) L ( d EM | θ EM ) . (17)The vector θ now contains both the parameters neededfor the GW waveform generation that affect the lightcurves, θ GW , EM , and the parameters to compute an EMlight curve, θ EM . For the GW likelihood, we assume thatit is proportional to the posterior distribution obtainedfrom the GW analysis described in the previous section: L ( d GW | θ GW , EM ) ∝ (cid:90) P ( θ | d GW ) d θ nuisance , (18)We marginalize over all parameters in the GW anal-ysis that are not connected to the light curve model( θ nuisance ) so that the vector θ GW , EM only contains {M c , q, χ BH , Λ NS , D } and {M c , q, Λ , Λ , D } for aBHNS and BNS system, respectively. One implicationof Equation (18) is that the prior distribution on theseparameters is already incorporated into the likelihood,which we need to take into account when choosing ourpriors. Since some of the parameters in θ EM depend onthe parameters θ GW , EM through Equations 2 - 11, wecan rewrite Equation 16 using a hierarchical approach: P ( θ | d ) ∝ π ( θ GW , EM ) π ( θ EM | θ GW , EM ) (19) × L ( d GW | θ GW , EM ) L ( d EM | θ EM ) . Lastly, we assume a Gaussian likelihood for the EMdata: L ( d EM | θ EM ) ∝ exp (cid:18) − ( d EM − L ( θ EM )) σ (cid:19) , (20)where d EM are the photometric measurements of thelight curve and L ( θ EM ) the computed light curve us-ing the model presented in this work. Initially we set5 . . . q Λ Λ Λ Λ BNS
GWGW+EM, σ EM = 0 . σ EM = 0 . ξ known . . . q − . . . . χ Λ Λ . . χ BHNS
GWGW+EM, σ EM = 0 . σ EM = 0 . ξ known Figure 9.
Posterior distribution on binary parameters given the GW data (red) using the high-spin BNS prior, and whenincluding an EM counterpart with either unknown parameter ξ (blue) or known (green). The grey dashed lines indicate thevalues of the binary properties used to simulate the GW and EM data (note that since Λ = 0 the underlying value is not visiblehere) . σ EM = 0 . σ EM = 0 .
01 mag.For our prior distributions, we choose a uniform priorover the parameters θ GW , EM in π ( θ GW , EM ) to minimizethe effect of imposing prior information twice. For theparameters θ EM that do not depend on θ GW , EM , we useuniform priors in the ranges described in Table 1. Forparameters in the distribution θ EM that do depend on θ GW , EM , i.e. the dynamical ejecta mass and velocity andthe disk wind mass, we take a uniform distribution in theuncertainty range defined for that parameter accordingto Equations 2 - 11. We assume that the detection ofthe EM counterpart leads to an identification of the hostgalaxy, which we simplify to fixing the distance to thetrue value D = 145 Mpc, both when computing the lightcurve magnitude and when evaluating L ( d GW | θ GW ).In Figure 8, we show the 90% credible region of the priorlight curves in the g-band for both the BHNS and BNSscenario, together with the simulated EM data.The posterior PDF distribution is again sampled with MultiNest and shown for the parameters of interest inFigure 9 when using the high spin BNS prior for thegravitational wave data and both the BNS prior (leftpanel) and BHNS prior (right panel) for the EM data.Three different analyses are shown where: (i) only theGW data is considered; (ii) both GW and EM data are considered with σ EM = 0 .
2; and (iii) both GW and EMdata are considered with σ EM = 0 . and ξ is fixed. Here-after we refer to these analyses as Case 1-3.For the posterior distributions using the BNS prior forthe EM data, we first note that including EM informa-tion (Case 2) puts significant constraints on the mass ra-tio of the system. This can be attributed to the strongcorrelation between the mass ratio and M BNSdyn , wherelower and higher values of q lead to too much or too littledynamical ejecta respectively. The tidal deformability isnot well constrained for the primary object, since thisonly has a small effect on the dynamical ejecta. How-ever, for the secondary object the tidal deformability isthe main parameter that relates to the disk winds andis therefore much more constrained. The double peakstructure is the result of the correlation with both thedisk wind and dynamical ejecta velocity. A higher tidaldeformability leads to a higher disk wind mass with lowvelocity, but to be consistent with the early EM data,this requires a higher dynamical ejecta velocity. On theother hand a lower tidal deformability leads to a lowerdisk wind mass with high velocity, and thus requires alower dynamical ejecta velocity. Fixing the ξ parameter(Case 3) only strengthens the double peak structure inΛ , as the higher values of Λ are more consistent withthe now better measured disk mass.In the right panel of Figure 9, where a BHNS prior isused to analyze the EM data, we note that for Case 26 Raaijmakers et al. . . . q − . . . . χ Λ Λ . . χ BHNS
GWGW+EM, σ EM = 0 . ξ known, tight priorGW+EM, σ EM = 0 . ξ known, extra EM dataGW+EM, σ EM = 0 . ξ known Figure 10.
Same as Figure 9, but now showing posteriordistributions where we decrease the prior volume for the EMparameters (see Table 1) (red), add EM data in the i- andz-band (blue) and decrease the uncertainty in the EM datapoints to σ EM = 0 .
01 (green). The dashed, grey lines indicatethe GW posterior for the high-spin BNS prior. and 3, where σ EM = 0 .
2, there are only slightly betterconstraints on the mass ratio than for Case 1, due tothe degeneracies between parameters in θ EM . The mostsignificant improvement is made in the Λ parameter.This is the result of both M disk and M dyn dependingsimilarly on Λ , meaning that low values will result intoo little ejected total mass, while high values will leadto too much ejected mass to be consistent with the EMdata. Lastly, we note that the EM data considered herehas very little additional constraining power on the di-mensionless spin of the black hole.To investigate whether we can improve on the con-straints on the binary parameters, we study three morecases where the BHNS prior on the EM data is used.These are: (i) a scenario where we put a tighter prior onall parameters in θ EM (see Table 1); (ii) a scenario whereon top of the EM data in the g- and r-band we also havedata in the i- and z-band at the same epochs; and (iii)a scenario where the EM data have σ EM = 0 .
01. Here-after, these are referred to as Case 4-6. Case 4 couldbe seen as a hypothetical future scenario where thereis a better understanding of the relation between thebinary properties and all outflow properties on top ofthe parameter ξ . Case 5 and 6 could help answer thequestion of whether the constraints improve more fromtaking data in other photometric bands or decreasing the uncertainty in the photometric points one alreadyhas.The posterior distributions corresponding to Case 4-6are shown in Figure 10. For all cases, the constraintson the binary parameters are tighter than in Case 2,although for Case 5, where data in more photometricbands are considered, the peak of the posterior still hasmore support for near-equal mass ratios away from theinjected parameter value. This is a result of both theposterior for the GW data alone having more supportfor near-equal mass ratios (see Figure 7) and the sys-tem considered here having a small amount of dynamicalejecta compared to disk wind ejecta. That means thatthe EM data can be well-fitted by models with only diskwind ejecta. This is illustrated in Figure 11, where theposterior distributions on the outflow properties showvery little support for dynamical ejecta larger than zero.Such models, with small dynamical ejecta compared todisk wind ejecta, are associated with more equal masssystems, which is why there is more posterior supportfor larger values of q . For Case 4 and 6, the data is lesswell fitted by models with only disk wind ejecta, causingthe true underlying values to be within the 68% credibleregions. For Case 4, this is due to some of the degener-acy between the dynamical and disk wind ejecta beingalleviated while for Case 6, the small error on the EMdata allows for the small amount of dynamical ejecta tobe distinguished from the disk wind ejecta. DISCUSSIONBased on the results presented in Figures 9, 10, and12, we conclude that a detection of an EM counterpartin the future can lead to improved constraints on thebinary progenitor parameters for systems similar to thelow-SNR, BHNS/BNS system considered here. How-ever, even in this idealized scenario where the EM datacan be exactly modeled within this work’s framework,the high-dimensionality of the problem and the exist-ing degeneracies between parameters impact how wellone can measure the fundamental properties of the pro-genitor. In this Section, we will first briefly summarizeour results and then discuss the important caveats andlimiting assumptions of our work.On one hand, when modeling the EM lightcurve, ifwe assume that the GW190425-like signal (Section 4) isa BHNS, the injected parameter values are within the95% credible region for all posterior distributions (rightpanel of Figure 9 and Figure 10). However, the posteriordistributions on the mass ratio have more support fornear-equal mass ratios when only considering GW data(case 1), which does not change when adding EM data(Cases 2,3 and 5). This is because the small amount of7 .
007 0 . M dyn κ w i nd , . . . . κ w i nd , . . . v w i nd . . M w i nd . . . M d i s k κ d y n . . . v d y n .
05 0 . v dyn κ dyn .
10 0 . M disk .
02 0 . M wind .
15 0 . v wind . . κ wind , κ wind , GW+EM, σ EM = 0 . σ EM = 0 . ξ knownGW+EM, σ EM = 0 . Figure 11.
Posterior distributions for outflow properties for Case 1, 2 and 6. The grey dashed lines indicate the injected valuesof outflow properties used to generate the EM data. We note that for Case 1 and 2 (red and blue lines) the dynamical ejectamass and velocity are mostly consistent with being zero. dynamical ejecta in our simulated model allows for theEM data to be well fitted with models that have enoughdisk wind ejecta, i.e., models with higher values of q .This is resolved in Cases 4 and 6, where the posteriordistribution on the mass ratio peaks around the injectedvalue (see also Figure 12). In all cases, the constraintson the tidal deformability improve the most, while con-straints on the spin of the BH improve the least, com-pared to Case 1 (only GW data). This implies thatcurrently the addition of EM data to GW190425-like events, where the tidal deformability is unconstrained(Abbott et al. 2020a), could give valuable informationon the dense matter EOS. From the results of Cases 5and 6, it is not yet clear whether it is preferred to takemore data in different photometric bands or to improvethe accuracy of the datapoints, which would need a moredetailed study of the full parameter space and is left tofuture work. As expected, the constraining power, how-ever, does increase in both cases.8 Raaijmakers et al. . . . . q Λ C a s e C a s e C a s e C a s e C a s e C a s e − . . . χ Figure 12.
Marginalized one-dimensional posterior distri-butions for the three binary parameters that connect to out-flow properties. We show all six cases considered here, de-scribed in the main text of Section 4.2, and conclude thatwith EM observations (Cases 2 - 5) the constraining powerincreases when comparing to using only GW data (Case 1).
On the other hand, when modeling the EM lightcurves, if we assume that the GW190425-like signal is aBNS (Section 4), the posterior distributions suggest thatthe injected binary parameters are not always within the68% credible region. Although this work focuses on theBNS framework in less detail, we note that parameterinference is highly dependent on the assumption of thenature of the system. GRMHD simulations by Mostet al. (2020) show that high-mass BNS systems lead toa fast dynamical ejecta component, whereas this compo-nent is absent for systems with similar parameters barthe primary object being a BH. In principle, this couldhelp distinguish the nature of the system. However, ifsuch a signal is not detected, both the GW data andEM data might look very similar (as shown in Figures5 and 6). In this case, one needs to be careful wheninterpreting results from parameter estimation analysesthat assume the nature of the system.Furthermore, there are simplifying assumptions in theframework, detailed below, that could have a signifi-cant impact on multimessenger parameter inference, andneed to be investigated in more detail in future work.Firstly, throughout this work, we assume that ejectaoutflows are spherically symmetric, which is a poorassumption especially for BHNS systems (see, e.g., Kawaguchi et al. 2016; Foucart et al. 2017). Dependingon the viewing angle with which the system is observed,the light curve can be brighter or fainter in differentphotometric bands (Bulla 2019; Darbha & Kasen 2020;Korobkin et al. 2020).Secondly, information on the viewing angle is poten-tially obtainable from the observation of a relativisticjet, which can improve parameter estimation on sys-tems with asymmetric outflow geometries (see, e.g.,Heinzel et al. 2021). Currently, the relativistic jet thatis launched in some mergers (see, e.g. Ruiz et al. 2021) isnot taken into account in the model presented here. Nei-ther is the interaction between the jet and the dynamicaland disk wind ejecta which can alter the kilonova lightcurve (Nativi et al. 2021; Klion et al. 2021; Nicholl et al.2021).Thirdly, although we have made the simplifying as-sumption that the nuclear heating rate is fixed to a spe-cific model derived by Korobkin et al. (2012), there existmultiple models for nuclear heating rates (see, e.g., Ho-tokezaka & Nakar 2020; Kasen & Barnes 2019). Theuncertainty in the nuclear physics input can lead to upto one order of magnitude uncertainty in the bolometricluminosity (Zhu et al. 2020).Fourthly, we note that the fitting formulae presentedin Section 2 do not cover the entire parameter space thatour parameter inference spans (Section 4). This meansthat the results presented here are dependent on inter-polations and extrapolations of existing numerical sim-ulations, which might be subject to change when moresimulations are available in the (near) future. One couldalso increase the number of free parameters in these fit-ting formulae, by making use of numerical simulationswhere, e.g., the spin of the neutron star (Ruiz et al. 2020)or the orbital eccentricity of the binary (Papenfort et al.2018) is taken into account.Fifthly, we stress that the framework developed here ismostly focused on BHNS systems, with an approximateapplicability to high-mass BNS systems when there isno short- or long-lived merger remnant.Finally, a more extensive parameter study is neededto investigate how EM observations of a larger diver-sity of binary mergers can constrain fundamental binaryparameters. In order to better understand modellingerrors, our kilonova model cross-comparison should beextended to incorporate parameter inference (see, e.g.Dietrich et al. 2020; Heinzel et al. 2021).In conclusion, EM observations add constrainingpower to the GW data for the low SNR, BHNS sys-tem considered in this work. Our results suggest thatthis constraining power will improve when numericalmodels of merger outflows, telescope sensitivity and9GW detectors progress (see, e.g., Cowperthwaite et al.2019; Chornock et al. 2019; Palmese et al. 2019). Cur-rently, significant effort is underway into adding moresophisticated treatments of microphysical processes tonumerical models, as well as increasing the coverage ofthe parameter space studied with these models (Baiotti& Rezzolla 2017; Shibata & Hotokezaka 2019; Diet-rich et al. 2020; Ciolfi 2020; Foucart 2020; Ruiz et al.2021). The Vera Rubin Observatory, an optical, wide-field telescope expected to begin operations in 2021(Ivezi´c et al. 2019), will have a large impact on tak-ing high-accuracy EM data and detecting GW eventsat larger distances, while telescopes as ZTF (Dekanyet al. 2020) and (near-future) dedicated GW follow-uptelescopes such as BlackGEM (Bloemen et al. 2016),GOTO (Gompertz et al. 2020) will increase the num-ber of detected EM counterparts. Concurrently, nextgeneration GW detectors, such as the Einstein Tele-scope (Maggiore et al. 2020) and the Cosmic Explorer(Reitze et al. 2019), will increase the SNR of GW detec-tions and widen the frequency range in comparison withLIGO/Virgo, allowing higher-precision measurementsof binary properties.We thank Kyohei Kawaguchi for providing us withtheir light curves and Andrew Williamson, MatthewLiska, Doosoo Yoon, Koushik Chatterjee, PhilippMoesta, Om Salafia,, Masaomi Tanaka and KentaKuichi for useful discussions. We thank Barbara Patri-celli and Leo Singer for their careful reading of thepaper and their comments. We are very grateful tothe the LIGO Scientific, Virgo and KAGRA Collab-orations for public access to their data products ofGW170817, GW190425. We also thank the GROWTHcollaboration for public access to their observationaldata products. GR, SMN, TH, KL are grateful for fi- nancial support from the Nederlandse Organisatie voorWetenschappelijk Onderzoek (NWO) through the Pro-jectruimte and VIDI grants (Nissanke). TH also ac-knowledges funding from the NWO sectorplan. FFand AH gratefully acknowledge support from NASAthrough grant number 80NSSC18K0565, from the DOEthrough Early Career Award DE-SC0020435, and fromthe NSF through grant number PHY-1806278. MMKacknowledges the GROWTH project funded by theNational Science Foundation under PIRE Grant No1545949. MMK acknowledges generous support fromthe David and Lucille Packard Foundation. MB ac-knowledges support from the Swedish Research Council(Reg. no. 2020-03330). RF acknowledges supportfrom the Natural Sciences and Engineering ResearchCouncil of Canada (NSERC) through Discovery GrantRGPIN-2017-04286, and from the Faculty of Science atthe University of Alberta. TV acknowledge support bythe John Bahcall Fellowship at the Institute for Ad-vanced Study and by the National Science Foundationunder Grant No. 2012086. M.C. acknowledges supportfrom the National Science Foundation with grant num-ber PHY-2010970. T.E. acknowledges support by theVetenskapsr˚adet (Swedish Research Council) throughcontract No. 638-2013-8993 and the Oskar Klein Centrefor Cosmoparticle Physics.
Software:
Python/C language (Oliphant 2007),NumPy (van der Walt et al. 2011), Cython (Behnel et al.2011), SciPy (Jones et al. 2001–), MPI (Forum 1994), MPIfor Python (Dalc´ın et al. 2008), Matplotlib (Hunter 2007;Droettboom et al. 2018), IPython (Perez & Granger 2007),Jupyter (Kluyver et al. 2016),
MultiNest (Feroz et al.2009),
PyMultiNest (Buchner et al. 2014), GetDist(Lewis 2019).REFERENCES
Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016,PhRvL, 116, 061102—. 2017a, PhRvL, 119, 161101—. 2017b, ApJL, 848, L12—. 2017c, ApJL, 848, L13—. 2017d, Nature, 551, 85—. 2018, PhRvL, 121, 161101—. 2019, Physical Review X, 9, 031040—. 2020a, ApJL, 892, L3—. 2020b, Classical and Quantum Gravity, 37, 055002Abbott, R., Abbott, T. D., Abraham, S., et al. 2020c, arXive-prints, arXiv:2010.14527 —. 2020d, ApJL, 896, L44Ackley, K., Amati, L., Barbieri, C., et al. 2020, A&A, 643,A113Anand, S., Coughlin, M. W., Kasliwal, M. M., et al. 2020,Nature Astronomy, arXiv:2009.07210Andreoni, I., Goldstein, D. A., Kasliwal, M. M., et al. 2020,ApJ, 890, 131Antier, S., Agayeva, S., Almualla, M., et al. 2020, MNRAS,497, 5518Baiotti, L., & Rezzolla, L. 2017, Rep. Prog. Phys., 80,096901 Raaijmakers et al.