AT2017gfo: Bayesian inference and model selection of multi-component kilonovae and constraints on the neutron star equation of state
Matteo Breschi, Albino Perego, Sebastiano Bernuzzi, Walter Del Pozzo, Vsevolod Nedora, David Radice, Diego Vescovi
MMNRAS , 000–000 (0000) Preprint 8 January 2021 Compiled using MNRAS L A TEX style file v3.0
AT2017gfo: Bayesian inference and model selection of multi-componentkilonovae and constraints on the neutron star equation of state
Matteo Breschi , Albino Perego , , Sebastiano Bernuzzi , Walter Del Pozzo , ,Vsevolod Nedora , David Radice , , , Diego Vescovi , , Theoretisch-Physikalisches Institut, Friedrich-Schiller-Universit¨at Jena, Fr¨obelstieg 1, 07743, Jena, Germany Dipartimento di Fisica, Universit´a di Trento, Via Sommarive 14, 38123, Trento, Italy INFN-TIFPA, Trento Institute for Fundamental Physics and Applications, via Sommarive 14, 38123, Trento, Italy Dipartimento di Fisica “Enrico Fermi”, Universit´a di Pisa, Largo B. Pontecorvo 14, 56127, Pisa, Italy INFN, Sezione di Pisa, Largo B. Pontecorvo 14, 56127, Pisa, Italy Institute for Gravitation & the Cosmos, The Pennsylvania State University, University Park, PA 16802, USA Department of Physics, The Pennsylvania State University, University Park, PA 16802, USA Department of Astronomy & Astrophysics, The Pennsylvania State University, University Park, PA 16802, USA Gran Sasso Science Institute, Viale F. Crispi 7, 67100, L’Aquila, Italy INFN, Sezione di Perugia, Via A. Pascoli 23, 06123, Perugia, Italy INAF, Observatory of Abruzzo, Via M. Maggini, 64100, Teramo, Italy
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
The joint detection of the gravitational wave GW170817, of the short γ -ray burst GRB170817A and of the kilonova AT2017gfo,generated by the the binary neutron star merger observed on August 17, 2017, is a milestone in multimessenger astronomyand provides new constraints on the neutron star equation of state. We perform Bayesian inference and model selection onAT2017gfo using semi-analytical, multi-components models that also account for non-spherical ejecta. Observational data favoranisotropic geometries to spherically symmetric profiles, with a log-Bayes’ factor of ∼ , and favor multi-component modelsagainst single-component ones. The best fitting model is an anisotropic three-component composed of dynamical ejecta plusneutrino and viscous winds. Using the dynamical ejecta parameters inferred from the best-fitting model and numerical-relativityrelations connecting the ejecta properties to the binary properties, we constrain the binary mass ratio to q < . and thereduced tidal parameter to < ˜Λ < . Finally, we combine the predictions from AT2017gfo with those from GW170817,constraining the radius of a neutron star of . (cid:12) to . ± . ( σ level). This prediction could be further strengthenedby improving kilonova models with numerical-relativity information. Key words: transients: neutron star mergers – methods: data analysis– equation of state
On August 17, 2017, the ground-based interferometers of LIGO andVirgo (Abbott et al. 2018a; Aasi et al. 2015; Acernese et al. 2015)detected the first gravitational-wave (GW) signal coming from abinary neutron star (BNS) merger, known as GW170817 (Abbottet al. 2017a). GW170817 was followed by a short gamma-ray burst(GRB) GRB170817A (Abbott et al. 2017c; Savchenko et al. 2017),which reached the space observatories Fermi (Ajello et al. 2016) andINTEGRAL (Winkler et al. 2011) ∼ . after coalescence time.Eleven hours later, several telescopes started to collect photometricand spectroscopical data from AT2017gfo, an unprecedented elec-tromagnetic (EM) kilonova transient (Coulter et al. 2017; Chornocket al. 2017; Nicholl et al. 2017; Cowperthwaite et al. 2017; Pianet al. 2017; Smartt et al. 2017; Tanvir et al. 2017; Tanaka et al.2017; Valenti et al. 2017) coming from a coincident region of thesky. Kilonovae (kNe) are quasi-thermal EM emissions interpreted asdistinctive signature of r -process nucleosynthesis in the neutron-richmatter ejected from the merger and from the subsequent BNS rem- nant evolution (Smartt et al. 2017; Kasen et al. 2017; Rosswog et al.2018; Metzger 2020; Kawaguchi et al. 2019). The follow up of thesource lasted for more than a month and included also non-thermalemission from the GRB170817A afterglow (e.g., Nynka et al. 2018;Hajela et al. 2019).The combined observation of GW170817, GRB170817A andAT2017gfo decreed the dawn of multimessenger astronomy withcompact binaries (Abbott et al. 2017b). From these multimessen-ger observations it is possible to infer unique information on theunknown equation of state (EOS) of neutron star (NS) matter, (e.g.Radice et al. 2017; Margalit & Metzger 2017; Bauswein et al. 2017;Radice et al. 2018b; Dietrich et al. 2018). Indeed, the EOS deter-mines the tidal polarizability parameters that describe tidal inter-actions during the inspiral-merger and characterize the GW sig-nal (Damour et al. 2012; Bernuzzi et al. 2014). It also determinesthe outcome of BNS mergers (e.g. Shibata et al. 2005; Bernuzziet al. 2015a, 2020) and the subsequent postmerger gravitational sig-nal from the remnant (e.g. Bauswein et al. 2014; Bernuzzi et al.2015b; Zappa et al. 2019; Agathos et al. 2020; Breschi et al. 2019). © 0000 The Authors a r X i v : . [ a s t r o - ph . H E ] J a n M. Breschi et al.
At the same time, the amount of mass, the velocity, and the composi-tion of the ejecta are also strongly dependent on the EOS, that has animprint on the kN signature, e.g. (Hotokezaka et al. 2013; Bausweinet al. 2013; Radice et al. 2018d,a).The spectrum of AT2017gfo was recorded from ultraviolet (UV)to near infrared (NIR) frequencies (e.g., Pian et al. 2017; Nakar et al.2018), and the observations showed several characteristic features.At early stages, the kN was very bright and its spectrum peakedin the blue band 1 day after the merger (blue kN). After that, thepeak of the spectrum moved towards larger wavelengths, peaking atNIR frequencies between five to seven days after merger (red kN).Minimal models that can explain these features require more thanone component. In particular, minimal fitting models assume spher-ical symmetry and include a lathanide-rich ejecta responsible for thered kN, typically interpreted as dynamical ejecta, and another ejectawith material partially reprocessed by weak interaction, responsiblefor the blue component (e.g., Villar et al. 2017b). Numerical rela-tivity (NR) simulations show that the geometry profiles of the ejectaare not always spherically symmetric and their distributions are nothomogeneous (Perego et al. 2017a). Moreover, NR simulations alsoindicate the presence of multiple ejecta components, from the dy-namical to the disk winds ejecta (Rosswog et al. 2014; Fern´andezet al. 2015; Metzger & Fern´andez 2014; Perego et al. 2014; Nedoraet al. 2019). Therefore, this information has to be taken into accountduring the inference of the kN properties.The modeling of kNe is a challenging problem, due to the com-plexity of the underlying physics, which is affected by a diverse in-teractions and scales (see Metzger 2020, and references therein).Together with the choice of ejecta profiles, the lack of a reliabledescription of the radiation transport is a relevant source of uncer-tainties in the modeling of kNe, due to the incomplete knowledgeon the thermalization processes (Korobkin et al. 2012; Barnes et al.2016) and on the energy-dependent photon opacities in r -processmatter (Tanaka et al. 2019; Even et al. 2020). Current kN mod-els often use either simplistic ejecta profiles or simplistic radiationschemes, (e.g., Grossman et al. 2014; Villar et al. 2017b; Coughlinet al. 2017; Perego et al. 2017a). Given the challenges and uncertain-ties associated to the theoretical prediction of kN features, Bayesianinference and model selection of the observational data can provideimportant insights on physical processes hidden in the kN signature.In this work, we explore model selection in geometrical and ejectaproperties using simplified light curve (LC) models, that nonethe-less capture the key features of the problem. The inference resultsare then employed to derive constraints on the neutron star EOS. InSec. 2, we describe the semi-analytical model and the ejecta compo-nents used in our analysis. In Sec. 3, we recall the Bayesian frame-work for model selection, highlighting the choices of the relevantstatistical quantities, such as likelihood function and prior distribu-tions. In Sec. 4, we discuss the inference on AT2017gfo, criticallyexamining the posterior samples in light of targeted NR simula-tions (Perego et al. 2019; Nedora et al. 2019; Endrizzi et al. 2020;Nedora et al. 2020a; Bernuzzi et al. 2020) and previous analyses. InSec. 5, we discuss new constraints on the NS EOS focusing first onmass ratio and reduced tidal parameter for the source of GW170817,and then on the neutron star radius R . . We conclude in Sec. 6. In this section, we first summarize basic analytical results and scal-ing relations that characterize the kN emission, and then describe in detail the models we employ for the ejecta components and LCcalculations.
Let us consider a shell of ejected matter characterized by a massdensity (cid:37) , with total mass m and gray opacity κ (mean cross sectionper unit mass). The shell is in homologous expansion symmetricallywith respect to the equatorial plane at velocity v , such that its meanradius is R ∼ vt after a time t following the merger. Matter opac-ity to EM radiation can be expressed in terms of the optical depth, τ , which is estimated as τ (cid:39) (cid:37)κR . After the BNS collision, whenmatter becomes unbound and r -process nucleosynthesis occurs, theejecta are extremely hot, T ∼ K (e.g. de Jes´us Mendoza-Temiset al. 2015; Wu et al. 2016; Perego et al. 2019). However, at earlytimes the thermal energy is not dissipated efficiently since the envi-ronment is optically thick ( τ (cid:29) ) and photons diffuse out only onthe diffusion timescale until they reach the photosphere ( τ = 2 / ).As the outflow expands, its density drops ( ρ ∝ t − ) and the opticaldepth decreases.The key concept behind kNe is that photons can contribute to theEM emission at a given time t if they diffuse on a timescale compa-rable to the expansion timescales, i.e., if they escape from the shellsoutside R diff , where R diff is the radius at which the diffusion time t diff (cid:39) Rτ /c equals the dynamical time t (Piran et al. 2013; Gross-man et al. 2014; Metzger 2020) . In the previous expression, c is thespeed of light. Since t diff ∝ t − , a larger and larger portion of theejecta becomes transparent with time. The luminosity peak of thekN occurs when the bulk of matter that composes the shell becomestransparent. As first approximation, the characteristic timescale atwhich the light curve peaks is commonly estimated (Arnett 1982)as: t peak = (cid:114) mκ πβvc , (1)where the dimensionless factor β depends on the density profileof the ejecta. For a spherical symmetric, homologously expandingejecta ( β (cid:39) ) with mass m = 10 − M (cid:12) , velocity v = 0 . c and opacity in the range κ (cid:39) −
50 cm g − , which are typicalvalues respectively for lanthanide-free and for lanthanide-rich mat-ter (Roberts et al. 2011; Kasen et al. 2013), Eq. (1) predicts a char-acteristic t peak in the range –
10 days (Abbott et al. 2017d).In the absence of a heat source, matter would simply cool downthrough adiabatic expansion. However, the ejected material is con-tinuously heated by the radioactive decays of the r -process yields,which provide a time dependent heating rate of nuclear origin. Anadditional time dependence is introduced by the thermalization effi-ciency, i.e. the efficiency at which this nuclear energy, released in theform of supra-thermal particles (electrons, daughter nuclei, photonsand neutrinos), thermalizes within the expanding ejecta (see, e.g.,Metzger & Berger 2012; Korobkin et al. 2012; Barnes et al. 2016;Hotokezaka et al. 2018). The kN LCs in our work are computed using the multicompo-nent, anisotropic semi-analytical
MKN model first introduced inRef. (Perego et al. 2017a) and largely based on the kN models pre-sented in Refs. (Grossman et al. 2014) and (Martin et al. 2015). Theejecta are either spherical or axisymmetric with respect to the rota-tional axis of the remnant, and symmetric with respect to the equa-
MNRAS000
MNRAS000 , 000–000 (0000)
T2017gfo: model selection of multi-component kNe torial plane. The viewing angle ι is measured as the angle betweenthe rotational axis and the line of sight of the observer.For each component the ejected material is described through theangular distribution of its ejected mass, m , root-mean-square (rms)radial velocity, v rms , and opacity, κ . In axisymmetric models, thelatter quantities are functions of the polar angle θ , measured fromthe rotational axis and discretized in N θ = 30 angular bins evenlyspaced in cos θ . Additionally, within each ray, matter is radially dis-tributed with a stationary profile in velocity space, ξ ( v ) such that ξ ( v ) ∝ (1 − ( v/v max ) ) , where ξ ( v )d v is the matter contained inan infinitesimal layer of speed [ v, v + d v ] , and v max = v max ( v rms ) is the maximum velocity at the outermost edge of the component.The characteristic quantities (cid:37) , v and κ are then evaluated for ev-ery bin according to the assumed input profiles. For every bin, weestimate the emitted luminosity using the radial model describe inRef. (Perego et al. 2017a), in §4 of Ref. (Barbieri et al. 2020). Inparticular, the model assumes that the luminosity is emitted as ther-mal radiation from the photosphere (of radial coordinate R ph ), andthe luminosity and the photospheric surface determine the effectiveemission temperature, T eff through the Stefan-Boltzmann law. Weexpect this assumption to be well verified at early times (with a fewdays after merger), while deviations from it are expected to becomemore and more relevant for increasing time.The time-dependent nuclear heating rate (cid:15) nuc entering these cal-culations is approximated by an analytic fitting formula, derivedfrom detailed nucleosynthesis calculations (Korobkin et al. 2012), (cid:15) nuc ( t ) = (cid:15) (cid:15) th ( t )0 . (cid:15) nr ( t ) (cid:20) − π arctan (cid:18) t − t σ (cid:19)(cid:21) α , (2)where σ = 0 .
11 s , t = 1 . , α = 1 . and (cid:15) th ( t ) is the thermal-ization efficiency tabulated according to Ref. (Barnes et al. 2016).The heating factor (cid:15) nr ( t ) is introduced as in Ref. (Perego et al.2017a) to roughly improve the behavior of Eq. (2) in the regimeof mildly neutron-rich matter (characterized by an initial electronfraction Y e (cid:38) . ), (see, e.g. Martin et al. 2015): (cid:15) nr ( t, κ ) = [1 − w ( κ )] + w ( κ ) (cid:15) Y e ( t ) , (3)where w ( κ ) is a logarithmic smooth clump function such that w ( κ < g − ) = 1 and w ( κ >
10 cm g − ) = 0 and thefactor (cid:15) Y e ( t ) encodes the dependence on Y e : if Y e < . , then (cid:15) Y e ( t ) = 1 , otherwise, when Y e ≥ . , (cid:15) Y e ( t ) = (cid:15) min + (cid:15) max (cid:104) e t/t (cid:15) − (cid:105) − , (4)where t (cid:15) = 1 day , (cid:15) min = 0 . and (cid:15) max = 2 . .Furthermore, in order to improve the description in the high-frequency bands (i.e., V , U , B and g ) within the timescale of thekilonova emission, and following Ref. (Villar et al. 2017a), we in-troduce a floor temperature, i.e. a minimum value for T eff . This isphysically related to the drop in opacity due to the full recombina-tion of the free electrons occurring when for the matter temperaturedrops below T floor (Kasen et al. 2017; Kasen & Barnes 2018). Un-der these assumptions, the condition T eff = T floor becomes a goodtracker for the photosphere location. Since kNe are powered by theradioactive decay of different blends of atomic species, we introducein our model two floor temperatures, T Nifloor and T LAfloor , that character-ize respectively the recombination temperature of lanthanides-freeand of lanthanide-rich ejecta.Eventually, the emissions coming from the different rays are com-bined to obtain the spectral flux at the observer location: F ν ( n , t ) = (cid:90) n Ω · n > (cid:18) R ph (Ω , t ) D L (cid:19) B ν ( T eff (Ω , t )) n · d Ω (5) Remnant
Isotropic ejecta Anisotropic ejecta
Orbital plane Dynamical component (D) Viscous component (V) ν -driven wind (N) Figure 1.
Graphic representation of the analyzed ejecta profiles for isotropicand anisotropic cases from an azimuthal perspective and for a fixed momentof time. The black dot represents the remnant and the dashed line is the pro-jected orbital plane of the binary. The shadowed areas describe the ejectaprofiles: the shape characterizes the mass distribution, while the colors referto the prior assumptions on the opacity parameter. In particular, blue regionsdenote opacities lower than g − , red regions refer to opacities greaterthan g − , and oranges areas indicate a broadly distributed opacity.All shells are isotropically expanding with a constant velocity. where n is the unitary vector along the line of sight, n Ω is the uni-tary vector spanning the solid angle Ω , D L is the luminosity dis-tance, R ph is the local radial coordinate of the photospheric surface,and B ν ( T eff ) is the spectral radiance at frequency ν for a surface oftemperature T eff . Lastly, from Eq. (5), it is possible to compute theapparent AB magnitude m b in a given photometric band b as: mag b ( n , t ) = − . ( F ν b ( n , t )) − . , (6)where ν b is the effective central frequency of band b . In order to describe the different properties of AT2017gfo it is nec-essary to appeal to a multi-component structure for the ejecta pro-ducing the kN. Different components are characterized by differentsets of intrinsic parameters, m , v rms and κ , and by their angular dis-tributions with respect to θ .Given the angular profiles of the characteristic parameters, thephysical luminosity produced by each component inside a ray iscomputed by using the model outlined in the previous section. Then,the total bolometric luminosity of the ray is given by the sum of thesingle contributions, i.e. L ( t ) = (cid:80) k L ( k ) ( t ) where k runs over thecomponents. The outermost photosphere is the one that determinesthe thermal spectrum of the emission Once R ph and T eff have beendetermined, the spectral flux and the AB magnitudes are computedaccording to Eqs. (5) and (6).We perform the analysis using two different assumptions on theprofiles of the source. Initially, we impose completely isotropic pro-files for every parameter of every ejecta component. These casesare labeled as isotropic, ‘ISO’. Subsequently, we introduce angularprofiles as functions of the polar angle for the mass and opacity pa-rameters, while we keep v rms always isotropic. This second case is MNRAS , 000–000 (0000)
M. Breschi et al. labeled as anisotropic, ‘ANI’. In parallel, we explore models with adifferent number of components. We always assume the presence ofthe dynamical ejecta, while we add to them one or two qualitativelydifferent disk-wind ejecta components.In the following paragraphs, we describe the physical assumptionson each component and the choice of the prior distributions (seeTab. 1). Fig. 1 shows a graphical representation of the employedejecta components.
Dynamical ejecta (D) . The BNS collision ejects unbound matteron the dynamical timescale, whose properties strongly depend onthe total mass of the BNS, on the mass ratio and on the EOS (e.g.Hotokezaka et al. 2013; Rosswog et al. 2013; Bauswein et al. 2013;Radice et al. 2016; Bovard et al. 2017; Radice et al. 2018c,d). Thisejection is due to tidal torques and shocks developing at the con-tact interface between the merging stars, when matter is squeezedout by hydrodynamical processes (Oechslin et al. 2006; Hotokezakaet al. 2013). The expansion of this ejecta component has a velocityof roughly v rms ∼ . c . Moreover, this phenomenon generates adistribution of ejected mass denser in the regions across the orbitalplane with respect to the region along its orthogonal axis, charac-terized by larger opacities at lower latitudes. In particular, neutrinoirradiation (if significant), increases the ejecta Y e and prevents theformation of lanthanides. For the anisotropic analyses, the mass pro-file is taken to be (cid:37) ( θ ) ∝ sin θ , and the opacity profile is take asa step function in the polar angle characterized by the parameters ( κ low , κ high ) , respectively for low- and high-latitudes, with a stepangle θ step = π/ (see Sec. 3.3). In terms of emitted LC, the de-scribed ejecta is characterized by a red equatorial component and ablue contribution at higher latitudes. Neutrino-driven wind (N) . Simulations of the remnant evolutionin the aftermath of a BNS merger reveal the presence of other ejec-tion mechanisms happening over the thermal and viscous evolutiontimescales (e.g. Metzger et al. 2008; Fern´andez & Metzger 2013;Perego et al. 2014, 2017b; Decoene et al. 2020). If the ejection hap-pens while the remnant is still a relevant source of neutrinos, neu-trino irradiation has enough time to increase Y e above 0.25, pre-venting full r -process nucleosynthesis, especially close to the polaraxis. Detailed simulations (Perego et al. 2014; Martin et al. 2015;Fujibayashi et al. 2018, 2020) show that a relatively small fractionof the expelled disk contributes to this component and its velocity isexpected to be v rms (cid:46) . c . For anisotropic analyses, the mass pro-file is taken to be uniform in the range θ ∈ [0 , π/ and negligibleotherwise, while the opacity profile is takes as as step function in thepolar angle, with a step angle θ step = π/ . Disk’s viscous ejecta (V) . In addition to neutrinos, viscous torquesof dynamical and magnetic origin can unbind matter from the diskaround massive NSs or black holes (Metzger et al. 2010; Metzger& Fern´andez 2014; Just et al. 2015). This viscous component isexpected to unbound a large fraction of the disk matter on longertimescale, reaching m (cid:46) − M (cid:12) , with a relatively low velocity, v rms (cid:46) . c . The corresponding ejecta are more uniformly dis-tributed over the polar angle than the dynamical ejecta and the ν -driven wind ejecta. The presence or the lack of a massive NS in thecenter can influence the Y e of these ejecta. Then, all angular pro-files are assumed to be isotropic for this component (Wu et al. 2016;Siegel & Metzger 2018).We conclude this section by recalling that the main motivation be-hind the usage of the semi-analytic model presented above is the op-timal compromise between its robustness and adaptability, essential to model the non-trivial structure of the ejecta, and the reduced com-putational costs, necessary to perform parameter estimation studies.However, it has been showed that simplified models that avoid thesolution of the radiation transport problem can suffer from system-atic uncertainties (Wollaeger et al. 2018). In particular, the analyticalmodel presented in (Grossman et al. 2014), on which ours is based,produces significantly lower light curves. The comparison with ob-served kN light curves and more detailed kN models showed howlarger nuclear heating rates (cid:15) systematically reduce this discrep-ancy. In this section, we recall the basic concepts of model selection asthey are stated in the Bayesian theory of probability. Then, we de-scribe the statistical technique used for the computations of theBayes’ factors. As convention, the symbol ‘ log ’ denotes the naturallogarithm while a logarithm to a different base is explicitly writtenwhen it is used.
Given some data d and a model H (hypothesis) described by a set ofparameters θ , the posterior probability is given by the Bayes’ theo-rem: p ( θ | d, H ) = p ( d | θ , H ) p ( θ | H ) p ( d | H ) , (7)where p ( d | θ , H ) is the likelihood function, p ( θ | H ) is the priorprobability assigned to the parameters and p ( d | H ) is the evidence.The latter value plays the role of normalization constant and it canbe computed by marginalizing the likelihood function, p ( d | H ) = (cid:90) Θ p ( d | θ , H ) p ( θ | H )d θ , (8)where the integral is computed over the entire parameters’ space Θ .In the framework of Bayesian theory of probability, we can com-pare two models, say A and B , by computing the ratios of the re-spective posterior probabilities, also known as Bayes’ factor , B AB = p ( A | d, H A ) p ( B | d, H B ) . (9)Using Eq. (7) and assuming that the data does not depend on thedifferent hypothesis, we get: B AB = p ( d | A, H A ) p ( d | B, H B ) . (10)Now suppose that the two models A, B are respectively describedby two sets of parameters θ A , θ B . Using the marginalization rulewe can write: p ( d | I, H I ) = (cid:90) Θ I p ( d | θ I , I, H I ) p ( θ I | I, H I ) d θ I , (11)for I = A, B . The integral in Eq. (11) represents the evidence com-puted for the hypotheses H (cid:48) I = { H I , I } , for I = A, B (i.e. theinvolved model becomes part of the background hypothesis). Then,we obtain that the Bayes’ factor B AB can be computed as B AB = p ( d | H (cid:48) A ) p ( d | H (cid:48) B ) . (12)From the previous results, we understand that if B AB > thenthe model A will be favored by the data, viceversa if B AB < . It MNRAS000
Given some data d and a model H (hypothesis) described by a set ofparameters θ , the posterior probability is given by the Bayes’ theo-rem: p ( θ | d, H ) = p ( d | θ , H ) p ( θ | H ) p ( d | H ) , (7)where p ( d | θ , H ) is the likelihood function, p ( θ | H ) is the priorprobability assigned to the parameters and p ( d | H ) is the evidence.The latter value plays the role of normalization constant and it canbe computed by marginalizing the likelihood function, p ( d | H ) = (cid:90) Θ p ( d | θ , H ) p ( θ | H )d θ , (8)where the integral is computed over the entire parameters’ space Θ .In the framework of Bayesian theory of probability, we can com-pare two models, say A and B , by computing the ratios of the re-spective posterior probabilities, also known as Bayes’ factor , B AB = p ( A | d, H A ) p ( B | d, H B ) . (9)Using Eq. (7) and assuming that the data does not depend on thedifferent hypothesis, we get: B AB = p ( d | A, H A ) p ( d | B, H B ) . (10)Now suppose that the two models A, B are respectively describedby two sets of parameters θ A , θ B . Using the marginalization rulewe can write: p ( d | I, H I ) = (cid:90) Θ I p ( d | θ I , I, H I ) p ( θ I | I, H I ) d θ I , (11)for I = A, B . The integral in Eq. (11) represents the evidence com-puted for the hypotheses H (cid:48) I = { H I , I } , for I = A, B (i.e. theinvolved model becomes part of the background hypothesis). Then,we obtain that the Bayes’ factor B AB can be computed as B AB = p ( d | H (cid:48) A ) p ( d | H (cid:48) B ) . (12)From the previous results, we understand that if B AB > thenthe model A will be favored by the data, viceversa if B AB < . It MNRAS000 , 000–000 (0000)
T2017gfo: model selection of multi-component kNe is important to observe that the Bayes’ factor implicitly takes intoaccount the so called Occam’s razor, i.e. if two models are both ableto capture the features of the data, then the one with lower number ofparameters will be favored (Sivia & Skilling 2006). In our analysis,this is a crucial point since different models have different numbersof parameters. In a realistic scenario, the form of the likelihood function is ana-lytically indeterminable and the parameter space has a non-trivialnumber of dimensions. For these reasons, the estimation of Eq. (11)is performed resorting to statistical computational techniques: weemploy the nested sampling Bayesian technique introduced inRef. (Skilling 2006) and designed to compute the evidence and ex-plore the full parameter space. The uncertainties associated with theevidence estimations are computed according to Ref. (Skilling 2006)and increasing the result by one order of magnitude, in order to con-servatively take into account systematics. The latter are expectedsince the model considered for our analyses (as many others) can-not capture all the physics processes involved in kNe, and it suffersof large uncertainties in the atomic physics and radiative processesimplementation.We perform inference with cpnest (Pozzo & Veitch Pozzo& Veitch), a parallelized nested sampling implementation. We use1024 live points and, for every step, we set a maximum numberof 2048 Markov-chain Monte Carlo (MCMC) iterations for the ex-ploration of the parameter space. The proposal step method usedin the MCMC is the same as the one implemented as default in cpnest software. It corresponds to a cycle over four different pro-posal methods: a random-walk step (Goodman & Weare 2010), astretch move (Goodman & Weare 2010), a differential evolutionmethod (Nelson et al. 2013) and a proposal based on the eigenvec-tors of the covariance matrix of the ensemble samples, as imple-mented in Ref. (Veitch et al. 2015).
In our analysis we assume the sky position of the source to beknown and the time of coalesce to be the same of the trigger time ofGW170817 (Abbott et al. 2017a). Furthermore, we do not take intoaccount the redshift contribution, given the larger systematic uncer-tainties in the model. We employ the parameters shown in Tab. 1,that can be divided in three subsets: the intrinsic ejecta parameters θ (D , V , N)ej , the intrinsic global parameters θ glob , and the extrinsic pa-rameters θ ext .The intrinsic ejecta parameters, θ ( k )ej for k = D , V , N , character-ize the properties of each ejecta component and they are the amountof ejected mass, m , the rms velocity of the fluid, v rms , and their greyopacity, κ . Under the assumption of isotropic geometry, the intrin-sic ejecta parameters θ ( k )ej are defined by a single value for everyshell, i.e. a single number characterizes the entire profile of the pa-rameter of interest, since it is spherically symmetric. However, foranisotropic cases, we have to introduce more than one independentparameters to describe an angular profile for a specific variable: thisis the case of the opacity parameter of the dynamical component,where the profile is chosen as step functions characterized by twodifferent parameters, κ low and κ high , respectively at low and highlatitudes. In such a cases, the angle θ step is introduced to denote theangle at which the profile changes value, as mentioned in Sec. 2.3.The intrinsic global parameters, θ glob , represent the properties of Table 1.
List of intrinsic and extrinsic parameters involved in the analysisand the respective prior bounds for the cases of anisotropic geometry. Forisotropic geometry cases, the bounds are identical except for the opacity κ ofdynamical component (D), where the low-latitude and high-latitude boundsare joined together.Intrinsic Ejecta Parameters θ (D , V , N)ej
Comp. m v rms κ high κ low θ step [10 − M (cid:12) ] [ c ] [cm g − ] [rad] D [0.1, 10] [0.15,0.333] [0.1,5] [5,30] π/ N [0.01,0.75] [0.05,0.15] [0.01,5] π/ V [1,20] [0.001,0.1] [0.01,30] –Intrinsic Global Parameter θ glob T Nifloor [K] [500, 8000] T LAfloor [K] [500, 8000] (cid:15) [ erg g − s − ] [2 × , × ] Extrinsic Parameters θ ext D L [Mpc] [15,50] ι [deg] [0,70] the source common to every component, such as the floor tempera-tures, T Nifloor and T LAfloor , and the heating rate constant (cid:15) . In principle,the latter is a universal property which defines the nuclear heatingrate as expressed in Eq. 2. The whole set of intrinsic parameters, θ glob and θ ( k )ej , determines the physical dynamics of the system and,therefore, they determine the properties of the kN emission, irrespec-tively of the observer location.The extrinsic parameters, θ ext , are the luminosity distance of thesource, D L and the viewing angle ι . These parameters do not dependon the physical properties of the source and they are related with theobserved signal through geometrical argumentation.The prior distributions for all the parameters are taken uniformin their bounds, except for the followings. For the extrinsic param-eters θ ext = { D L , ι } , we set the priors equal to the marginalizedposterior distributions coming from the low-spin-prior measurementof GW170817 (Abbott et al. 2019b); For the heating rate factor (cid:15) ,we use a uniform prior distribution in log (cid:15) , i.e. p ( (cid:15) | H ) ∝ (cid:15) − ,since this parameter strongly affects the LC and it is free to vary ina wide range. Moreover, we adopt a prior range according with theestimation given in Ref. (Korobkin et al. 2012). Tab. 1 shows theprior bounds used for the analysis of the anisotropic cases. For theisotropic studies, the bounds are identical except for the opacity κ of dynamical component, where the low-latitude and high-latitudebounds are joined together. The data { d b,i ± σ b,i } are the apparent magnitudes observed fromAT2017gfo, with their standard deviations. They have been collectedfrom (Villar et al. 2017b), where all the precise reference to the orig-inal works and to the data reduction techniques can be found. Theindex b runs over all considered photometric bands, covering a widephotometric range from the UV to the NIR, while for each band b MNRAS , 000–000 (0000)
M. Breschi et al. the index i runs over the corresponding sequence of N b temporalobservations. Additionally, the magnitudes have been corrected forGalactic extinction (Cardelli et al. 1989). We introduce a Gaussianlikelihood function in the apparent magnitudes with mean and vari-ance, d b,i , σ b,i , from the observations of AT2017gfo, log p ( θ | d, H ) ∝ − (cid:88) b N b (cid:88) i =1 (cid:12)(cid:12) d b,i − mag b,i ( θ ) (cid:12)(cid:12) σ b,i , (13)where mag b,i ( θ ) are the magnitudes generated by the LC model,of Sec. 2, which encodes the dependency on the parameters θ , forevery band b at different times i . The likelihood definition Eq. (13)is in accordance with the residuals introduced in Ref. (Perego et al.2017a) and it takes into account the uncertainties due to possibletechnical issues of the instruments and generic non-stationary con-tributions, providing a good characterization of the noise . For bothgeometric configurations, isotropic (ISO) and anisotropic (ANI), weperform Bayesian analyses using different combinations of compo-nents, testing the capability to fit the data. In this section, we present the results gathered from the Bayesiananalysis. In Sec. 4.1 we describe the capability of the synthetic LCsto fit the observed data. After that, in Sec. 4.2, we discuss the es-timated evidence inferring the preferred model. Finally, in Sec. 4.3,we discuss the interpretation of the recovered posterior distributions.
Figure 2 shows the LCs computed from the recovered maximum-likelihood parameters for each discussed models. The estimated LCsare compared with AT2017gfo data for six representative photomet-ric bands. Moreover, Fig. 3 shows the uncertainties associated withthe estimated LCs, computed over the recovered posterior samples,for each considered model. Generally, the errors associated withthe near UV (NUV) magnitudes are larger compared with the otherbands, reflecting the lower number of data points in this photometricregion. Furthermore, none of the considered model is able to fullycapture the trend described by the observed data in the Ks band fortime larger then 10 days, within the provided prior bounds. This isexpected from the simplified treatment of the radiation transport andthe approximated heating rate in our models.The isotropic models (ISO-D and ISO-DV) give a good fittingto the data for early times and their LCs capture the general trendsof the data. However, for times larger than ∼ days, these modelsdo not capture all the features of the data within the provided priorbounds. This inaccuracy is particularly evident in the NIR, where theLCs predicted by the ISO-D and the ISO-DV models do not recoverthe correct slopes of the data.The anisotropic single-component case, ANI-D, is apt at adapt-ing the model to the different features present in the data, evenfor large time-scales. However, it overestimates the kN emissionin the blue band. This inconsistency could be reduced allowing thehigh latitude opacity parameter κ high to lower values. Regarding the Also the work presented in Ref. (Villar et al. 2017b) employs a Gaussianlikelihood, with the inclusion of an additional uncertainty parameter; while,in Ref. (Coughlin et al. 2017), the authors proposed a likelihood distributedas a χ . Table 2.
Estimated log-evidences for the analyzed kNe models. The re-ported uncertainties refers to the standard deviations estimated according toRef. (Skilling 2006).Profile Components log p ( d | Model)
ISO D − ± ISO D+V − ± ANI D − ± ANI N+V − ± ANI D+V − ± ANI D+N+V − ± anisotropic two-components models, the ANI-VN gives a good fit-ting for early times, but the model largely underestimates the data attimes (cid:38) days. This is due to the absence of a fast blue component.The anisotropic ANI-DV model gives LCs similar to ANI-D exceptfor a slight excess of power for time (cid:38) days, especially in the NIRregion, i.e. z , K and K s bands. This behavior could be mitigatedby reducing the lower bound on the T LAfloor parameter. However, itcould also indicate a significant deviation from the black-body emis-sion adopted in our model at late times. Furthermore, the ANI-DVmodel overshoots the data in the NUV, as it is for the respectivesingle-component case ANI-D. This can be explained looking at therecovered value of dynamical ejected mass, which exceeds theoreti-cal expectations estimated from NR simulations (Perego et al. 2019;Endrizzi et al. 2020; Nedora et al. 2019; Bernuzzi et al. 2020; Ne-dora et al. 2020a)(see Sec. 4.3.4). Similar considerations hold forthe anisotropic three-component case ANI-DNV. However, the un-certainties on the estimated LCs for this model are narrower withrespect to the ones obtained from the ANI-DV, corresponding toan improvement in the capability of constraining the measurement.The main improvement of the three-component ANI-DVN modelover the two-component ANI-DV model lies in its ability to betterfit early-times data due to the inclusion of a third component.
The logarithmic evidences estimated for the considered models areshown in Tab. 2. The evidence increases with the number of models’components. This is consistent with the hierarchy observed in theLC residuals, and the better match to the data for multi-componentmodels. The only exception is the ANI-NV case, for which the fea-tures of the data at late times are not well captured due to the absenceof a fast equatorial component. Furthermore, for a fixed number ofcomponents, the anisotropic geometries are always favored with re-spect to isotropic geometries, with a log B ANIISO of the order of .The preferred model among the considered cases is the anisotropicthree-component, in agreement with previous findings, e.g. (Cow-perthwaite et al. 2017; Perego et al. 2017a; Villar et al. 2017b). In the following paragraphs, we discuss the properties of the pos-terior distributions for each model and their physical interpretation.Table 3 and Tab. 4 show the mean values of the parameters, andtheir 90% credible regions, extracted from the recovered posteriordistributions. A general fact is that the marginalized posterior for the
MNRAS000
MNRAS000 , 000–000 (0000)
T2017gfo: model selection of multi-component kNe m ag B band
ISO-DISO-DVANI-DANI-DVANI-VNANI-DVN g band m ag r band AT2017gfo z band t [days] m ag i band t [days] Ks band
Figure 2.
Apparent magnitudes computed using the maximum-likelihood parameter for each considered model; ISO-D in blue, ISO-DV in yellow, ANI-D ingreen, ANI-DV in red, ANI-VN in purple and ANI-DVN in brown. The different panels refer to different photometric bands, respectively B , g , r , z , i and Ks .The black squares are the observed data of AT2017gfo for the corresponding photometric band with the respective standard deviations. ejected mass of the viscous component is always constrained againstthe lower bound − M (cid:12) , when this component is involved. More-over, for the majority of the analyses, the distance parameter is bi-ased towards larger values, inconsistently with the estimates fromRef. (Abbott et al. 2017a,b), and the heating rate parameter (cid:15) isgenerally overestimated comparing with the estimates from nuclearcalculations (Korobkin et al. 2012; Barnes et al. 2016; Kasen &Barnes 2018; Barnes et al. 2020; Zhu et al. 2020). This behaviorcan be explained from Eqs. (2), (5) and (6): D L and (cid:15) are largely degenerate and both concur to determine the brightness of the ob-served LCs. Thus, the correlations between these parameters inducebiases in the recovered values. The physical explanation of this ef-fect can be motivated with the poor characterization of the modelin the NIR bands: this lack of knowledge generates a fainter kN inthis photometric region and, in order to match the observed data, therecovered heating rate are larger. Note that this bias concurs in theoverestimation of the LC in the high-frequency bands (i.e. U , B and MNRAS , 000–000 (0000)
M. Breschi et al. − . . . D e v i a t i o n B band
ISO-DISO-DVANI-DANI-DVANI-VNANI-DVN g band − . . . D e v i a t i o n r band z band t [days] − . . . D e v i a t i o n i band t [days] Ks band
Figure 3.
Deviations from the maximum-likelihood template of the LCs computed from the whole set of posterior samples. The solid lines represent the medianvalues and the shadowed areas are the 90% credible regions. Different color refers to a different model; respectively, ISO-D in blue, ISO-DV in yellow, ANI-Din green, ANI-DV in red, ANI-VN in purple and ANI-DVN in brown. The different panels show different photometric bands, respectively B , g , r , z , i and Ks . Table 3.
Recovered values from the posterior distributions of the of the intrinsic ejecta parameters. The reported quantities are the means with the 90% credibleregions. The conventions (cid:38) , (cid:46) denote marginalized posterior distributions constrained respectively around the upper and the lower prior bounds. We remarkthat κ low and κ high refer respectively to the gray opacity parameters for low and high latitudes. Model Dynamical ejecta Viscous ejecta ν -driven wind m v rms κ high κ low m v rms κ m v rms κ (cid:2) − M (cid:12) (cid:3) [ c ] (cid:2) cm g − (cid:3) (cid:2) − M (cid:12) (cid:3) [ c ] (cid:2) cm g − (cid:3) (cid:2) − M (cid:12) (cid:3) [ c ] (cid:2) cm g − (cid:3) ISO-D . +0 . − . . +0 . − . . +0 . − . – – – – – – –ISO-DV . +0 . − . . +0 . − . . +0 . − . (cid:46) (cid:38) . . +0 . − . – – –ANI-D . +0 . − . . +0 . − . (cid:46) . (cid:38) – – – – – – –ANI-DV . +0 . − . . +0 . − . (cid:46) . . +0 . − . (cid:46) . +0 . − . . +0 . − . – – –ANI-VN – – – – (cid:46) . +0 . − . . +0 . − . (cid:38) .
75 0 . +0 . − . . +0 . − . ANI-DVN . +0 . − . . +0 . − . (cid:46) . . +0 . − . (cid:46) . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . V ), where the number of measurements is lower with respect to theother employed bands. We start considering the simplest employed model, the isotropicone-component model labelled as ISO-D. Fig. 4 shows the marginal-ized posterior distribution in the ( m, v rms ) plane. The velocity isconstrained around ∼ . c while the ejected mass lies around × − M (cid:12) , both in agreement with the observational results re-covered in Ref. (Villar et al. 2017b; Cowperthwaite et al. 2017;Abbott et al. 2017d; Coughlin et al. 2018). Moreover, the opacity posterior peaks in proximity of κ ∼ g − , consistently withRef. (Cowperthwaite et al. 2017).Regarding the extrinsic parameters, the posterior for the inclina-tion angle ι is coincident with the imposed prior, since the employedprofiles do not depend on this coordinate. The model is not able toconstrain the value of T Nifloor , which returns a posterior identical tothe prior, while T LAfloor is recovered around . The obtained flatposterior distribution for the T Nifloor parameter highlights the unsuit-ability of this model in capturing the features of the observed data.
MNRAS000
MNRAS000 , 000–000 (0000)
T2017gfo: model selection of multi-component kNe Table 4.
Recovered values from the posterior distributions of the of the globalintrinsic parameters and of the extrinsic parameters. The reported quantitiesare the means with the 90% credible regions. The conventions (cid:38) , (cid:46) denotemarginalized posterior distributions constrained respectively around the up-per and the lower prior bounds. Model T Nifloor T LAfloor (cid:15) ι D L [K] [K] (cid:2) erg g − s − (cid:3) [deg] [Mpc] ISO-D +3157 − +450 − . +1 . − . +27 − (cid:38) ISO-DV +778 − +243 − . +0 . − . +24 − . +0 . − . ANI-D +47 − +219 − +3 − . +0 . − . (cid:38) ANI-DV +105 − +175 − . +0 . − . . +0 . − . (cid:38) ANI-VN +56 − (cid:46)
500 8 . +0 . − . +1 − . +0 . − . ANI-DVN +105 − (cid:46)
500 30 . +0 . − . +1 − (cid:38) For the anistropic single-component model ANI-D, the value of theejected mass agrees with the one coming from the ISO-D case.However, in order to fit the data, ANI-D requires a larger velocity, ∼ . c , as shown in Fig. 4. The high-latitude opacity is constrainedaround the lower bound . g − while the low-latitude con-tribution exceeds above
30 cm g − , that largely differs from therespective isotropic case, ISO-D. In practice, that is due to the lackof ejected mass that is balanced with a more opaque environment.Nevertheless, according to the estimated evidences, this model ispreferred with respect to the isotropic case. The reason is clear fromFig. 2: the anisotropic model is able to characterize the late-timesfeatures of the data.The posterior distribution for viewing angle ι peaks around 44 de-grees, inconsistently with the estimations coming from the GRBanalysis (Abbott et al. 2017c; Savchenko et al. 2017; Ghirlanda et al.2019). Moreover, unlike the ISO-D case, both temperature parame-ters T Nifloor and T LAfloor are well constrained for the ANI-D analysis:these parameters affect mostly the late-times model, modifying theslope of the recovered LCs. Thus, these terms are responsible of theimprovement in the fitted LCs.
For both components, the individual most-likely value for ejectedmass parameter lies around ∼ − M (cid:12) , in agreement with themeasurement presented in Ref. (Abbott et al. 2017d). This range ofvalues is slightly overestimating the expectations coming from NRsimulations for the dynamical component (Perego et al. 2019; Ne-dora et al. 2019; Endrizzi et al. 2020; Nedora et al. 2020a; Bernuzziet al. 2020). This could be explained by considering the effect ofthe spiral-wave wind (Nedora et al. 2019), that constitute a massiveand fast ejecta on timescales of − ms. The spiral-wave windis not considered as components in our models because it would behighly degenerate with the dynamical ejecta. The recovered opacityparameters are roughly − g − . The velocity of the dynamicalcomponent is greater than secular velocity, accordingly with the the-oretical expectations. Comparing with other fitting models, the re-covered ejected masses m (D) result smaller with respect to the anal-ogous analysis of Ref. (Villar et al. 2017b), while the results roughlyagree with the estimations coming from Ref. (Coughlin et al. 2018).However, it is not possible to perform an apple-to-apple comparisonbetween these results, due to the systematic differences in model-ing between the semi-analytical model (used in this work) and theradiative-transport methods employed in Ref. (Villar et al. 2017b;Coughlin et al. 2018). The temperature parameters, T Nifloor and T LAfloor , are much moreconstrained comparing with the respective isotropic single compo-nent case ISO-D, and this is reflected in the improvement of fittingthe different trends of the data in the high-frequency bands. Themarginalized posterior distribution of the inclination angle is coin-cident with the prior, according with the isotropic description. Fur-thermore, the biases on the distance D L and the heating parameter (cid:15) are reduced with respect to the ISO-D, since two-component caseaccounts for a larger amount of total ejected mass.According with the estimated evidences, the isotropic two-components ISO-DV model is disfavored with respect to theanisotropic single-component ANI-D. The main difficulty of ISO-DV is, again, to fit the data at late-times. The ANI-DV model is the second best fitting model to AT2017gfoamong the considered cases. Fig. 4 shows the posterior distributionfor the mass and the velocity parameters for the dynamical com-ponent. The ejected mass value lies around ∼ − M (cid:12) , in agree-ment with previous estimates (Abbott et al. 2017d). On the otherhand, the recovered mass slightly overestimates the results comingfrom targeted NR simulations (Perego et al. 2019; Nedora et al.2019; Endrizzi et al. 2020; Nedora et al. 2020a; Bernuzzi et al.2020), similarly to ISO-DV (see Sec. 4.3.3). The velocity is wellconstrained around ∼ . c . The recovered low-latitude opacity cor-responds roughly to
12 cm g − and high-latitude opacity is con-strained around the lower bound, . g − . This result can beexplained by considering that the mass of the dynamical componentslightly overshoots the NR expectations (Perego et al. 2019; Nedoraet al. 2019; Endrizzi et al. 2020; Nedora et al. 2020a; Bernuzzi et al.2020) (of a factor ∼ . ), and by noticing that the ejected masscorrelates with the luminosity distance and the heating factor (thatare generally biased). This combination generates the overestima-tion of the data in the NUV region. In order to improve the fitting tothe observed data, the model tries to compensate this effect and thehigh-latitude opacity tends to move towards lower values.Concerning the viscous component, its velocity results an orderof magnitude smaller than the one of the dynamical ejecta, in agree-ment with the expectations. This enforce the hypothesis for whichthe viscous ejecta contributes mostly to the red kN. The the pos-terior distribution of opacity parameter peaks around ∼ g − ,denoting a medium opaque environment.The temperatures T Nifloor and T LAfloor are well constrained respec-tively around ∼ and ∼
700 K . The posterior for inclinationangle results similar to the ANI-D case, according with the fact thatthe viscous component, as we have defined it, does not introducefurther information on the inclination.
According to Tab. 2, this ANI-VN is the least likely model amongall anisotropic cases. As previously mentioned, the reason for thisis clear from the LCs. The parameters of the viscous component arecharacterized by a slow velocity of ∼ × − c and a low opacityenvironment, κ ∼ . g − . On the other hand, the neutrino-driven wind mass is overestimated compared with aftermath compu-tations presented in Ref. (Perego et al. 2017a), in order to compen-sate the lack of overall ejected mass due to the absence of a dynami-cal component. Moreover, the neutrino-driven wind is characterizedby a realistic velocity of ∼ . c , and by a low-opaque environment, κ ∼ g − . MNRAS , 000–000 (0000) M. Breschi et al.
ISO − DANI − D .
76 0 .
78 0 .
80 0 .
82 0 . m (D) [10 − M (cid:12) ] . . . . . . . . v ( D ) [ c ] ANI − DVANI − DVN . . . . . m (D) [10 − M (cid:12) ] . . . . . v ( D ) [ c ] Figure 4.
Marginalized posterior distribution of ejected mass m and velocity v rms of dynamical component. Left panel: comparison for the one-componentstudies, ISO-D and ANI-D. The anisotropic case requires larger velocities in order to fit the observed data. Right panel: comparison for the two preferredmodels, ANI-DV and ANI-DVN. Regarding the extrinsic parameters, the ANI-VN model is the casethat gives the best agreement with Ref. (Korobkin et al. 2012) interms on heating factor. The distance, instead, is recovered around ∼
20 Mpc , underestimating the GW distance (Abbott et al. 2017a).The T Nifloor parameter takes lower values ( ∼ K) comparing withthe ANI-DV case ( ∼ K), since the model has to fit the dataemploying a polar geometry (N) instead of an equatorial ejecta (D).The viewing angle is biased toward larger values, roughly ∼ deg,inconsistent with GRB expectations (Abbott et al. 2017c; Savchenkoet al. 2017). This is the model that gives the largest evidence, within the providedprior bounds. Regarding the dynamical and viscous ejecta compo-nents, the general features are similar to the one of the ANI-DVcase. The dynamical ejected mass is slightly overestimated compar-ing with NR simulations (Perego et al. 2019; Nedora et al. 2019;Endrizzi et al. 2020; Nedora et al. 2020a; Bernuzzi et al. 2020) ofa factor ∼ . The dynamical component is described by a low opac-ity environment for high-latitudes ( κ high ∼ . g − ) and highopacity for low-latitudes ( κ low ∼
11 cm g − ), in agreement withNR simulations (Perego et al. 2019; Nedora et al. 2019; Endrizziet al. 2020; Nedora et al. 2020a; Bernuzzi et al. 2020). These resultsapproximately agree also with other observational estimations (e.g.,Villar et al. 2017b; Cowperthwaite et al. 2017; Abbott et al. 2017d;Coughlin et al. 2018) Furthermore, the ‘D’ component results thefasted ejected shell, validating the interpretation that this contribu-tion is generated at dynamic time-scales. On the other hand, the vis-cous ejecta is characterized by an average opacity ∼ g − andby low velocity ∼ × − c , an order of magnitude smaller then theone of the dynamical ejecta. These results agrees with the studiespresented in Ref. (Radice et al. 2018c) and they contribute to theLCs in the optical band. Regarding the neutrino-driven wind, the posterior distribution forits ejected mass m (N) shows a bimodality and this degeneracy cor-relates with the heating rate parameter (cid:15) . This behavior can be seenin Fig. (5), that shows the marginalized posterior distribution for (cid:15) and for the total ejected mass M ej , defined as M ej = (cid:88) k =D , N , V m ( k ) , (14)where the index k runs over all the involved components. Themarginalized posterior distribution for m (N) has its dominat peak inproximity of . × − M (cid:12) , while the secondary mode is locatedslightly below × − M (cid:12) . Despite the bimodality, the recoveredvalues of m (N) are smaller compared with the same parameter ex-tracted from the ANI-VN analysis. These results are largely consis-tent with aftermath computations (Perego et al. 2014) and with the-oretical expectations (Perego et al. 2017a), as it is for the recoveredvelocity and opacity parameters.Furthermore, also for the ANI-DVN case, the viewing angle isbiased toward larger values, roughly ∼ deg. The same trendis shown by the anisotropic three-component model employed inRef. (Villar et al. 2017b). The posterior distribution for the T Nifloor parameter peaks around ∼ , while, the temperature T LAfloor isconstrained around the lower bound,
500 K . In this section, we apply the information coming from NR fittingformulae (Nedora et al. 2020b,c) to the posterior distribution of thepreferred kN model (ANI-DVN), in order to infer the mass ratio andthe reduced tidal parameter of the BNS source. Subsequently, wecombine the kN and GW results to derive constraints on the radius R . of an irrotational NS of 1.4 M (cid:12) . MNRAS000
500 K . In this section, we apply the information coming from NR fittingformulae (Nedora et al. 2020b,c) to the posterior distribution of thepreferred kN model (ANI-DVN), in order to infer the mass ratio andthe reduced tidal parameter of the BNS source. Subsequently, wecombine the kN and GW results to derive constraints on the radius R . of an irrotational NS of 1.4 M (cid:12) . MNRAS000 , 000–000 (0000)
T2017gfo: model selection of multi-component kNe ISO − DVANI − DVANI − DVN . . . . . log (cid:0) (cid:15) / erg s − g − (cid:1) . . . . . . . . . M e j (cid:2) − M (cid:12) (cid:3) Figure 5.
Marginalized posterior distribution of heating parameter (cid:15) andtotal ejected mass M ej for three selected cases: ISO-DV (blue), ANI-DV(yellow) and ANI-DVN (green). The heating parameter (cid:15) is plotted usingthe logarithm to base 10 in order to evince the recovered orders of magni-tude. The total mass M ej is computed extending the sum to all the involvedcomponents. A BNS is characterized by the masses of the two objects, m and m , and by the tidal quadrupolar polarizability coefficients, Λ i = 23 k ,i C − i , (15)where k ,i is the quadrupolar Love number, C i = Gm i / ( R i c ) the compactness of star, G the gravitational constant, R i the radiusof the star and i = 1 , . Furthermore, we introduce the mass ratio q = m /m ≥ and the reduced tidal parameter ˜Λ as: ˜Λ = 1613 ( q + 12) q Λ + (1 + 12 q )Λ (1 + q ) . (16)The NR fits presented in Ref. (Nedora et al. 2020c) use simulationstargeted to GW170817 (Perego et al. 2019; Endrizzi et al. 2020; Ne-dora et al. 2019; Bernuzzi et al. 2020; Nedora et al. 2020a) and givethe mass m (D) and velocity v rms(D) of the dynamical ejecta as func-tions of the BNS parameters ( q, ˜Λ) . In order to recover the posteriordistribution of the latter, we adopt a resampling method, similar tothe procedure presented in Ref. (Coughlin et al. 2017, 2018): a sam-ple ( q, ˜Λ) is extracted from the prior distribution , exploiting theranges q ∈ [1 , and ˜Λ ∈ [0 , . Subsequently, the tuple ( q, ˜Λ) is mapped into the dynamical ejecta parameters ( m (D) , v rms(D) ) us-ing the NR formulae presented in Ref. (Nedora et al. 2020c). The The prior distribution is taken uniformly distributed in the tidal parameters ˜Λ ; while, regarding the mass ratio q , we employ a prior distribution uniformin the mass components, that corresponds to a probability density propor-tional to [(1 + q ) /q ] / , analogously to GW analyses (Abbott et al. 2017a;Gamba et al. 2020a). AT2017gfoGW170817Combined ˜Λ . . . . . . . . q ALF2APR4BLhDD2H4LS220SFHoSLy4
Figure 6.
Posterior distribution in the (˜Λ , q ) plane. The blue solid lines referto the resampled values extracted from the kN analysis (ANI-DVN). Theorange solid lines refer to the GW results, where the samples have beenreweighted over a flat prior in ˜Λ . The green solid lines are the combinedinference. The contours represent the 90% credible regions The plot showsthe expectations of some representative EOS. likelihood is estimated in the dynamical ejecta parameter space us-ing a kernel density estimation of the marginalized posterior dis-tribution recovered from the preferred model (ANI-DVN). Further-more, since NR relations have non-negligible uncertainties, we in-troduce calibration parameters α , α , such that log m (D) = (1 + α ) · log m (D)fit ( q, ˜Λ) ,v rms(D) = (1 + α ) · v rms(D)fit ( q, ˜Λ) . (17)The calibrations parameters α , are sampled along the other pa-rameters using a normally distributed prior with vanishing meansand standard deviations prescribed by the relative uncertainties ofNR fits equal to . for both. The resampled posterior distributionis marginalized over the calibration parameters. The BNS parameterspace is explored using a Metropolis-Hasting technique. Note thata correct characterization of the fit uncertainty is crucial, since thiscontribution is the largest source of error in the inference of ( q, ˜Λ) .The posterior distribution in the ( q, ˜Λ) plane as obtained fromthe dynamical ejecta properties fitted to AT2017gfo data is shownin Fig. 6. The measurement of the tidal parameter leads to ˜Λ =900 +310 − , with a bimodality in the marginalized posterior distribu-tion, due to the quadratic nature of the employed NR formulae, withmodes ˜Λ ∼ and ˜Λ ∼ . The mass ratio is constrained tobe lower then . at the 90% confidence level. The uncertainties ofthese estimations are larger than those of the GW analyses (Abbottet al. 2017a, 2019a; Gamba et al. 2020a) and the principal source oferror is the uncertainty of the NR fit formulae.Fig. 6 shows also the results coming from the GW170817 anal-ysis extracted from Ref. (Gamba et al. 2020a). The GW parame-ter estimation has been performed on LIGO-Virgo data strain (Ab-bott et al. 2017a, 2019a,b) with the nested sampling provided by MNRAS , 000–000 (0000) M. Breschi et al.
Table 5.
Estimated values of mass ratio q , reduced tidal parameter ˜Λ and NSradius R . measured from the analyses of AT2017gfo and GW170817. The R . are estimated using the relation proposed in Ref. (De et al. 2018; Radice& Dai 2019) and employing the chirp mass posterior distribution comingfrom the GW analysis (Gamba et al. 2020a).Data q ˜Λ R . [km]AT2017gfo ≤ .
54 900 +310 − . +0 . − . GW170817 ≤ .
33 510 +350 − . +1 . − . Combined ≤ .
27 460 +210 − . +0 . − . the bilby pipeline (Ashton et al. 2019) employing the effective-one-body waveform approximant TEOBResumSPA (Nagar et al.2018; Gamba et al. 2020a) and setting the maximum-frequency cut-off to 1 kHz . Furthermore, the GW posterior samples have beenreweighted with a rejection sampling to the prior distributions em-ployed in the kN study, in order to assume the same prior informa-tion for both analyzes.Under the assumption that GW170817 and AT2017gfo are gen-erated by the same physical event, the ( q, ˜Λ) posterior distributionscoming from the two independent analyses can be combined, in or-der to constrain the estimation of the inferred quantities. The jointprobability distribution is computed as the product of the singleterms, p (cid:0) q, ˜Λ (cid:12)(cid:12) d kn , d gw (cid:1) = p (cid:0) q, ˜Λ (cid:12)(cid:12) d kn (cid:1) · p (cid:0) q, ˜Λ (cid:12)(cid:12) d gw (cid:1) , (18)and the samples are extracted from the with a rejection sampling.The combined inference, shown in Fig. 6, leads to a constraint onthe mass ratio of (cid:46) . and on the tidal parameter ˜Λ = 460 +210 − ,at the 90% confidence levels. Imposing these bounds, stiff nuclearEOS, such as DD2, are disfavored. Using the universal relation presented in Ref. (De et al. 2018; Radice& Dai 2019), it is possible to impose a constraint on the radius R . of a NS of . (cid:12) . We employ the marginalized posterior distribu-tion for the (source-frame) chirp mass M = ( m m ) / / ( m + m ) / coming from the GW170817 measurement (Gamba et al.2020a) and the posterior on the tidal parameter ˜Λ obtained withthe joint analyses AT2017gfo+GW170817. We adopt a resamplingtechnique to account for the uncertainties in the universal relation,introducing a Gaussian calibration coefficient with variance pre-scribed by Ref. (De et al. 2018; Radice & Dai 2019). We estimate R . = 12 . +0 . − . km . The presented measurement agrees with theresults coming from literature (Annala et al. 2018; De et al. 2018;Radice & Dai 2019; Coughlin et al. 2019; Abbott et al. 2018b; Raai-jmakers et al. 2020; Capano et al. 2020; Essick et al. 2020; Dietrichet al. 2020) and its overall error at σ level corresponds roughly to
500 m .In Fig. 7, the R . estimation is compared with the mass-radiuscurves from a sample of nuclear EOS. Our bounds impose observa-tional constraints on the nuclear EOS, excluding both very stiff EOS, On the one hand, this choice minimizes waveform systematics (Gambaet al. 2020b). On the other hand, it implies slightly larger statistical uncertain-ties on the reduced tidal parameters. Hence, our results are more conservativethan previous multimessenger analyses in the treatment of uncertainties ofGW data. R [km] . . . . . . . . . M [ M (cid:12) ] ALF2APR4 BLhDD2 H4LS220 SFHoSLy4
AT2017gfo + GW170817+ Y e + Y e + M disk Figure 7.
Posterior distribution of the radius R . estimated with the joinedinference of AT2017gfo and GW170817 plotted on top of the mass-radiusrelations coming from a sample of nuclear EOS (dashed lines). The blue solidline is computed using the mass and velocity information of the dynamicalcomponent, the orange solid curve takes into account also the contribution ofthe electron fraction and the green solid line is the result with the additionalinclusion of the disk mass information. such as DD2, BHB Λ φ and MS1b, and very soft equations, such as2B. We conduct two further analyses, in order to show that the contri-bution of additional NR information can improve the previous esti-mation. In the first case, we take into account the contribution of theelectron fraction; while, in the second, we include the informationon the disk mass. These studies are discussed in the following para-graphs and they are intended to represent proofs-of-principle anal-yses, since they involve extra assumptions on the ejecta parametersand their relation with the EOS properties. A more accurate mappingbetween these quantities will be discussed in a further study.
From NR simulations, it is possible to estimate the average elec-tron fraction, Y e , of the dynamical ejecta (Nedora et al. 2020b,c).This quantity is the ratio of the net number of electrons to the nu-mer of baryons and it is strictly related with the opacity of theshell (Lippuner & Roberts 2015; Miller et al. 2019; Perego et al.2019), since it mostly determines the nucleosynthesis yields in lowentropy, neutron-rich matter. We compute the average opacity ¯ κ of ashell as the integral of the opacity over the polar angle weighted onthe mass distribution, ¯ κ = 1 m (cid:90) π (cid:37) ( θ ) κ ( θ ) sin θ d θ . (19) MNRAS000
From NR simulations, it is possible to estimate the average elec-tron fraction, Y e , of the dynamical ejecta (Nedora et al. 2020b,c).This quantity is the ratio of the net number of electrons to the nu-mer of baryons and it is strictly related with the opacity of theshell (Lippuner & Roberts 2015; Miller et al. 2019; Perego et al.2019), since it mostly determines the nucleosynthesis yields in lowentropy, neutron-rich matter. We compute the average opacity ¯ κ of ashell as the integral of the opacity over the polar angle weighted onthe mass distribution, ¯ κ = 1 m (cid:90) π (cid:37) ( θ ) κ ( θ ) sin θ d θ . (19) MNRAS000 , 000–000 (0000)
T2017gfo: model selection of multi-component kNe AT2017gfoGW170817Combined ˜Λ . . . . . . . . q ALF2APR4BLhDD2H4LS220SFHoSLy4
Figure 8.
Posterior distribution in the (˜Λ , q ) plane, analogously to Fig. 6,including the contribution of the electron fraction Y e . Imposing the assumptions on the profiles of the dynamical ejecta,we get ¯ κ (D) = (cid:18)
12 + 1 π (cid:19) κ (D)low + (cid:18) − π (cid:19) κ (D)high . (20)Thanks to this definition, it is possible to map the opacity ¯ κ intothe electron fraction Y e , using the relation presented in Ref. (Tanakaet al. 2019). Subsequently, the Y e can be related with the BNS pa-rameters ( q, ˜Λ) , using NR fit formulae (Nedora et al. 2020c). Weintroduce an additional calibration parameter α , such that Y e = (1 + α ) · Y e fit ( q, ˜Λ) , (21)with a Gaussian prior with mean zero and standard deviation of . .In this way it is possible to take into account also the contribution ofthe opacity posterior distribution, introducing additional constraintson the inference of the NS matter.The results are shown in Fig. 8. This further contribution has astrong effect on the mass ratio, constraining it to be (cid:46) . . Thiseffect is motivated by the fact that high-mass-ratio BNS merg-ers are expected to have Y e (cid:46) . (Bernuzzi et al. 2020; Ne-dora et al. 2020b). The recovered electron fraction correspond to Y e = 0 . +0 . − . . Regarding the tidal parameter, the Y e informa-tion affects the importance of the modes, improving the agreementwith GW estimations (Abbott et al. 2017a, 2019a; Gamba et al.2020a), and it reduces the support of the posterior distribution, lead-ing to an estimation of ˜Λ = 480 +550 − . Combining kN and GW pos-terior distribution, we estimate an upper bound on the mass ratioof . and a tidal parameter ˜Λ = 465 +175 − , that corresponds to R . = 12 . +0 . − . km , at the 90% confidence level. The employed kN model contains information also on the bary-onic wind ejecta. These components are expected to be generated
AT2017gfoGW170817Combined ˜Λ . . . . . . . . q ALF2APR4BLhDD2H4LS220SFHoSLy4
Figure 9.
Posterior distribution in the (˜Λ , q ) plane, analogously to Fig. 6,including the contributions of electron fraction Y e and disk mass M disk . by the disk that surrounds the remnant (Kasen et al. 2015; Metzger& Fern´andez 2014; Just et al. 2015), if present. The disk mass canbe estimated from NR simulations as function of the BNS param-eters ( q, ˜Λ) , albeit with large uncertainties (Radice et al. 2018b,d;Nedora et al. 2020c). We map a fraction ξ of the disk mass M disk into the mass of the baryonic wind components, m (V) + m (N) = ξ · M disk . (22)The mass fraction ξ is sampled along the other parameters with auniform prior in the range [0 . , . . We include the disk mass infor-mation together with the electron fraction contribution, previouslydiscussed.The results are shown in Fig. 9. The disk mass contributionslightly reinforces the constraint on the mass ratio posterior, giv-ing the 90% confidence level for q = 1 . . The distribution of thetidal parameter ˜Λ is sparser with respect to the case discussed inSec. 5.3.1, due to the correlations induced by the M disk formula.The electron fraction results Y e = 0 . +0 . − . ; while, the mass frac-tion corresponds to ξ = 0 . +0 . − . . The joined inference with theGW posterior leads to a mass ratio (cid:46) . and a tidal parameter of ˜Λ = 430 +180 − , at the 90% confidence. This result can be translatedin a radius of R . = 11 . +0 . − . km . In this paper, we have performed informative model selection onkN observations within a Bayesian framework applied to the caseof AT2017gfo, the kN associated with the BNS merger GW170817.We have then combined the posteriors obtained from the kN obser-vation with the ones extracted from the GW signal and with NR-based fitting formulae on the ejecta and remnant properties to settight constraints on the NS radius and EOS.From the analysis of AT2017gfo, the anisotropic description of
MNRAS , 000–000 (0000) M. Breschi et al. the ejecta components is strongly preferred with respect to isotropicprofiles, with a logarithmic Bayes’ factor of the order of ∼ .Moreover, the favored model is the three-component kN constitutedby a fast dynamical ejecta (comprising both a red-equatorial and ablue-polar portion), a slow isotropic shell and a polar wind. For thebest model, the dynamical ejected mass overestimates of a factortwo the theoretical expectation coming from NR simulations (Peregoet al. 2019; Nedora et al. 2019; Endrizzi et al. 2020; Nedora et al.2020a; Bernuzzi et al. 2020). This biases can be explained by con-sidering the effect of the spiral-wave wind (Nedora et al. 2019) andtaking into account the correlations between the extrinsic parame-ters. The recovered velocity of the dynamical component agrees withNR simulations (Perego et al. 2019; Nedora et al. 2019; Endrizziet al. 2020; Nedora et al. 2020a; Bernuzzi et al. 2020), reinforcingthe interpretation of this ejecta component. The intrinsic propertiesof the dynamical ejecta component are in agreement with previousresults (Villar et al. 2017b; Coughlin et al. 2019). Regarding the sec-ular winds, the neutrino-driven mass and velocity are compatiblewith the calculations of Ref. (Perego et al. 2014, 2017a). The vis-cous component is the slowest contribution and is broadly compati-ble with the estimates of Ref. (Radice et al. 2018c), that are inferredfrom NR and other disc simulations. The viewing angle resultingfrom the preferred kN model is larger than the one deduced fromindependent analysis (Abbott et al. 2017c; Savchenko et al. 2017;Ghirlanda et al. 2019), and also different from the one obtained byprevious application of the same kN model (Perego et al. 2017a). Inthe latter case, and differently from the present analysis, the profileof the viscous ejecta was assumed to be mostly distributed across theequatorial angle. This discrepancy confirms the non-trivial depen-dence of the light curves from the ejecta geometry and distributions.Under a modeling perspective, current kN description containslarge theoretical uncertainties, such as thermalization effects, heat-ing rates and energy-dependent photon opacities, e.g. (Zhu et al.2020). These effects propagate into systematic biases in the globalparameters of the model, as shown in the posterior distributions forluminosity distance D L and heating rate parameter (cid:15) . Hence, thedevelopment and the improvements of kN templates is an urgenttask in order to conduct reliable and robust analyses in the future.We use of the preferred kN model to constrain the properties ofthe progenitor BNS and the EOS of dense, cold matter. Combiningthe kN measurement with the information coming from NR simu-lations, the ejecta properties are mapped in terms of mass ratio andreduced tidal deformability of the binary progenitor. Subsequently,this information is combined with the measurements of the GW data.The joint kN+GW analysis constrains the reduced tidal parameter to ˜Λ = 460 +210 − and the mass ratio of the BNS system to be lower than . , at the 90% credible level. Furthermore, the joint analysis pre-dicts a radius for a NS of . (cid:12) approximately of R . ≈ . with an uncertainty of ∼
500 m at one- σ level. The R . estimationcan be further improved including additional physical informationextracted from the kN model in the inferred model, such as the elec-tron fraction of the dynamical ejecta and the mass of the disk aroundthe merger remnant. Figure 10 summarizes ours and the current esti-mations of R . extracted from literature (Annala et al. 2018; Radiceet al. 2018b; De et al. 2018; Radice & Dai 2019; Coughlin et al.2019; Abbott et al. 2018b; Raaijmakers et al. 2020; Capano et al.2020; Essick et al. 2020; Dietrich et al. 2020).In addition to the kN modeling uncertainties discussed above, an-other source of error of our estimates is the accuracy of the NR for-mulae. The relations employed here used exclusively targeted dataand simulations with state-of-art treatment of microphysical EOSand neutrino treatment (Perego et al. 2019; Nedora et al. 2019; En- R a d i cee t a l . ( ) A nn a l a e t a l . ( ) D ee t a l . ( ) A bb o tt e t a l . ( ) R a d i ce & D a i ( ) C o u g h li n e t a l . ( ) C a p a n o e t a l . ( ) R aa i j m a k e r s e t a l . ( ) E ss i c k e t a l . ( ) D i e t r i c h e t a l . ( ) T h i s w o r k T h i s w o r k T h i s w o r k R . M (cid:12) [ k m ] AT2017gfo + GW170817+ Y e + Y e + M disk Figure 10.
Summary plot of the current estimations of R . . The reportedvalues are the means and the 90% credible regions extracted from Refs. (An-nala et al. 2018; Radice et al. 2018b; De et al. 2018; Radice & Dai 2019;Coughlin et al. 2019; Abbott et al. 2018b; Raaijmakers et al. 2020; Ca-pano et al. 2020; Essick et al. 2020; Dietrich et al. 2020). The dashedline and the shadowed area are respectively the average over all the cur-rent estimations and the respective 90% credible region, corresponding to R . = 12 . +1 . − . km. drizzi et al. 2020; Nedora et al. 2020a; Bernuzzi et al. 2020). How-ever, the simulation sample is limited to about hundrends of simu-lations, with fitting errors that could be reduced by considering dataat even higher grid resolutions (Nedora et al. 2020c). For example,assuming all the fit formulae to be exact (i.e. removing all calibra-tion terms), it will be possible to infer the ˜Λ parameter from a kNobservation with an accuracy of the order of 10, that corresponds toa constraint on the radius R . of roughly
100 m . ACKNOWLEDGEMENTS
M.B. and S.B. acknowledges support by the European Union’sH2020 under ERC Starting Grant, grant agreement no. BinGraSp-714626. D.R. acknowledges support from the U.S. Department ofEnergy, Office of Science, Division of Nuclear Physics under AwardNumber(s) DE-SC0021177 and from the National Science Foun-dation under Grant No. PHY-2011725. The computational exper-iments were performed on the ARA cluster at Friedrich SchillerUniversity Jena supported in part by DFG grants INST 275/334-1 FUGG and INST 275/363-1 FUGG, and ERC Starting Grant,grant agreement no. BinGraSp-714626. Data postprocessing wasperformed on the Virgo “Tullio” server at Torino supported byINFN. This research has made use of data, software and/or webtools obtained from the Gravitational Wave Open Science Center( ), a service of LIGOLaboratory, the LIGO Scientific Collaboration and the Virgo Col-laboration. LIGO Laboratory and Advanced LIGO are funded by theUnited States National Science Foundation (NSF) as well as the Sci-
MNRAS000
MNRAS000 , 000–000 (0000)
T2017gfo: model selection of multi-component kNe ence and Technology Facilities Council (STFC) of the United King-dom, the Max-Planck-Society (MPS), and the State of Niedersach-sen/Germany for support of the construction of Advanced LIGO andconstruction and operation of the GEO600 detector. Additional sup-port for Advanced LIGO was provided by the Australian ResearchCouncil. Virgo is funded, through the European Gravitational Ob-servatory (EGO), by the French Centre National de Recherche Sci-entifique (CNRS), the Italian Istituto Nazionale della Fisica Nucle-are (INFN) and the Dutch Nikhef, with contributions by institutionsfrom Belgium, Germany, Greece, Hungary, Ireland, Japan, Monaco,Poland, Portugal, Spain. DATA AVAILABILITY
The observational data underlying this article were provided by (Vil-lar et al. 2017b) under license. The posterior samples presented inthis work will be shared on request to the corresponding author.
REFERENCES
Aasi J., et al., 2015, Class. Quant. Grav., 32, 074001Abbott B. P., et al., 2017a, Phys. Rev. Lett., 119, 161101Abbott B. P., et al., 2017b, Astrophys. J., 848, L12Abbott B. P., et al., 2017c, Astrophys. J., 848, L13Abbott B. P., et al., 2017d, Astrophys. J., 850, L39Abbott B. P., et al., 2018a, Living Rev. Rel., 21, 3Abbott B. P., et al., 2018b, Phys. Rev. Lett., 121, 161101Abbott B. P., et al., 2019a, Phys. Rev., X9, 011001Abbott B. P., et al., 2019b, Phys. Rev., X9, 031040Acernese F., et al., 2015, Class. Quant. Grav., 32, 024001Agathos M., Zappa F., Bernuzzi S., Perego A., Breschi M., Radice D., 2020,Phys. Rev., D101, 044006Ajello M., et al., 2016, Astrophys. J., 819, 44Annala E., Gorda T., Kurkela A., Vuorinen A., 2018, Phys. Rev. Lett., 120,172703Arnett W. D., 1982, ”Astrophys. J.”, 253, 785Ashton G., et al., 2019, Astrophys. J. Suppl., 241, 27Barbieri C., Salafia O. S., Perego A., Colpi M., Ghirlanda G., 2020, Eur.Phys. J., A56, 8Barnes J., Kasen D., Wu M.-R., Martinez-Pinedo G., 2016, Astrophys. J.,829, 110Barnes J., Zhu Y., Lund K., Sprouse T., Vassh N., McLaughlin G.,Mumpower M., Surman R., 2020Bauswein A., Goriely S., Janka H.-T., 2013, Astrophys.J., 773, 78Bauswein A., Stergioulas N., Janka H.-T., 2014, Phys.Rev., D90, 023002Bauswein A., Just O., Janka H.-T., Stergioulas N., 2017, Astrophys. J., 850,L34Bernuzzi S., Nagar A., Balmelli S., Dietrich T., Ujevic M., 2014,Phys.Rev.Lett., 112, 201101Bernuzzi S., Nagar A., Dietrich T., Damour T., 2015a, Phys.Rev.Lett., 114,161103Bernuzzi S., Dietrich T., Nagar A., 2015b, Phys. Rev. Lett., 115, 091101Bernuzzi S., et al., 2020, Mon. Not. Roy. Astron. Soc.Bovard L., Martin D., Guercilena F., Arcones A., Rezzolla L., Korobkin O.,2017, Phys. Rev., D96, 124005Breschi M., Bernuzzi S., Zappa F., Agathos M., Perego A., Radice D., NagarA., 2019, Phys. Rev., D100, 104029Capano C. D., et al., 2020, Nature Astron., 4, 625Cardelli J. A., Clayton G. C., Mathis J. S., 1989, ApJ, 345, 245Chornock R., et al., 2017, Astrophys. J., 848, L19Coughlin M., Dietrich T., Kawaguchi K., Smartt S., Stubbs C., Ujevic M.,2017Coughlin M. W., et al., 2018, Mon. Not. Roy. Astron. Soc., 480, 3871 Coughlin M. W., Dietrich T., Margalit B., Metzger B. D., 2019, Mon. Not.Roy. Astron. Soc., 489, L91Coulter D. A., et al., 2017, ScienceCowperthwaite P. S., et al., 2017, Astrophys. J., 848, L17Damour T., Nagar A., Villain L., 2012, Phys.Rev., D85, 123007De S., Finstad D., Lattimer J. M., Brown D. A., Berger E., Biwer C. M.,2018, Phys. Rev. Lett., 121, 091102Decoene V., Gu´epin C., Fang K., Kotera K., Metzger B., 2020, JCAP, 04,045Dietrich T., Bernuzzi S., Br¨ugmann B., Tichy W., 2018, in 2018 26th Eu-romicro International Conference on Parallel, Distributed and Network-based Processing (PDP). pp 682–689 ( arXiv:1803.07965 ),doi:10.1109/PDP2018.2018.00113, https://inspirehep.net/record/1663472/files/1803.07965.pdf
Dietrich T., Coughlin M. W., Pang P. T., Bulla M., Heinzel J., Issa L., TewsI., Antier S., 2020Endrizzi A., et al., 2020, Eur. Phys. J. A, 56, 15Essick R., Tews I., Landry P., Reddy S., Holz D. E., 2020, Phys. Rev. C, 102,055803Even W., et al., 2020, Astrophys. J., 899, 24Fern´andez R., Quataert E., Schwab J., Kasen D., Rosswog S., 2015, Mon.Not. Roy. Astron. Soc., 449, 390Fern´andez R., Metzger B. D., 2013, Mon. Not. Roy. Astron. Soc., 435, 502Fujibayashi S., Kiuchi K., Nishimura N., Sekiguchi Y., Shibata M., 2018,Astrophys. J., 860, 64Fujibayashi S., Wanajo S., Kiuchi K., Kyutoku K., Sekiguchi Y., Shibata M.,2020Gamba R., Bernuzzi S., Nagar A., 2020aGamba R., Breschi M., Bernuzzi S., Agathos M., Nagar A., 2020bGhirlanda G., et al., 2019, Science, 363, 968Goodman J., Weare J., 2010, Communication in Applied Mathematics andComputational Science, 5, 65–80Grossman D., Korobkin O., Rosswog S., Piran T., 2014, Mon. Not. Roy.Astron. Soc., 439, 757Hajela A., et al., 2019, Astrophys. J. Lett., 886, L17Hotokezaka K., Kiuchi K., Kyutoku K., Okawa H., Sekiguchi Y.-i., et al.,2013, Phys.Rev., D87, 024001Hotokezaka K., Beniamini P., Piran T., 2018, Int. J. Mod. Phys., D27,1842005Just O., Obergaulinger M., Janka H. T., 2015, Mon. Not. Roy. Astron. Soc.,453, 3386Kasen D., Barnes J., 2018Kasen D., Badnell N. R., Barnes J., 2013, Astrophys. J., 774, 25Kasen D., Fern´andez R., Metzger B., 2015, Mon. Not. Roy. Astron. Soc.,450, 1777Kasen D., Metzger B., Barnes J., Quataert E., Ramirez-Ruiz E., 2017, NatureKawaguchi K., Shibata M., Tanaka M., 2019, ] 10.3847/1538-4357/ab61f6Korobkin O., Rosswog S., Arcones A., Winteler C., 2012, Mon. Not. Roy.Astron. Soc., 426, 1940Lippuner J., Roberts L. F., 2015, Astrophys. J., 815, 82Margalit B., Metzger B. D., 2017, Astrophys. J., 850, L19Martin D., Perego A., Arcones A., Thielemann F.-K., Korobkin O., RosswogS., 2015, Astrophys. J., 813, 2Metzger B. D., 2020, Living Rev. Rel., 23, 1Metzger B., Berger E., 2012, Astrophys.J., 746, 48Metzger B. D., Fern´andez R., 2014, Mon.Not.Roy.Astron.Soc., 441, 3444Metzger B., Piro A., Quataert E., 2008, Mon.Not.Roy.Astron.Soc., 390, 781Metzger B. D., Arcones A., Quataert E., Martinez-Pinedo G., 2010, Mon.Not. Roy. Astron. Soc., 402, 2771Miller J. M., et al., 2019, Phys. Rev., D100, 023008Nagar A., et al., 2018, Phys. Rev., D98, 104052Nakar E., Gottlieb O., Piran T., Kasliwal M. M., Hallinan G., 2018, Astro-phys. J., 867, 18Nedora V., Bernuzzi S., Radice D., Perego A., Endrizzi A., Ortiz N., 2019,Astrophys. J., 886, L30Nedora V., et al., 2020cNedora V., et al., 2020aNedora V., et al., 2020b MNRAS , 000–000 (0000) M. Breschi et al.
Nelson B., Ford E. B., Payne M. J., 2013, The Astrophysical Journal, Sup-plement Series, 210, 11Nicholl M., et al., 2017, Astrophys. J., 848, L18Nynka M., Ruan J. J., Haggard D., Evans P. A., 2018, Astrophys. J. Lett.,862, L19Oechslin R., Janka H.-T., Marek A., 2006, Astron.Astrophys.Perego A., Rosswog S., Cabezon R., Korobkin O., Kaeppeli R., et al., 2014,Mon.Not.Roy.Astron.Soc., 443, 3134Perego A., Radice D., Bernuzzi S., 2017a, Astrophys. J., 850, L37Perego A., Yasin H., Arcones A., 2017b, J. Phys., G44, 084007Perego A., Bernuzzi S., Radice D., 2019, Eur. Phys. J., A55, 124Pian E., et al., 2017, NaturePiran T., Nakar E., Rosswog S., 2013, Mon. Not. Roy. Astron. Soc., 430,2121Pozzo W. D., Veitch J., , https://github.com/johnveitch/cpnest , doi:10.5281/zenodo.835874Raaijmakers G., et al., 2020, Astrophys. J. Lett., 893, L21Radice D., Dai L., 2019, Eur. Phys. J., A55, 50Radice D., Galeazzi F., Lippuner J., Roberts L. F., Ott C. D., Rezzolla L.,2016, Mon. Not. Roy. Astron. Soc., 460, 3255Radice D., Bernuzzi S., Del Pozzo W., Roberts L. F., Ott C. D., 2017, Astro-phys. J., 842, L10Radice D., Perego A., Bernuzzi S., Zhang B., 2018a, Mon. Not. Roy. Astron.Soc., 481, 3670Radice D., Perego A., Zappa F., Bernuzzi S., 2018b, Astrophys. J., 852, L29Radice D., Perego A., Hotokezaka K., Bernuzzi S., Fromm S. A., RobertsL. F., 2018c, Astrophys. J. Lett., 869, L35Radice D., Perego A., Hotokezaka K., Fromm S. A., Bernuzzi S., RobertsL. F., 2018d, Astrophys. J., 869, 130Roberts L. F., Kasen D., Lee W. H., Ramirez-Ruiz E., 2011, Astrophys.J.,736, L21Rosswog S., Piran T., Nakar E., 2013, Mon. Not. Roy. Astron. Soc., 430,2585Rosswog S., Korobkin O., Arcones A., Thielemann F. K., Piran T., 2014,Mon. Not. Roy. Astron. Soc., 439, 744Rosswog S., Sollerman J., Feindt U., Goobar A., Korobkin O., Wollaeger R.,Fremling C., Kasliwal M. M., 2018, Astron. Astrophys., 615, A132Savchenko V., et al., 2017, Astrophys. J., 848, L15Shibata M., Taniguchi K., Uryu K., 2005, Phys. Rev., D71, 084021Siegel D. M., Metzger B. D., 2018, Astrophys. J., 858, 52Sivia D. S., Skilling J., 2006, Data Analysis - A Bayesian Tutorial, 2nd edn.Oxford Science Publications, Oxford University PressSkilling J., 2006, Bayesian Anal., 1, 833Smartt S. J., et al., 2017, NatureTanaka M., et al., 2017, Publ. Astron. Soc. Jap.Tanaka M., Kato D., Gaigalas G., Kawaguchi K., 2019Tanvir N. R., et al., 2017, Astrophys. J., 848, L27Valenti S., et al., 2017, Astrophys. J., 848, L24Veitch J., et al., 2015, Phys. Rev., D91, 042003Villar V. A., Berger E., Metzger B. D., Guillochon J., 2017a, Astrophys. J.,849, 70Villar V. A., et al., 2017b, Astrophys. J., 851, L21Winkler C., Diehl R., Ubertini P., Wilms J., 2011, Space Science Reviews,161, 149–177Wollaeger R. T., et al., 2018, Mon. Not. Roy. Astron. Soc., 478, 3298Wu M.-R., Fern´andez R., Mart´inez-Pinedo G., Metzger B. D., 2016, Mon.Not. Roy. Astron. Soc., 463, 2323Zappa F., Bernuzzi S., Pannarale F., Mapelli M., Giacobbo N., 2019, Phys.Rev. Lett., 123, 041102Zhu Y., Lund K., Barnes J., Sprouse T., Vassh N., McLaughlin G.,Mumpower M., Surman R., 2020de Jes´us Mendoza-Temis J., Wu M.-R., Martinez-Pinedo G., Langanke K.,Bauswein A., Janka H.-T., 2015, Phys. Rev., C92, 055805MNRAS000