Tight multi-messenger constraints on the neutron star equation of state from GW170817 and a forward model for kilonova light curve synthesis

We present a rapid analytic framework for predicting kilonova light curves following neutron star (NS) mergers, where the main input parameters are binary-based properties measurable by gravitational wave detectors (chirp mass and mass ratio, orbital inclination) and properties dependent on the nuclear equation of state (tidal deformability, maximum NS mass). This enables synthesis of a kilonova sample for any NS source population, or determination of the observing depth needed to detect a live kilonova given gravitational wave source parameters in low latency. We validate this code, implemented in the public MOSFiT package, by fitting it to GW170817. A Bayes factor analysis overwhelmingly ($B>10^{10}$) favours the inclusion of an additional luminosity source in addition to lanthanide-poor dynamical ejecta during the first day. This is well fit by a shock-heated cocoon model, though differences in the ejecta structure, opacity or nuclear heating rate cannot be ruled out as alternatives. The emission thereafter is dominated by a lanthanide-rich viscous wind. We find the mass ratio of the binary is $q=0.92\pm0.07$ (90% credible interval). We place tight constraints on the maximum stable NS mass, $M_{\rm TOV}=2.17^{+0.08}_{-0.11}$ M$_\odot$. For a uniform prior in tidal deformability, the radius of a 1.4 M$_\odot$ NS is $R_{1.4}\sim 10.7$ km. Re-weighting with a prior based on equations of state that support our credible range in $M_{\rm TOV}$, we derive a final measurement $R_{1.4}=11.06^{+1.01}_{-0.98}$ km. Applying our code to the second gravitationally-detected neutron star merger, GW190425, we estimate that an associated kilonova would have been fainter (by $\sim0.7$ mag at one day post-merger) and declined faster than GW170817, underlining the importance of tuning follow-up strategies individually for each GW-detected NS merger.


INTRODUCTION
The first binary neutron star (NS) merger detected through its gravitational wave emission, GW170817 (Abbott et al. 2017a), was accompanied by both a short gamma-ray burst (GRB; Goldstein et al. 2017;Savchenko et al. 2017) and an optical counterpart (Coulter et al. 2017;Soares-Santos et al. 2017;Valenti et al. 2017;Arcavi et al. 2017;Tanvir et al. 2017;Lipunov et al. 2017). Thermal transients from these mergers had long been predicted in the form of a kilonova (Li & Paczyński 1998;Rosswog et al. 1999;Metzger et al. 2010): a non-relativistic outflow heated by the decays of heavy nuclei formed through rapid neutron captures (the so-called 'r-process'; Lattimer & Schramm 1974;Eichler et al. 1989;Davies et al. 1994;Freiburghaus et al. 1999). Extensive optical and infrared follow-up of GW170817 by many groups showed photometric and spectroscopic properties that were remarkably consistent with the expectations for an r-process kilonova (Andreoni et al. 2017;Arcavi et al. 2017;Chornock et al. 2017;Cowperthwaite et al. 2017;Díaz et al. 2017;Drout et al. 2017;Evans et al. 2017;Hu et al. 2017;Kasliwal et al. 2017;McCully et al. 2017;Nicholl et al. 2017a;Pian et al. 2017;Pozanenko et al. 2018; ★ Contact e-mail: m.nicholl.1.bham.ac.uk † NASA Einstein Fellow Shappee et al. 2017;Smartt et al. 2017;Tanvir et al. 2017;Troja et al. 2017;Utsumi et al. 2017;Valenti et al. 2017). Studying these sources is crucial for understanding cosmic chemical evolution. The total ejected mass powering the GW170817 kilonova was inferred to be ∼ 0.05 M , with at least two components of different compositions and hence opacities (e.g. Villar et al. 2017;Kasen et al. 2017;Margutti & Chornock 2020): one containing light rprocess elements (atomic mass number 140), and one containing the heavier lanthanides, which absorb optical photons in the forest of transitions in their open -shells ; . The large total mass and broad range of mass numbers covering all three r-process abundance peaks (Burbidge et al. 1957) suggested that NS mergers may be the dominant production site of heavy nuclei in the Universe.
Another central goal of multi-messenger astronomy is to determine the NS equation of state (EoS), i.e. the relation between pressure and density. This provides unique tests of nuclear physics beyond the saturation density of nuclear matter (e.g. Lattimer & Prakash 2016). Because the EoS determines the NS radius and maximum stable mass, measuring these quantities from astrophysical observation directly constrains the EoS. For example, a viable EoS must support the existence of the most massive known pulsars, with mass 2 M (Demorest et al. 2010; Antoniadis et al. 2013;Cromartie et al. 2020). NS mergers provide an excellent laboratory for EoS measurements.
The degree of tidal deformation imprinted in the GW signal (Flanagan & Hinderer 2008;Damour et al. 2012) and the strength of shocks at the contact interface (Bauswein et al. 2013b;Hotokezaka et al. 2013) both depend on the radius, typically parameterised as 1.4 , the radius of a 1.4 M NS. Similarly, the lifetime of the massive merger product until collapse to a black hole (BH) depends on the maximum stable NS mass, the Tolman-Oppenheimer-Volkoff mass TOV (Shapiro & Teukolsky 1986;Margalit & Metzger 2017).
So far, most of the EoS constraints for GW170817 have come from the GW data (e.g. Abbott et al. 2018;De et al. 2018;Raithel et al. 2018). However, electromagnetic data have also provided important insight. Nicholl et al. (2017a) pointed out that the inferred combination of large mass and high velocity in the low-lanthanide ejecta was predicted only in the case of relatively small radii (Bauswein et al. 2013b;Shibata & Hotokezaka 2019). More quantitatively, Margalit & Metzger (2017) used the observed electromagnetic signatures of GW170817 to infer that the merged remnant must have collapsed to a BH within a short timescale (∼ 10 ms to ∼ 1 s), pointing to TOV 2.17 M . Bauswein et al. (2017) related this lifetime to the radius, finding 1.4 > 10.68 km.
A major caveat in these analyses is our ability to associate a measured ejecta component with a specific physical origin. NS mergers can eject mass dynamically during the merger, through tidal stripping or shock heating, and after the merger, via winds from an accretion disk or long-lived remnant (see review by Metzger 2019). Each mechanism is expected to produce a different combination of ejecta mass, velocity and composition, and the relative importance of each mechanism depends on the properties of the system before merger (e.g. Bauswein et al. 2013b;Hotokezaka et al. 2013;Sekiguchi et al. 2016;Radice et al. 2018;Shibata & Hotokezaka 2019;Ciolfi & Kalinani 2020).
Most model fits to the electromagnetic data from GW170817 have been specified in terms of the 'post-merger' parameters: i.e. the individual masses, velocities and compositions/opacities of some number of ejecta components. This introduces subjectivity into the origin of each component, and neglects physically-expected correlations between the properties of the different components. Specifying instead a forward model in terms of 'pre-merger' properties -binary masses and mass ratios, and the EoS -would allow one to trace explicitly the ejection mechanism of each component, and exclude models that fit the data well but require unphysical combinations of ejecta parameters. Moreover, analysing the data in terms of pre-merger properties provides a more direct link between the electromagnetic and GW signals.
Progress in combining light curve models with GW inference was first made by Coughlin et al. (2019a), who fit the light curve of the GW170817 kilonova and related the implied ejecta masses to the predictions of numerical relativity simulations for different system mass ratios and combinations of 1.4 and TOV . Including GW information in their priors, they measured 1.4 = 11.3 − 13.5 km. However, their model was unable to match the early blue part of the light curve, possibly because a sufficiently high mass of lowlanthanide ejecta was not predicted by merger simulations at the GW-inferred binary mass of GW170817. This may suggest the need to look at alternative luminosity sources for the early emission, such as cooling of matter shock-heated by the GRB jet (Kasliwal et al. 2017;Piro & Kollmeier 2018;Arcavi 2018). While our manuscript was in preparation, other efforts have built on these techniques. Dietrich et al. (2020) published an updated analysis, using a light curve model that reasonably reproduced the early flux without shock cooling, and found 1.4 = 11.75 ± 0.86 km. An even more recent analysis of GW170817 was published by Breschi et al. (2021), who found 1.4 = 12.2 ± 0.5 km using a three-component kilonova model. A better understanding of NS physics will come with larger samples of kilonovae, but finding these events has proven challenging. Despite subsequent GW detections of binaries hosting possible NSs, such as GW190425 (Abbott et al. 2020c) and GW190814 (Abbott et al. 2020d), and significant effort invested in counterpart searches (Hosseinzadeh et al. 2019;Lundquist et al. 2019;Coughlin et al. 2019b;Gomez et al. 2019;Andreoni et al. 2020b;Vieira et al. 2020;Ackley et al. 2020;Kasliwal et al. 2020;Gompertz et al. 2020;Paterson et al. 2020;Anand et al. 2021), at the end of the third LIGO-Virgo observing run (O3), GW170817 remains the only NS binary detected in both electromagnetic and gravitational waves (but see Graham et al. 2020, for a possible optical counterpart to a BH binary merger). Analysis of short GRB afterglows has yielded several kilonova candidates (Berger et al. 2013;Tanvir et al. 2013;Jin et al. 2016;Troja et al. 2018Troja et al. , 2019Gompertz et al. 2018;Lamb et al. 2019;Fong et al. 2020) and constraining non-detections (Rastinejad et al. 2021;O'Connor et al. 2021), while optical surveys are also looking for kilonovae in their data streams, though without an unambiguous detection to date (McBrien et al. 2020;Andreoni et al. 2020a). However, a step change in our ability to find kilonovae is anticipated, due to planned improvements in the sensitivity of GW detector networks (Abbott et al. 2020b), and the deepest ever time-domain sky survey: the Vera Rubin Observatory Legacy Survey of Space and Time (LSST; LSST Science Collaboration et al. 2009).
In this paper we provide a fast, analytic forward model for kilonova light curves that is completely specified by the binary configuration and EoS-dependent properties. Fitting such a model to observed kilonova data allows one to measure fundamental pre-merger properties even in the absence of a GW signal, or to use GW information when available as priors for multi-messenger inference. Forward modelling the light curves of kilonovae can increase the efficiency of GW-triggered optical searches, by predicting the electromagnetic light curve directly from the GW signal and providing the required observing depth to detect (or rule out) a counterpart at any time since merger (see also Salafia et al. 2017;Barbieri et al. 2020;Stachie et al. 2021). In the context of all-sky surveys like LSST, such a model will allow for kilonova population synthesis (see also Kashyap et al. 2019;Mochkovitch et al. 2021) from a given NS source population, informing strategies for selecting candidate kilonovae from the data stream or allowing one to constrain the NS merger rate even given non-detections.
GW170817 serves as an ideal testing ground for such a model. By fitting the rich observed data set, we will show the power of using GW results as priors, constrain the physical origins of the different emission components, and present new measurements of the mass ratio of the progenitor system, the viewing angle, and the EoS-dependent quantities 1.4 and TOV .
The structure of this paper is as follows. In section 2, we describe the model setup and assumptions. Section 3 shows its application to GW170817 and our main results. We briefly demonstrate the utility of our code for aiding future kilonova searches, using the example of GW190425, in section 4, before concluding in section 5. A set of kilonova model light curves produced using our code is available for download (see Data Availability section).

DESCRIPTION OF MODEL
We have built our model using the Modular Open Source Fitter for Transients ( ; Guillochon et al. 2018  A full description of these parameters is given in section 2. Left: Pre-merger properties of a NS binary. Bold, green font is used to indicate properties that can be measured from GW waveforms. Right: Ejecta properties probed by electromagnetic observations. Left arrows are used to mean 'depends on', in order to highlight the most important parameter sensitivities. The merged remnant can be a long-lived or short-lived NS or a BH, depending on the ratio ( 1 + 2 )/ TOV . The blue ejecta are from shocks at the collisional interface and have an opacity blue ∼ 0.5 cm 2 g −1 , the red ejecta are tidally stripped in the orbital plane and have red ∼ 10 cm 2 g −1 , and the purple ejecta are winds from an accretion disk and have an intermediate opacity that depends on the lifetime of the remnant. Low opacity ejecta may not be visible for large viewing angle . The shock cooling emission arises from a cocoon of width heated at its interface with the GRB jet.
velocities and opacities have already been implemented by Cowperthwaite et al. (2017) and Villar et al. (2017). Therefore our main task here is to define modules that provide these masses, velocities and opacities for a given binary configuration. We also provide new modules to take into account the effects of observer viewing angle and shock cooling. We have made these new modules and the parameter files to implement them publicly available via GitHub 1 , as a new model called ' '. In this section we will describe the main features. Figure 1 provides a visual overview that summarises the geometry and the relations between key parameters. The connection between our modelling framework and the GW observables is discussed in detail in section 2.6.

Dynamical ejecta
During the merger process, mass is ejected dynamically both by tidal forces in the orbital plane Sekiguchi et al. 2015) and by shocks at the contact interface between the two NSs (Bauswein et al. 2013b;Hotokezaka et al. 2013). To estimate the mass of this material, we use a fit by where is the mass of the th NS and is its compactness given a radius . The term 1 ↔ 2 represents exchanging the subscripts, and the best-fit values of the free parameters , , , , can be found in their study. Quantities denoted * are the baryonic masses (differing from the gravitational mass by the binding energy), which we calculate for a given (in solar masses) following Gao et al. (2020): * = + 0.08 2 . The simulations cover the ranges 1 < / < 2 and 0.1 < < 0.23. The typical fractional uncertainty on dyn using this parameterisation was found by  to be ≈ 72%. Similar expressions were provided for the equatorial and polar ejecta velocities, with smaller uncertainties of only 13 − 33%.
Inspection of equation 1 shows that the dynamical ejecta mass depends on terms of the form / . We adopt the convention that 1 is the heavier NS (for consistency with GW analyses, e.g. Abbott et al. 2017aAbbott et al. , 2018, such that the mass ratio of the system is = 2 / 1 < 1. The other important sensitivity in equation 1 is to the compactness , an EoS-dependent quantity. This will be discussed further in section 2.6. The relationship between binary and ejecta masses for NS mergers was recently revisited by Krüger & Foucart (2020) and Nedora et al. (2020), who used a different set of numerical data but found broadly consistent fits. Coughlin et al. (2019a) provided a fitting formula in terms of log( dyn ), which performed similarly well. However, we find this logarithmic form can lead to unrealistically large ejecta masses in some regions of parameter space, which is not optimal for population synthesis.
To better capture the multiple physical processes contributing to the dynamical ejecta, we divide it into two components. The tidal ejecta, concentrated in the equatorial plane, are expected to contain a very low electron fraction (the ratio of protons to total nucleons), < 0.25, allowing the complete r-process synthesis of very heavy nuclei. These ejecta have been termed 'red' due to the high opacity of the resulting lanthanides at optical wavelengths . The polar material has a higher > 0.25, due to the importance of shocks and neutrino heating at higher latitudes (e.g. Bauswein et al. 2013b;Sekiguchi et al. 2016;Metzger 2019), and is correspondingly 'blue' due to an incomplete r-process only up to mass number ∼ 140. Radice et al. (2018) find that the ratio of red to blue ejecta is largely insensitive to the total system mass or the EoS, and depends primarily on the mass ratio of the system.
To estimate the relative mass in the red and blue dynamical ejecta, we use simulations by Sekiguchi et al. (2016). They provide distributions in the dynamical ejecta as a function of for two equations of state, one stiff (large radius) and one soft (compact). We integrate these distributions (in their Figure 5) to determine red ( ) ≡ ( < 0.25)/ (total) and blue ≡ 1 − red . We find that red is insensitive to the EoS but is a steep function of , ranging from 20% for an equal-mass merger ( = 1) to 76% for = 0.86. For 0.8, red ≈ 1, because the less massive NS is tidally disrupted prior to contact, preventing strong shocks (e.g. Bauswein et al. 2013b;Hotokezaka et al. 2013;Lehner et al. 2016;. For > 0.8, we interpolate red between the available simulations using a smooth polynomial fit (shown in appendix Figure A1). The predicted ejecta masses in each component, for a range of NS binaries, are shown in Figure 2. Following Radice et al. (2018), we assume a fixed grey opacity for each component: blue = 0.5 cm 2 g −1 , and red = 10 cm 2 g −1 (Villar et al. 2017;Metzger 2019;Tanaka et al. 2020).

Post-merger ejecta
For most systems, additional mass will be ejected after the merger, driven by neutrino-heated winds from a viscous accretion disk around the remnant or from the surface of a massive remnant NS (Metzger & Fernández 2014). Simulations show that in many cases the mass of post-merger ejecta will exceed that of the dynamical ejecta (Ciolfi et al. 2017;Siegel & Metzger 2017Radice et al. 2018;Fernández et al. 2019).
The properties of the wind ejecta depend on the lifetime of the merger remnant, which may be a stable NS, a supramassive (hypermassive) NS temporarily supported by (differential) rotation, or may promptly collapse to a BH. To first order, longer-lived remnants produce a larger disk mass with a higher (Metzger & Fernández 2014;Lippuner et al. 2017).
We use the analytic expression for the disk mass from Here tot = 1 + 2 is the total binary mass, and thr is the threshold mass for prompt collapse to a BH, parameterised as i.e. proportional to the maximum stable mass of a non-rotating NS and the compactness (Bauswein et al. 2013a). An alternative fit by Radice et al. (2018) was parameterised in terms of the tidal deformability rather than TOV . Following Coughlin et al. (2019a), we prefer to retain the explicit dependence on TOV so that we can try to constrain it observationally. Recent fits by Nedora et al. (2020)  . We note however that Kiuchi et al. (2019) have shown that the mass ratio can play an important role and enhance the disk mass for small . A fraction disk ∼ 0.1 − 0.5 of the disk mass is expected to be ejected by viscously-driven winds, with disk ≈ 0.2 likely typical (Siegel & Metzger 2017;Radice et al. 2018). This fraction can be fixed in our model for population synthesis, or left free to infer its value from fits to observed KNe (although in principle, disk may change as a function of disk mass; De & Siegel 2020). Additional post-merger ejecta may arise through neutrino-or magnetically-driven winds (Radice et al. 2018;Metzger et al. 2018;Ciolfi & Kalinani 2020), which we will discuss in section 2.5.
To determine the velocity and composition of the disk ejecta, we first estimate the remnant properties as a function of tot , TOV and thr using table 3 of Metzger (2019). For tot < TOV , the remnant is indefinitely stable and we assume an ejecta velocity of ≈ 0.1 based on the simulations of Metzger & Fernández (2014). For tot > thr (prompt collapse), the same simulations predict a velocity of ≈ 0.03 . We linearly interpolate between these velocities for intermediate tot . Figure 2 shows the predicted disk ejecta masses and velocities alongside the dynamical ejecta from section 2.1.
To estimate the opacity we use the simulations from Lippuner et al. (2017), who present distributions in their Figure 2 for disks with various remnant lifetimes, ranging from indefinitely stable to prompt collapse. We integrate these distributions to determine the mass-weighted . We then convert this mean into an effective grey opacity, purple , using fits to the − relations from Tanaka et al. (2020) (shown in appendix Figure A2). The inferred range of 0.25 < < 0.38 corresponds to an opacity range of 5.5 purple /(cm 2 g −1 ) 1.5 for the disk ejecta. The 'purple' subscript refers to the fact that this lies intermediate between the red and blue dynamical ejecta.

Luminosity from r-process decay
A semi-analytic model that captures the essential light curve physics of kilonovae was developed by Metzger (2019). This model approximates the heating rate at time , averaged over all r-process radioactive decay chains, as ( ) ∝ th r −1.3 (Korobkin et al. 2012), where r is the r-process mass and th is a time-dependent thermalisation efficiency (Barnes et al. 2016). A module to calculate ( ) in this way was implemented in by Cowperthwaite et al. (2017) and Villar et al. (2017).
To determine the output luminosity, ( ), the ejecta are assumed to have a constant grey opacity and expand homologously with scale velocity ej , allowing for an analytic radiative diffusion solution to the first law of thermodynamics (Arnett 1982). The temperature of the photosphere is initially calculated using the Stefan-Boltzmann law, with radius = ej , until such time as the cooling ejecta start to recombine. At this point, the photosphere recedes into the ejecta while maintaining the recombination temperature, rec . The modules to solve for the radiative diffusion and the photospheric temperature and radius are described in Nicholl et al. (2017b). Cowperthwaite et al. (2017) and Villar et al. (2017) generalised this model to include the sums of multiple ejecta components: the luminosity scale of each component is proportional to its mass, with a diffusion timescale determined also by its velocity and opacity. They allowed rec to vary independently for each of 2-3 ejecta components, finding best fit values in the range 1000 K rec 4000 K. To reduce the number of free parameters (and restrict to binary-based quantities of interest) we assume a fixed rec = 2500 K, based on the theoretically-predicted recombination temperature of ionised r-process ejecta . This falls in the middle of the plausible range found by Villar et al. (2017). We tested the effect of leaving this parameter free, and found a best-fit value rec,free ≈ 2600 ± 200 K, consistent with theory. For completeness, we also tested the effect of using three independent temperatures in our base model; in this case they converged to similar values to Villar et al. (2017), and other parameters remained consistent with the fixed temperature results to better than 1 .

Geometry
Most analytic light curve models assume spherical symmetry. However, in kilonovae the separate ejecta components occupy distinct spatial regions, due to their different ejection mechanisms (as illustrated in Figure 1). The low-medium (i.e. blue and purple) ejecta are expected to escape mostly orthogonal to the orbital plane, where shock heating and weak interactions lower the free neutron flux (e.g. Wanajo et al. 2014;Goriely et al. 2015;Foucart et al. 2016). In contrast, the neutron-rich red ejecta escape parallel to the orbital plane since they form tidally. Simulations show that this geometry can be well approximated as a sphere with conical sections of half-opening angle open ∼ 45 • around the poles. The blue/purple ejecta reside within the conical sections and the red ejecta outside.
This results in a viewing angle dependence in the observed luminosity, colour and timescale of the kilonova Kawaguchi et al. 2018;Wollaeger et al. 2018;Bulla 2019;Korobkin et al. 2020, e.g.), which can be understood intuitively: from viewing angles close to the equator, emission from the low-lanthanide polar ejecta is obscured by the curtain of lanthanide-rich ejecta at lower latitudes, meaning that only a slow-evolving red kilonova is observed. Radiative transfer simulations show that this can reduce the peak luminosity by ≈ 1 − 2 magnitudes (Bulla 2019). Neglecting this effect is problematic both for population synthesis and for fitting observed kilonovae, since an offset in the luminosity introduces a systematic error in measuring the ejected mass. Another subtle viewing angle effect is associated with the polar cavity carved out by a GRB jet (whether failed or successful) that can accompany the merger, as recently studied by Nativi et al. (2021); Klion et al. (2021). This can expose interior regions of the ejecta which are hotter or lower opacity, resulting in brighter or bluer emission. Villar et al. (2017) approximated the viewing angle effect in by multiplying the luminosity of the equatorial ejecta by cos and the polar ejecta by 1 − cos , with the viewing angle a free parameter. Here we implement a new module that treats this more exactly, using the formalism of Darbha & Kasen (2020). This takes a sphere with conical polar caps of half-opening angle open , and scales the luminosity of each component using the projected area of the caps (blue/purple ejecta) or remaining sphere (red ejecta) as a function of and open . We have verified that the change in luminosity with is consistent with the radiative transfer results of Bulla (2019) for the same geometry. We note that this treatment does not account for the jet-interaction viewing angle effect discussed by Klion et al. (2021); Nativi et al. (2021), as this would depend on uncertain properties of the jet. For our fiducial model, we assume open = 45 • and vary the observer viewing angle .

Extensions to the base model
We include in our model two optional final ingredients that can contribute additional luminosity. Although they are tied less explicitly in our formalism to the initial binary parameters, neglecting these   Kasen et al. (2017) for the closest available ejecta masses (assuming lanthanide fractions of 10 −9 , 10 −2 and 10 −1 for the blue, purple and red components, respectively). These models have not been tuned to match our simulations and are intended only to guide the eye. The viewing angle effect is very roughly approximated by adding −0.5 mag to the band and +0.5 mag to the band.
contributions may lead to an underestimate in the model luminosity or to a bias in the inference of other parameters from kilonova data. These components can easily be turned off at run-time, allowing the user to recover the fully predictive baseline model described in the previous sections.
The first is an enhancement in the blue ejecta due to magneticallydriven winds Ciolfi & Kalinani 2020). This mechanism is possible only if the remnant avoids prompt collapse, hence we apply this effect only when the merger mass is below thr (equation 4). We adopt a simple parameterisation, inspired by Coughlin et al. (2019a): blue,tot = blue,dyn / , where 0 < < 1. Therefore fixing = 1 turns off this source of ejecta.
The second effect is shock-heating of the ejecta by a GRB jet. GW170817 was accompanied by a short GRB (as expected in NS mergers; Eichler et al. 1989;Narayan et al. 1992;Berger 2014), and some studies interpreted the early blue luminosity of the kilonova as arising from cooling of this shocked 'cocoon' (Kasliwal et al. 2017;Arcavi 2018). We include this effect using the formulae from Piro & Kollmeier (2018), specifically their equations 38 and 39. The breakout radius is approximated as the time delay between the GW and GRB signals (1.7 s for GW170817; Abbott et al. 2017b) multiplied by the speed of light (Piro & Kollmeier 2018). Following Nakar & Piran (2017), we assume that the shock deposits constant energy per decade of velocity in the ejecta. The isotropic-equivalent luminosity is proportional to the mass of shocked ejecta, which is a fraction 2 c /2 of the total ejecta. The shock can therefore be suppressed by setting the cocoon half-opening angle c equal to zero (or cos c = 1). We assume that the shocked mass is a fraction of the total dynamical ejecta, as it is unclear whether post-merger ejecta will precede the launching of the GRB jet. The output light curves are only weakly sensitive to this assumption, partially due to degeneracy between ej and c . Because the cocoon shock occurs primarily in the polar material, we assume this matter has the same opacity as the blue ejecta component. The luminosity of the shock is inversely proportional to this opacity, so we note that a lanthanide-rich polar component could significantly reduce the importance of the shock emission.
We show examples of light curves calculated for different mass ratios, inclinations and cocoon opening angles in Figure 3. Table 1. Parameters in the model and application to GW170817. Some parameters are relevant only when fitting to data. Priors in square brackets are flat over the specified ranges, otherwise Gaussian. Posterior values for tidal deformability are given in terms ofΛ to facilitate comparison to GW analyses. The 'Surface ejecta' model includes additional blue ejecta from a long-lived remnant; the 'Shocked cocoon' model includes emission from a GRB-heated cocoon; the 'Agnostic' model includes both. The stated best-fit values and uncertainties correspond to the medians and 16th/84th percentiles of the marginalised posteriors. Hydrogen column density in host galaxy (proportional to extinction); White noise in likelihood function; * Additional constraint that derived parameter Λ < 800, see section 2.6

Free parameters and relation to GW observables
After implementing the modules described above, our model depends on the masses of the constituent neutron stars, their compactness, the maximum stable mass TOV , and the viewing angle . The fraction of the disk mass ejected in the purple component can be varied or fixed to a fiducial value. Two further parameters allow one to increase the early blue flux through magnetic winds and/or shock cooling. Finally, we include line-of-sight extinction, and a white noise parameter in our likelihood function when fitting to data (see Nicholl et al. 2017b;Villar et al. 2017;Guillochon et al. 2018).
There is a choice to be made about how best to express some of these parameter dependencies. We opt to use the parameterisation that most closely matches the quantities constrained by GW observations. This allows the user to quickly predict the electromagnetic light curve for a given GW signal, or to use GW information to inform their priors when fitting to an observed multi-messenger source. We therefore parameterise the masses of the system in terms of the 'chirp' mass M = ( 1 2 ) 3/5 ( 1 + 2 ) −1/5 , the quantity most accurately measured from the GW signal, and = 2 / 1 ≤ 1, to which this signal is also somewhat sensitive.
The dependence on the NS compactness is more complicated. For population synthesis, it is simplest to specify a radius , equivalent to implicitly choosing an equation of state (since all NSs of comparable mass should have approximately the same radius). However, when fitting to observed data (i.e. when we want to measure, rather than impose, a NS radius), it is more useful to express in terms of the NS tidal deformability, which measures the responsiveness of the NS to an external tidal field: Λ = (2/3) 2 −5 , where 2 ( ) ≈ 0.05 − 0.15 is the quadrupolar tidal Love number (Hinderer 2008;Postnikov et al. 2010;Damour et al. 2012). The reason is that the observable GW phase is to leading order determined by a specific combination of Λ 1 and Λ 2 (Flanagan & Hinderer 2008;Wade et al. 2014;Favata 2014),

Λ = 16 13
( 1 + 12 2 ) 4 1 Λ 1 + ( 2 + 12 1 ) 4 2 Λ 2 ( 1 + 2 ) 5 , called the binary or effective tidal deformability, which can therefore be constrained with GW data. The relation between Λ and allows us to include GW information in our priors to better measure and hence . A similar approach was adopted by Coughlin et al. (2019a), though they used a different relation proposed by De et al. (2018) to relateΛ to without calculating explicitly. For consistency with GW analyses (Abbott et al. 2018) and to reduce systematic error, we instead sample Λ in terms of the so-called 'symmetric' tidal deformability: Using 'quasi-universal relations' or QURs (largely independent of the EoS, and hence of complicating factors such as 2 ), the antisymmetric tidal deformability Λ a ≡ (Λ 2 − Λ 1 )/2 can be derived from Λ s to ≈ 5% precision (Yagi & Yunes 2016), allowing one to reconstruct Λ 1 , Λ 2 andΛ directly from Λ s . One may then impose the constraint that Λ be consistent with the GW results, using the class in (see Guillochon et al. 2018). The compactness of each NS is then derived from Λ using another QUR: = 0.360 − 0.0355 ln (Λ ) + 0.000705 ln (Λ ) 2 This relation is accurate to ≈ 6.5% across a variety of equations of state (Yagi & Yunes 2017). These are then used in equation 1 to estimate the dynamical ejecta mass. Following a fit to observed data, the NS radii can be derived from the posterior distribution of Λ s following the same procedure and using the constituent NS masses in equation 2.
Our final parameter set therefore consists of 5-11 parameters, which differ slightly between the generative and fitting versions of the model. We list all of these parameters and suggest priors in Table 1. For comparison, the parameterisation of the original kilonova model by Villar et al. (2017) also had 11 free parameters.

APPLICATION TO GW170817
As with any model of NS mergers, GW170817 provides an ideal testing ground for the formalism described above. In this section, we fit the observed optical data to demonstrate (i) how using a binary (rather than ejecta) based parameter set provides additional insight into the physical origin of each luminosity component; (ii) the utility  Best-fit kilonova model for GW170817 without surface ejecta enhancement or GRB shock cooling. Right: Best-fit model with these effects included. Both models provide a good match to the data beyond 1.5 days post-merger, but the fit without shock cooling is too faint by ≈ 0.7 magnitudes (almost a factor 2) during the first day. The model on the right is preferred to the base model with a Bayes factor > 10 19 . Light curve fits for the surface ejecta and shock-only models are shown in the appendix. of this model in probing the NS EoS with kilonova data; and (iii) the extent to which the use (or not) of GW information in the model priors affects our posteriors.
We use data from the Open Kilonova Catalog (Guillochon et al. 2017), compiled by Villar et al. (2017). Given the high level of overlap between the observations from different groups, we fit the following subset, chosen to cover the full range of bands with a comparable (approximately nightly) density of sampling in each: Our priors are given in Table 1 We run our fits on the University of Birmingham cluster, and use (Speagle 2020) to integrate the model evidence and sample the posteriors of the model parameter space. We use the likelihood function where and are the set of observed and model magnitudes, are the errors on the data, and is a white noise free parameter to account for additional uncertainty in the data or model; for a good fit within the observed errors, one therefore finds . Speagle (2020) gives an extensive account of how the evidence is evaluated directly in ; we refer the reader there for details (and a good overview of Bayes' Theorem). In short, random draws from the model priors are used to integrate numerically the hypervolume between shells of constant likelihood. The products of likelihood and volume for these points yield a set of weights that can be summed to determine the total evidence. Points are evolved by replacing those with lowest likelihood; the weights of the final set are proportional to the model posterior (Speagle 2020).

Model comparison and importance of shock cooling
We begin by fitting GW170817 with four variants of the model to determine the relative importance of different effects. These are (i) the baseline model ( = cos c = 1), (ii) a model with blue ejecta enhancement from e.g. magnetic surface winds ( ≤ 1), (iii) a model with shock-cooling of a GRB cocoon (cos c ≤ 1), and (iv) an agnostic model allowing for both effects ( , cos c ≤ 1). The model evidence and marginalised posteriors are given in Table 1. We show resulting light curves compared to the GW170817 data in Figure 4.
We use the Bayes factor, ≡ 1 / 2 where is the Bayesian evidence for model , to compare these models. Generally, > 10 is taken to indicate a strong preference, based on the available data, for one model over another, and > 100 to indicate a decisive preference. Compared to the base model, the model with enhanced surface wind ejecta is preferred with > 10 8 .
However, our fits overwhelmingly prefer models that include shock cooling compared to either of the other models, with Bayes factors > 10 10 compared to the surface ejecta model and > 10 19 . Posterior distributions of free parameters for the preferred agnostic model (including shock cooling and enhanced surface wind ejecta) fit to GW170817, using multi-messenger priors (Table 1).
compared to the base model. The agnostic model (with both surface winds and shocks) is weakly favoured over the shock-only model, with = 2.2. The decisive preference for shock cooling can be understood from Figure 4. Models without shocks, generating their early luminosity only from radioactive decay in the blue ejecta, are too faint by ≈ 0.7 magnitudes (almost a factor 2 in luminosity) at ∼ 0.5 days. Increasing the blue ejecta mass through the parameter improves the fit at ∼ 1.5 days, but is still too faint during the first night (see appendix). However, a cocoon with an opening angle c ∼ 20 − 30 • provides the excess luminosity at = 0.5 days required to match the data. This supports previous work that has argued for the importance of shock cooling in the early light-curve (Kasliwal et al. 2017;Gottlieb et al. 2018;Piro & Kollmeier 2018).
As an important caveat, we note that there are other models that may be degenerate with cocoon emission. Free neutron decays, if present, could heat the outermost ejecta during the first few hours (Kulkarni 2005;Metzger et al. 2015), and produce a signal degenerate with shock heating . Nativi et al. (2021); Klion et al. (2021) recently showed that the bright early light-curve can alternatively be explained by accounting for an asymmetric ejecta distribution induced by the GRB jet. Finally, any differences in the nuclear heating rates from those assumed here would also affect the peak luminosity, as recently shown by Zhu et al. (2021). These effects are not modeled within our Bayesian analysis, and therefore remain a potential alternative to the shock-cooling interpretation.

Fit parameters and the origin of each luminosity component
Taking the agnostic version as our preferred model, we show the two-dimensional posteriors for this fit in Figure 5. As with all model variants considered, the chirp mass posterior is essentially the prior, since this parameter is constrained to such high precision by the GW data. The electromagnetic data tighten the mass ratio of the system to > 0.84 (about 50% of the GW-inferred prior volume) due to the requirement for blue ejecta. Viewing angles close to ∼ 30 • are preferred, more in line with afterglow modelling than with VLBI (Hajela et al. 2019;Nakar & Piran 2020). This is also consistent with the GW-inferred viewing angle at the distance of the GW170817 host .
The largest degeneracies in the model parameters are between the mass ratio and the surface ejecta enhancement 1/ , since equal mass binaries also produce more blue ejecta; between TOV and disk , since both affect the mass of purple disk ejecta; and between c and Λ, with a larger shock heating contribution to the luminosity relaxing the preference for low tidal deformability (compact EoS). We will discuss the EoS-dependent quantities, TOV and 1.4 (Λ), in more detail in section 3.3.
The best-fitting models produce typical ejecta masses blue ≈ 0.005 M , red ≈ 0.001 M , and purple ≈ 0.02 M . As we have established, the blue ejecta can only account for ∼ 50% of the flux at ∼ 0.5 days, and the rest is attributed to cooling of a GRB-shocked cocoon. The large mass of purple ejecta dominates the luminosity at ∼ 2 − 10 days, with the low mass of red ejecta contributing significantly only in the extended tail. These results are in agreement with Villar et al. (2017); however in our forward model we can uniquely associate the purple component with remnant disk winds, the blue component with shock-driven dynamical ejecta (possibly but not necessarily boosted by magnetic surface winds), and the red component with tidal dynamical ejecta, while identifying the importance of the cocoon contribution to the luminosity at early times.

Constraints on the NS equation of state
In this section, we examine the posteriors for EoS-dependent quantities in our fit, including estimates of the systematic uncertainty. While TOV is treated explicitly as a model parameter, the NS radii must be derived from the posterior of the symmetric tidal deformability, Λ s = 240 +94 −74 , and the masses of the two NSs, obtained from M and . The major source of systematic uncertainty in 1.4 and TOV is the scatter in relations 1 and 3, which are respectively ∼ 70% for the dynamical ejecta mass  and ∼ 50% of our inferred disk mass (Coughlin et al. 2019a). We test two methods to simultaneously model these systematic errors. First we run a suite of fits where equations 1 and 3 are modified each time by a random factor drawn from a Gaussian distribution with mean equal to 1 and width equal to the calibration uncertainties. Comparing the posteriors of 50 such runs, we find a scatter of ±40 in Λ s , and negligible scatter in TOV . As a complementary method, we run a single realisation of the fit but with two additional free parameters with Gaussian priors, dyn and disk (where the dynamical ejecta mass is dyn = dyn dyn and    Breschi et al. (2021). Top: NS radii from the posterior of Λ s . Middle: NS radii after re-weighting using equations of state that support an TOV within our 90% confidence interval. Bottom: construction of these prior weights using mass-radius curves from Dietrich et al. (2020). The red curves satisfy our constraint on TOV . Inset: Gaussian fit to allowed radii. the disk ejecta is treated analogously). In the former case, we broaden the posterior of Λ s using random draws from a Gaussian of width 40, whereas in the latter case we can use the posterior distribution of Λ s directly. We find that the two approaches yield indistinguishable posteriors for Λ s (and hence 1.4 ) and TOV .
Since the systematic error is negligible for TOV , our estimate of this quantity is equal to the posterior distribution shown in Figure 5. We measure a 90% credible interval TOV = 2.17 +0.08 −0.11 M . Models without shock cooling (strongly disfavoured by the Bayes factor analysis) prefer a more massive TOV ∼ 2.3 M . For comparison, the analysis of Margalit & Metzger (2017) found that TOV 2.2 M , which is consistent with our preferred shock cooling models but in potential tension with the shock-free models (although  have extended the Margalit & Metzger 2017 analysis and found TOV 2.3 M ). Other recent works Lucca & Sagunski 2020;Shao et al. 2020) have also favoured an EoS with TOV 2.2M . It is worth noting that the TOV constraints obtained here are complementary to these other works, which were based on arguments that are independent from the kilonova modeling adopted here. It is therefore interesting (and non-trivial) that the different approaches yield consistent results.
We next derive 1.4 from the posteriors of Λ s , M and , using equations 2, 6 and 7. During this conversion we include the systematic uncertainty introduced by the QURs (Yagi & Yunes 2016 by adding random scatter to the derived values for (6.5%) and Λ a (5%), drawn from Gaussian distributions. The masses of the primary and secondary, which enter equation 2, are 1 = 1.41 ± 0.04 M and 2 = 1.31 ± 0.03 M ; therefore their measured radii (especially for the primary) can be taken as very close to 1.4 .
The top panel of Figure 6 shows the derived probability density functions for the NS radii evaluated in this way. We find the primary has a radius 1 = 10.72 +1.67 −1.48 km (90% credible interval). The models without shock cooling prefer a smaller radius, ∼ 10 km, due to their lower values for Λ s . The most direct comparison for our radius is with the GW-inferred value; Abbott et al. (2018) find GW = 10.8 +2.0 −1.7 km (without imposing any EoS constraints a priori). Our measurement is in excellent agreement with this value, and has a 90% credible bound that is smaller by ∼ 0.5 km. Abbott et al. (2018) find a larger value for the radius and a tighter posterior distribution, GW = 11.9 ± 1.4 km, if they impose a constraint that the radius must be compatible with equations of state that support a TOV mass larger than the most massive known NS. In a similar vein, we can use our constraint on TOV to tighten our posteriors on the NS radii. This assumes that the posteriors of Λ and TOV are uncorrelated ( Figure 5 shows this to be a good approximation). We download the mass-radius curves for a sample of 4000 equations of state from Dietrich et al. (2020), constructed using chiral effective theory (up to 1.5 times nuclear saturation density, and a speed-of-sound parameterisation at higher densities) -an approach that is motivated by nuclear theory (Tews et al. 2018;Capano et al. 2020). As shown in the bottom panel of Figure 6, we select only those that support a TOV mass within our 90% confidence interval. The inset histogram shows the distribution of 1.4 for this EoS subset. We take this as the a priori probability of a given radius in order to re-weight our posteriors for 1 and 2 (after first transforming to a prior flat in 1.4 following Raithel et al. 2020), resulting in the middle panel of Figure 6. Our final measurement using Λ s and ensuring consistency with TOV is 1 = 11.06 +1.01 −0.98 km. Similar to Abbott et al. (2018), imposing a constraint on the EoS results in an increase in 1.4 and a tightening of the posterior. Our final value is consistent with theirs at around the 1 level, and narrower by ∼ 0.8 km. We also compare to multimessenger constraints from the literature that used the kilonova data for GW170817. Our median value for the radius lies just below the [11.3,13.5] km credible region obtained by Coughlin et al. (2019a), though there is significant overlap between their bounds and our 90% confidence interval. Our measurement is in excellent agreement with the 11.0 +0.9 −0.6 km found by Capano et al. (2020). It is also consistent with the 11.75 +0.86 −0.81 km found by Dietrich et al. (2020), at around the 1 level. Breschi et al. (2021) and Most et al. (2018) report slightly larger values for 1.4 , finding it to be 12.16 +0.89 −1.11 km or 12.39 +1.06 −0.39 km, respectively, which are consistent with our results only at the ∼ 2 level.

The importance of multi-messenger data
In this section we emphasise the importance of the multi-messenger constraints we used to fit GW170817. Table 1 lists two sets of priors: one set of suggested default priors for fitting an arbitrary NS merger, and one specifically tuned for GW170817. This latter set of priors, which we employed in the previous section, makes use of the GW constraints on M, , 0 andΛ, GRB afterglow constraints on (the GRB also constrains 0 ), and constraints on extinction from analysis of the host galaxy (Blanchard et al. 2017;Levan et al. 2017;Pan et al. 2017).
Taking our preferred model for GW170817 (the agnostic model including both shock cooling and an enhanced surface wind ejecta), we change the priors to the default set and re-run the fit to GW170817. Figure 7 shows the difference in posteriors when using the default priors, compared to the results using the multi-messenger priors. We find a fit of similar quality (ln = 79.2), but with much broader posterior distributions for several parameters.
The posterior distribution of M obtained using the broad priors has a median M broad = 1.2 M , close to the known value from the GW data. However, the credible region is ∼ 30 times wider than when using the very tight GW constraints. Another notable difference is in the viewing angle, where much larger viewing angles broad ∼ 45 • are preferred without the GRB constraints. The distribution of mass ratio spans a similar range as in the multi-messenger case, but is skewed towards more equal systems ( broad ∼ 0.94).
The posteriors of parameters that use the same priors in the default and multi-messenger models are (unsurprisingly) affected less severely, e.g. the posteriors of and cos c are very similar in the two cases, as are the posteriors ofΛ. However, the weaker constraints on the total system mass lead to a broadening of the posterior distribution of TOV , and a skew towards larger values.
The takeaway message is that our model can do a reasonable job of measuring important system properties for a GW170817-like merger even without a GW detection. However, having information from electromagnetic, GW, GRB and host galaxy analyses leads to tighter constraints, and this effect is perhaps most important when attempting to probe the NS EoS, as in section 3.3.

LIGHT CURVE PREDICTIONS: EXAMPLE OF GW190425
Having validated our model using GW170817, we now (briefly) demonstrate its utility for generating simulated kilonova data. While the broader applications in terms of population synthesis will be explored in future works, an important use case that we discuss here is to predict the electromagnetic observables of a particular kilonova given a detected GW signal (see e.g. , in order to make follow-up observations more efficient. This is especially relevant given the possibility of target-of-opportunity GW follow-up . Posterior distributions for model fit to GW170817 when using the broad priors (blue), compared to the posteriors obtained using the multi-messenger priors (black; same as Figure 5). Both sets of priors are listed in Table 1. with the upcoming Vera Rubin Observatory (Margutti et al. 2018a;Cowperthwaite et al. 2019;Smith et al. 2019;Chen et al. 2020).
During the LIGO-Virgo O3 run, the GW detectors identified several binary NS candidates. One event is contained in the latest GW source catalog, GWTC-2 (Abbott et al. 2020a), with the others discovered in the second half of O3 such that source parameters have not yet been released. The published event is GW190425, which was notable for its large total mass 1 + 2 = 3.4 M (Abbott et al. 2020c), implying that (unlike GW170817) at least one of the constituent NSs was more massive than the canonical ≈ 1.4 M .
Various groups of astronomers followed up the GW discovery of GW190425 with electromagnetic observations in an attempt to detect a counterpart, though without success on this occasion (Hosseinzadeh et al. 2019;Coughlin et al. 2019b;Lundquist et al. 2019;Antier et al. 2020). This may have been due in part to the large distance to the source, 159 +69 −71 Mpc, and in part because the signal was significant only in one GW detector, leading to a wide sky localization of > 8000 deg 2 (Abbott et al. 2020c). . The shaded regions correspond to the 90% credible ranges in magnitude for each band. Top: assuming a GW170817-like kilonova shifted to the GW-inferred distance of GW190425 (following Hosseinzadeh et al. 2019). Bottom: predicting the light curve directly from the GW constraints on the chirp mass, mass ratio, distance and inclination. The wider uncertainty bands and faster decline in this case show that assuming a GW170817-like event leads to unrealistically tight constraints and observations that may be too shallow to confidently detect this source.
pared to the light curves of GW170817 (shifted to the distance of GW190425). In Figure 8, we show these limits compared to our best-fit model of GW170817 in the left panel. In agreement with Hosseinzadeh et al. (2019), we see that the deep galaxy-targeted searches rule out a GW170817-like counterpart (in those particular galaxies) over most of the 90% plausible distance range. With our new model, we can use the chirp mass and mass ratio constraints obtained by LIGO and Virgo for GW190425 to predict a light curve more directly applicable to this source. A similar approach was adopted by Barbieri et al. (2020) for both NS-NS and NS-BH binaries. This prediction is not only tuned specifically to GW190425, but also folds in the uncertainties on the binary parameters as well as the distance uncertainty. We use a Gaussian distribution for the chirp mass M = 1.44 ± 0.013, and flat distributions 0.8 < < 1, 0 < cos < 1, 0.71 < cos c < 1, and a fiducial disk = 0.15 based on GW170817. The result is shown in the right panel of Figure 8.
Two features of this prediction are noteworthy. One is the increased width of the 90% distribution of credible light curves, allowing more space for a faint kilonova to go undetected at the observed depths. More significantly, the median light curve is fainter (by ∼ 0.7 mag in band at 1 day post-merger) and declines more rapidly than GW170817. The main reason for this is that the more massive remnant in this system is expected to collapse promptly to a BH, reducing the mass of the purple ejecta component by an order of magnitude compared to GW170817. With this faster decline, even many of the galaxy-targeted observations obtained 1 day after merger rule out only ∼ 50% of the plausible light curves. Our goal here is not to calculate precisely the fraction of the localization volume probed by optical observations of GW190425, but simply to make the general point that how constraining a detection limit is depends on the choice of comparison model, and that this in turn depends on the binary source parameters. A better strategy in future would be to tune the observed depths to a level necessary to cover the 90% credible range of kilonova luminosities for a given GW detection. At present, this is not possible, as the GW-inferred binary parameters needed to simulate the light curve are not released at the time of merger. We suggest that greater cooperation between the GW and electromagnetic astronomy communities could be very useful in this regard: if the GW source parameters (M, ,Λ) were released in low latency, we could use our model to immediately predict the required observing depth at any wavelength to find or rule out a kilonova counterpart at high significance.

CONCLUSIONS
We have presented a simple forward model that combines a range of literature results on ejecta masses, r-process opacities, and the NS EoS, to predict kilonova light curves directly from GW-accessible binary parameters. This model can also be used to derive limits on the kilonova rate from optical surveys, given a model for the BNS source population; this will be the subject of future work. Our code has been developed within the framework and is publicly available. We validated the model by fitting it to the well-observed kilonova associated with GW170817, incorporating GW and GRB afterglow results in our priors. We found that a model including a shock cooling cocoon component was overwhelmingly favoured by the Bayes factor, compared to models with r-process radioactive decay as the only luminosity source. Data obtained during the first day after merger was essential to differentiating these models (Arcavi 2018). Our approach directly associates each of the r-process ejecta components with a physical origin, where the low opacity material is shock-heated dynamical ejecta, and the higher opacity material that dominates the emission from ∼ 2 − 10 days after merger is a viscous disk wind. The very lanthanide-rich tidal ejecta makes a significant contribution only at very late times.
The posteriors of our model fit contain important information on the GW170817 progenitor system. We find that this system likely had a small mass ratio, ≈ 0.9, and was viewed ≈ 30 • off-axis. By including the maximum stable NS mass TOV and the symmetric tidal deformability Λ s in our free parameters, we constrained these EoS-dependent parameters. We measure TOV = 2.17 +0.08 −0.11 M (90% credible interval) and derive 1.4 =11.06 +1.01 −0.98 km using the posteriors of Λ, M, and TOV , carefully taking into account the systematic errors on all relations used. This radius constraint is consistent and competitive with others in the literature.
For high-significance detections in O3, the LIGO-Virgo collaboration provided publicly the distance and source classification probabilities (as NS, BH, NS-BH, or mass-gap binary) along with the sky localisation, giving astronomers a rough guide to sources they may wish to follow up electromagnetically. However, even for a source confidently classified as a NS merger, the binary parameters -especially the chirp mass, mass ratio and viewing angle -have a substantial impact on the likelihood of detecting the counterpart, since they control the peak luminosity, colour and decline rate ( Figures  2 and 3). If such parameters (which may be estimated even in low latency; Biscoveanu et al. 2019;Finstad & Brown 2020;Krastev et al. 2020) are made available to the astronomical community, this can aid real-time kilonova searches following binary NS merger discoveries from the GW detector network (also advocated for by Barbieri et al. 2020, for example). We demonstrated this by synthesising light curves appropriate to GW190425, for which optical searches did not uncover a counterpart, finding that many were not deep enough to fully rule out a kilonova. For future observing runs with a higher binary NS detection rate, the ability to better predict the required observing depth on a case-by-case basis would enable more efficient allocation of telescope time (which is of course a finite resource) for GW follow-up.
We have not discussed in this work the signature of a merger between a NS and a companion BH; these are also expected to produce kilonovae as long as the mass ratio is modest (though without shock-driven, low-lanthanide ejecta). Our next step will be to implement an analogous model to simulate light curves for this class of events (see also Barbieri et al. 2020).
Finding a larger sample of kilonovae in the coming years, as GW detectors and optical telescopes increase their survey power, will be crucial to improve our understanding of the NS EoS and the source population of NS binaries. Hopefully, the code we present here will be useful to the community both in helping to plan observations that increase this sample and in deriving important physics from the data obtained. Opacity (cm 2 g -1 ) 5000 K 7500 K 10000 K Figure A2. Opacity of r-process ejecta as a function of electron fraction and temperature from Tanaka et al. (2020). The broken polynomial fit is used to estimate the opacity of the disk ejecta, using the average predicted by Lippuner et al. (2017) for stable, short-lived or prompt-collapse merger remnants. In practice, the average > 0.25 for all reasonable input parameters. i-1 z-2 y-3 J-4 H-5 K-6 Figure B1. Left: Best-fit kilonova model for GW170817 with surface ejecta enhancement. Right: Best-fit model with shock cooling. The shock cooling model is preferred with a Bayes factor > 10 10 . This figure is similar to Figure 4 in the main text, which shows the best fit light curves including either both or neither of these effects.