Modelling Populations of Kilonovae

The 2017 detection of a kilonova coincident with gravitational-wave emission has identified neutron star mergers as the major source of the heaviest elements, and dramatically constrained alternative theories of gravity. Observing a population of such sources has the potential to transform cosmology, nuclear physics, and astrophysics. However, with only one confident multi-messenger detection currently available, modelling the diversity of signals expected from such a population requires improved theoretical understanding. In particular, models which are quick to evaluate, and are calibrated with more detailed multi-physics simulations, are needed to design observational strategies for kilonovae detection, and to obtain rapid-response interpretations of new observations. We use grey-opacity models to construct populations of kilonovae, spanning ejecta parameters predicted by numerical simulations. Our modelling focuses on wavelengths relevant for upcoming optical surveys, such as the Rubin Observatory Legacy Survey of Space and Time (LSST). In these simulations, we implement heating rates that are based on nuclear reaction network calculations. We create a Gaussian-process emulator for kilonova grey opacities, calibrated with detailed radiative transfer simulations. Using recent fits to numerical relativity simulations, we predict how the ejecta parameters from BNS mergers shape the population of kilonovae, accounting for the viewing-angle dependence. Our simulated population of binary neutron star (BNS) mergers produce peak i-band absolute magnitudes $-20 \leq M_i \leq -11$. A comparison with detailed radiative transfer calculations indicates that further improvements are needed to accurately reproduce spectral shapes over the full light curve evolution.


INTRODUCTION
Astrophysical understanding of the multi-messenger signals from binary neutron star (BNS) mergers has advanced considerably with the discovery of GW170817/AT2017gfo (Abbott et al. 2017a,b;The LIGO Scientific Collaboration et al. 2017;Coulter et al. 2017;Cowperthwaite et al. 2017;Kasliwal et al. 2017;Lipunov et al. 2017;Tanvir et al. 2017;Smartt et al. 2017;Soares-Santos et al. 2017;Valenti et al. 2017).Extensive analysis of this merger has converged on a description of the optical and infrared (OIR) electromagnetic (EM) signal as being produced by ejected mass with a composition that varies with the polar angle in the frame of the merger (Cowperthwaite et al. 2017;Perego et al. 2017;Kawaguchi et al. 2018).
The first detection intensified the efforts to realistically model ★ E-mail: christian.setzer@fysik.su.se neutron star mergers.While an entirely realistic modelling of all physical aspects is still out of reach, major efforts have been undertaken to explore the parameter space of BNS mergers in terms of masses, mass ratios, equations of state, spins, and other orbital parameters, albeit with approximate physics.With hundreds of available simulations, empirical relations between the intrinsic binary properties and the characteristic kilonova (kN/e) ejecta have been constructed (Kawaguchi et al. 2016;Dietrich & Ujevic 2017;Coughlin et al. 2019;Krüger & Foucart 2020;Nedora et al. 2022).Dietrich & Ujevic (2017), following the work by Foucart (2012); Kawaguchi et al. (2016) for black hole-neutron star systems, provided the first fits of empirical formulae relating the properties of a BNS merger, i.e., the component gravitational masses  1,2 and neutron star compactness parameters  1,2 , to the dynamical ejecta parameters: ejecta mass  ej and ejecta velocity  ej .The fitting formulae presented in Dietrich & Ujevic (2017) have since been updated in the works by Coughlin et al. (2019); Radice et al. (2018b); Krüger & Foucart (2020); Nedora et al. (2022).With these relations it has become feasible to investigate the diversity of kNe resulting from a set of priors describing the progenitor BNS population.
Understanding the population of kNe will help resolve systematics in measurements of the Hubble constant from standard sirens (Mortlock et al. 2019;Chen 2020;Coughlin et al. 2020;Moresco et al. 2022), improve population level inference of the nuclear matter equation of state (EOS) from neutron stars (Lackey & Wade 2015; Wysocki et al. 2020;Ghosh et al. 2022), and improve understanding of the rate of heavy element enrichment of the Universe (Cowan et al. 2021).As the population becomes more understood, observational selection biases can be studied and incorporated into analyses aimed at determining population level parameters, or place informed priors on multi-messenger parameter estimation from such events.A complete population description and resulting observational selection will be ever more important if standard sirens, being independent of distance-ladder calibration uncertainties, are going to be the arbiter of the Hubble-tension (Mortlock et al. 2019).
We present a population model of kNe where the parameters of the progenitor BNS describe the resulting kN transient signal distribution.We extend the kN signal modelling used in our previous work, Setzer et al. (2019), in several ways.First, we augment the relations presented by Coughlin et al. (2019) in their multimessenger parameter estimation of GW170817 with expressions from Radice et al. (2018b,c) for additional components of the ejecta to obtain a total ejecta mass contributing to the kN.Further, using simulation data from Radice et al. (2018b), we derive relations for the mass-averaged composition of the ejecta material as a function of viewing-angle.Then, using a library of detailed nuclear heating rates, we calibrate a semi-analytic grey-opacity kN model with high-fidelity radiation transport simulations; these steps constitute improved physical modelling with respect to recent population models presented by Nicholl et al. (2021) and Colombo et al. (2022).However, in contrast to these works, we simplify the computation of our model by neglecting composition differences from off-axis contributions of the ejecta.We then train a Gaussian process emulator to create draws from this calibrated model for a broad parameter space spanning the ejecta properties of the merger.Finally, we study the resulting distribution of kN light curves to get an understanding of the range and diversity of signals that are possible.This work represents a significant step towards improved modelling for rapidly simulating populations of kNe consistent with BNS merger progenitors and detailed physics, which is necessary to make accurate predictions for detection prospects (Rosswog et al. 2017;Scolnic et al. 2017;Setzer et al. 2019;Almualla et al. 2021;Mochkovitch et al. 2021;Sagués Carracedo et al. 2021;Andreoni et al. 2022;Chase et al. 2022;Colombo et al. 2022;Just et al. 2022).
In Sec. 2 we detail the choices of the priors on the population of neutron star binaries.Section 3 presents the updates to modelling of the optical-NIR kN signal resulting from BNS mergers.We present the results from these simulations in Sec. 4 and also discuss the dependence of the resulting kN signals on the binary parameters and the impact of modelling uncertainties on the resulting distributions.Finally, we conclude in Sec. 5 with a summary of the implications for future detections of BNS kNe and, discuss the additional astrophysical modelling needed to improve our understanding of these objects.We begin with the typical range considered for BNS mergers by the LVC searches, modifying it with additional restrictions.The feature in the lower-right corner arises from a limit on the maximum ejecta mass which correlates strongly with mass-ratio as seen in Fig. 2, rather than the minimum mass-ratio limit.This results in an effective mass-ratio limit of  0.55.

POPULATION OF KILONOVA-PRODUCING BNS
We will now define priors on the set of parameters which will allow us to simulate kNe for each BNS merger.While all mergers of neutron stars are predicted to produce some level of EM emission (Metzger & Berger 2012), not all will produce a kN signal that is bright enough to be observed by current and near-term instruments (Setzer et al. 2019;Sagués Carracedo et al. 2021).This is expected due to the magnitude limit in the context of surveys, and also due to the intrinsic variability in the population.To model the population of BNS mergers in the Universe, we use the broad prior on BNS masses used by the Laser Interferometer Gravitational-wave Observatory and Virgo Collaboration (LVC) template bank-based searches for compact binary coalescence signals (Canton & Harry 2017).For our population of BNS, this is a uniform prior on the component masses of the binary from one to three solar masses.However, from this starting point we make three modifications to the component mass prior.
We set the maximum neutron star mass to the maximum Tolman-Oppenheimer-Volkov (TOV) mass given by a choice of equation of state (Tolman 1939;Oppenheimer & Volkoff 1939),  TOV ≈ 2.0 − 2.5 M , see Godzieba et al. (2021); Drischler et al. (2021) and references therein.We also place a cut on the massratio,  =  min / max .We choose to simulate binaries in the range  ∈ [0.4,1].This cut is chosen to accommodate the range of massratios inferred from the detections of GW170817 and GW190425 (Abbott et al. 2019(Abbott et al. , 2020)), while also including the ranges from current observational limits of  min,obs ∼ 0.75 (Martinez et al. 2015) and the prediction of mass ratios from stellar population synthesis models of  min ∼ 0.5 − 0.6 (Dominik et al. 2013;Tauris et al. 2017;Kruckow et al. 2018;Andrews & Mandel 2019).Notably, given an imposed upper limit on the ejecta mass, this effectively limits mass ratios of our population to the range  ∈ [0.55, 1].The neutron star component masses are drawn uniformly from this joint distribution, see Fig. 1.We note that a portion of our population, where the total remnant mass is less than 1.2 TOV , will produce a long-lived neutron star remnant and likely have a magnetar-driven kilonovae (Margalit et al. 2022;Sarin et al. 2022).We do not modify our modelling for this sub-population, though we find approximately 50% of our population have a predicted total remnant mass below 1.2 TOV .
Additionally, we assume the binary systems are comprised of non-rotating neutron stars (Bildsten & Cutler 1992;Kochanek 1992).Further assuming that there is no precession of the orbital angular momentum of the binary, we can promote the binary mergerframe polar angle, i.e. the angle of total angular momentum, to the observer viewing-angle,  obs .Given this we can equate the observer viewing-angle and the binary orbital inclination.With the further assumption that there is symmetry of the merger ejecta about the binary orbital plane we map viewing-angles greater than 90 degrees as, (1) Each binary is then oriented assuming no preferential direction in alignment of their inclination.This is realized by drawing the observer viewing-angle/binary inclination angle from a uniform distribution on the sphere.Given we are assuming azimuthal symmetry, i.e., symmetry about the binary orbital axis, this maps to where  obs is the aziumthal angle of inclination on the unit sphere, which factors out due to axisymmetry.

Equation of State
Assuming all neutron stars obey a single equation of state (EOS), we consider those used by the simulations of Radice et al. (2018b).These EOS are consistent with the tidal deformability constraints from GW170817 (Radice et al. 2018a).Though none of the EOS are strongly favoured relative to another; based on the analysis of Coughlin et al. (2019) and due to its general use in the neutron star simulation community, we choose to use the SFHo EOS, with  TOV ≈ 2.06 M (Steiner et al. 2013).This EOS was designed with consideration towards properties of observed neutron stars and results of near-saturation density nuclear experiments.
The mass-radius relation, combined with the solution of the TOV equations for our given EOS, allows us to calculate the stellar compactness, , given by where  is the gravitational constant,  is the speed of light, and  is the neutron star radius.Given the range of variation allowable from the choice of EOS we will consider only EOS without a neutron star crust model.With this specification of EOS we now have completed our priors on the kN-progenitor binary systems.We will now describe how these parameters are mapped to the inputs necessary for simulating the counterpart kN signal.

MODELLING OF BNS KILONOVAE
We now describe the process of modelling each kN, beginning with a brief summary of the ejecta components which contribute to the kN-emission.To model a kN dependent on the BNS parameters described in Sec. 2, we map those parameters to those of the kN model.Our kN model characterizes the ejecta with four parameters: the total amount of ejected matter,  ej, total , the median outflow velocity,  ej , the electron fraction of the material,  e , and the grey opacity of the ejecta,  grey .The first three parameters are determined via empirical mappings from numerical simulations, discussed in Sec.3.2.In Sec.3.3, we describe the nuclear heating rate prescription we use and the radiation transport simulations used for calibrating grey-body opacities, which we compute with S N (Wollaeger et al. 2013(Wollaeger et al. , 2018;;Even et al. 2020).The process for determining the final parameter, the grey opacity, is described in 3.4.The process of calibrating the grey-opacities of the kN model with the additional radiation transport simulations is described in Sec. 3.4.1. In Sec. 3.4.3,we use these results to train a Gaussian process emulator to predict grey-opacities over the entire ejecta parameter space.

Phenomenological Description
As the inspiral of a binary neutron star system progresses, the two stars approach each other and will eventually come under the influence of the tidal forces of their companion beginning the process of ejecting material from the system; then, in their collision, squeezed and shock-heated material will also be ejected from the gravitationally bound system (Oechslin et al. 2007;Rosswog 2015;Tanaka 2016;Metzger 2020).These processes happen on the dynamical time-scale of the merger and form what is called the dynamical ejecta.Additional components of the ejecta occur from neutrino winds, magneto-rotational instabilities in the accretion disc, and secular processes that occur during the formation of the post-merger remnant, i.e. a super-massive or hyper-massive neutron star or a stellar mass black hole (Rosswog 2015;Radice et al. 2018b,b;Metzger 2020;Sarin & Lasky 2021).Assuming that the BNS inspiral is circular, and the binary is comprised of two non-spinning neutron stars, the resulting ejecta is expected to be symmetric about the merger plane.However, along the polar angle in the merger-frame, the properties of the ejecta will vary.Related back to the observer, this creates a viewing-angle dependence of kNe.
While each of these components of the ejecta have their own properties, and likely different composition profiles, we make the assumption that the light curves will be dominated by the contribution from the outermost ejecta in the observer's line-of-sight.This simplifying assumption allows for significant advantages in terms of decreased model complexity and computational efficiency.Further, Kawaguchi et al. (2020) indicate the outer dynamical ejecta do occlude, at least partially, contributions from the other components, suggesting their effect on the light curves is secondary.Given these assumptions we can then model the kN, including viewing-angle dependencies, by solving a 1D homologous flow radiation transport model with inputs specified by the total ejecta mass, median ejecta velocity, and the line-of-sight composition of the outermost material.

Connecting Ejecta Parameters to BNS Parameters
We now describe how we determine the bulk ejecta properties and line-of-sight composition, which are inputs into our model.Sampling the parameters { 1 ,  2 ,  obs } from the population prior, with a fixed EOS, we construct a mapping to all other parameters describing the simulated kN-signals { ej ,  ej ,  e ,  grey }.

Mapping Binary Parameters to Kilonova Ejecta Properties
Fully self-consistent simulations of BNS mergers from inspiral, General Relativistic (GR) merger, to GR-Magneto-Hydrodynamic (GR-MHD) outflow with fully relativistic radiation transport are immensely computationally expensive, even for a single event.We approach modelling of the BNS merger population using empirical relations derived from a substantial number of merger simulations which relate the binary parameters to the parameters of the ejecta producing the kN signal.We adopt the relations from Coughlin et al. (2019) which present updates to the fitting formula of Radice et al. (2018b) for  ej and  ej using 259 numerical simulations.They have found decreased error by fitting log 10 ( ej ) instead of  ej and have simplified the relations by removing the need for solving for the component baryonic masses (Coughlin et al. 2019).The equations and parameters are with  = −0.0719, = 0.2116,  = −2.42, = −2.905, = −0.3090, = −1.879,and  = 0.657, where [1 ↔ 2] refers to repetition of the preceding fitting terms with the binary parameter indices interchanged, i.e.  1 ↔  2 (Coughlin et al. 2019).This corresponds to a fractional error of 36% in log 10 ( fit ej ), equation 4 (dominated by the low-mass end,  ej ≈ 10 −4 − 10 −5 , though significantly better for larger  ej ) and 18% in  ej (Coughlin et al. 2019).
To obtain the ejected matter contributions to the total ejecta, we include additional relations for the secular ejecta coming from Radice et al. (2018b) and Coughlin et al. (2019).These primarily include an estimate for the remnant disc mass,  disc given by: where  tot is the total mass of the merging binary,  thresh , is the mass threshold for prompt black hole collapse, and  = −31.335, = −0.9760, = 1.0474,  = 0.05957 are the fitting parameters to numerical data.The black hole prompt collapse threshold mass is computed following Bauswein et al. (2013), which incorporates the chosen EOS through specification of the TOV maximum mass,  TOV , and the neutron star radius at 1.6 solar masses,  1.6 , This contributes to the total ejecta given an efficiency of unbinding of the disc (Radice et al. 2018b,c),  disc , due to secular processes, i.e., neutrino winds, GRMHD instabilities, etc.Given numerical results which show that the unbinding of the disc can contribute between 10% − 40% of the total disc mass to the ejecta, we sample uniformly between these two bounds to obtain the percentage of unbound material contributing to the ejecta for each kNe powering the EM transient (Metzger et al. 2008;Siegel & Metzger 2018;Miller et al. 2019;Fernández et al. 2019).Thus the total ejecta mass for each kN is given by, The ejecta parameters resulting from these mappings are shown in Fig. 2 and produce a complex shape with a significant concentration of kNe at high velocity and low ejecta mass.As noted by Radice et al. (2018b), these fitting formulae do not capture all the effects that contribute to the ejecta.Additionally, the numerical rel-ativity simulations which comprise the basis for these fits do not all simulate significantly beyond the time of merger.Further, the detailed microphysics of neutrino transport, magnetohydrodynamic (MHD) turbulence, and viscous effects are not usually simulated simultaneously.However, given the lack of data for BNS mergers and remaining uncertainty in the physics of these mergers, the above relations are sufficiently robust for simulating signals from a cosmological population of BNS mergers.Further investigating results from numerical simulations will allow us to encode additional dependencies of BNS kN ejecta on properties of the binary merger, such as spins and tidal deformability, which we leave to future work.

Viewing-Angle Dependence
Using the data from the simulations of Radice et al. (2018b), we adopt the description of the dynamical kN-ejecta properties as a function of the binary merger-frame polar angle  from Perego et al. (2017); see fig. 2 from their work.This uses the electron fraction,  e , as the primary property describing the composition of the material.We sub-select simulations from this dataset which do not promptly form a black hole (Agathos et al. 2020;Kölsch et al. 2021), i.e.,  BH ≤ 2 ms as this would not expected to be accompanied by an observable EM-counterpart due to extremely low total ejecta mass (Margalit & Metzger 2017).These simulation outputs are divided into angular slices which each contain a distribution of matter with different properties, e.g., fig. 4 of Radice et al. (2018b).
We take the mass-weighted average of the electron fraction data along each angular slice to obtain a profile of  e vs. .Noting the findings of Perego et al. (2017), a simple trigonometric function is able to describe the ejecta mass as a function of polar angle, i.e., m ej () ≈ sin 2 ().We obtain a similar description of the massweighted electron fraction,  e , by fitting trigonometric functions to the data.The best fit function for those tested, by least-squares optimisation, is for where  = 0.22704 and  = 0.16147 are the fit parameters for the SFHo EOS.We improve this fit further by splitting the simulation data based on the underlying equation of state and also selecting the subset of physics which includes not only neutrino cooling, but also neutrino heating.This reduces fit deviations substantially to max(Δ  e ) ≈ 0.05, see Fig. 3.We note that including additional parameters, i.e., the total binary mass and the binary mass ratio, at linear order did not show measurable improvement.The electron fraction composition is the primary determining factor of the nuclear heating rates and grey opacity which is described next.

Nuclear Heating
The nuclear heating rate, (), in the ejecta consists of the energy released during decays of a large number of radioactive isotopes which are produced in the rapid neutron-capture process (-process).
Although it has been demonstrated that such heating can be well approximated by the power law () ∼  − (Metzger et al. 2010;Hotokezaka et al. 2017), with  ≈ 1.2 − 1.3, the accuracy of this approximation is insufficient for our purposes.We therefore created a library of nuclear heating rates, parameterized by the initial electron fraction e ∈ [0.05, 0.5] and ejecta velocity  ej ∈ [0.05, 0.4].The nuclear heating rate was computed on a grid of these parameters using a nucleosynthesis network W NET (Winteler et al. 2012;Korobkin et al. 2012), an update of the B N network (Thielemann et al. 2011).This is the same nucleosynthesis code which was used to compute the -process nucleosynthesis in Wollaeger et al. (2018), using 5831 isotopes and the reaction rates from the compilation of Rauscher & Thielemann (2000).The nuclear masses far from stability are not experimentally known and one has to resort to theoretical mass models.In our network, we use the finite range droplet model (FRDM; Moller et al. 1995).The weak interaction rates are taken from Arcones & Martínez-Pinedo (2011).For fission and neutron capture, the fission rates of Panov et al. (2010) and -delayed fission probabilities of Panov et al. (2005) were used.
The nucleosynthesis calculations are performed along ejecta trajectories whose density as a function of time is computed from the ejecta mass  ej = 0.05  and expansion velocity  ej on the grid.The initial temperature is computed from the equation of state using the initial density and entropy, which is adopted to be  ej = 15   /baryon.Later, during the evolution of abundances, the entropy is incremented according to the produced heat, and the temperature is computed accordingly following Freiburghaus et al. (1999).
For each of the 120 points on the { e ,  ej }-grid, the nuclear heating rate is fit with an approximate formula, that has 11 fitting coefficients.The latter are then interpolated to obtain the values of fitting coefficients for the ejecta parameters in between the values of the grid.We describe the fitting procedure and the nuclear heating rates library elsewhere (Rosswog & Korobkin 2022).Both the full blown radiative transfer simulations with S N and our simpler semi-analytic model use the same nuclear heating rates from this library for the kN light curves calculations.
Instead of computing the efficiencies over the entire density profile of the ejecta, we replace the density profile (, r) with an averaged density.The details of our implementation are given in App. A.

Grey Opacity Dependence on Ejecta Properties
In order to produce a computationally inexpensive model useful for parameter estimation and rapid interpretation of observations, we employ a grey-opacity model for the kN emission.In general, the opacity of a kN has contributions from a large number of lines due to the presence of heavy elements (Kasen et al. 2013).The number of these lines can be greater than 10 6 for ions of some lanthanides and actinides, such as Terbium, Erbium, or Protactinium.However, we intend to summarize this with a single grey-opacity.As this is a significant simplification of the physics, we expect the grey opacity we infer to have dependencies on other parameters of the kN, such as the density and expansion-rate of the material.Thus we approach the mapping to grey-opacity including all parameters of the ejecta derived to this point, i.e., { ej ,  ej ,  e }.

Grey-body Model and Opacity Fitting Scheme
To accomplish this we compute synthetic spectra from multi-group radiative transfer simulations using S N (Wollaeger et al. 2013(Wollaeger et al. , 2018;;Even et al. 2020), for a range of kN-like ejecta trajectories.We choose a grid of kN-parameters, { ejecta ,  ejecta ,  e }, that spans the range of ejecta values predicted for the population of BNS kNe we are simulating to produce our training set.This represents a grid of the following parameter values: Θ train ej = { ej ∈ (0.002, 0.01, 0.03, 0.05, 0.08),  ej ∈ (0.05, 0.10, 0.15, 0.20),  e ∈ (0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40)}.However, to reduce the number simulations necessary to create the training set, we remove a subset of the grid that falls significantly outside the region of the population's parameters as shown in Fig. 2.This is also augmented with 14 additional simulations offset from the above grid, resulting in a total of 230 simulations in our training set.The results of these S N simulations provide where : f is the spectral flux density, σ2 is the variance of the spectral flux density,  is the index for the wavelength bins running from 1, ...,   and  indexes the time-steps of the simulation from 1, ...,   .Note, we also allow the number of wavelength bins to be dependent on the time-step, i.e.,  ,  .While simulations directly record the total flux per wavelength bin, F, this is converted to spectral flux density, f , for the output by dividing the flux by the width of the corresponding wavelength bin, i.e., F,  = ( bin,+1 −  bin, ) f,  .
We fit the S N spectral time-series with our grey-body model.We model the kN signal from each BNS merger using a semianalytic eigenmode expansion (SAEE) model presented by Wollaeger et al. (2018); Rosswog et al. (2018) and previously used in Setzer et al. (2019), see appendix A of Rosswog et al. (2018) for a comprehensive summary of the radiation transport physics.For reference, the spectral time-series is given by a blackbody with effective temperature that evolves according to eq.A.25 of Rosswog et al. (2018).We make several modifications to the underlying model which solves the diffusion equation and the spectral generation scheme.We introduce an additional parameter, the electron fraction, to model the viewing-angle dependence and to determine the heating rates, see Sec. 3.2.2 and Sec.3.3.This model, similarly to S N , is parameterized by the ejecta mass, ejecta velocity, and the electron fraction of the ejecta.However, it also contains an additional parameter, the grey opacity,  grey , i.e. the spectral flux density time-series produced by the model is a function  =  ,  ( grey , θ ejecta ). (15) We fit the grey opacity of the model to the spectral data from S N using a weighted chi-sq.method, see Fig. 4 for a representative example of this fit. 1

Restriction to Observationally Relevant Data
We are concerned with fitting the kN-data most accurately near peak luminosity, as this is when the transient would most likely be detected, enabling follow-up observations to be triggered.Additionally, we are concerned with emulating the signals as they would be observed by optical and near-infrared surveys.For these reasons, we make the following modifications to the data: • We remove S N data prior to 0.25 days, as this timeperiod is undergoing numerical relaxation from the initial conditions to a stable evolution.
• We remove S N data in the far infrared, keeping only wavelength bins with  < 11000 [ Å].
• To prioritize detectability of the kN, we weight the time period closest to peak luminosity higher with respect to the contributions to the total chi-sq.per model, see below.
• We remove data after 5 days, as this time-period is beyond the range when LTE radiation transport is reliable for kNe (Pognan et al. 2022).
The weighting scheme we adopt places a weight on a given timestep based on the relative luminosity of that time-step with respect to the peak luminosity of the S N model being fit.The weights are where the luminosity of each time-step of the S N data,   , is defined from the provided data as This weighting is implemented as a re-scaling of the errors that enter the chi-sq.calculation.Explicitly the re-scaled flux density errors, σ,  , are given by which then modifies the chi-sq.in the following manner, Having explored several functional forms for the weighting, we   ,  e = 0.15.This model represents a typical event drawn from the population.The upper-left panel shows the spectrum closest to peak luminosity of the grey-body model.The other epochs shown are typical of the time scales at which it is expected observational followup of kNe to be made.We see good agreement in the relative color, i.e. shape, of the spectrum, although the model tends to systematically under-predict the overall luminosity, i.e., the amplitude.The agreement of the bolometric luminosity, as seen in the lower-right panel, is good over the range of evolution most likely detectable, i.e., less than 2 days, with deviations increasing at later times.Note the grey-body fit uses the same ejecta parameters, i.e.,  ej ,  ej ,  e , as those used in the S N simulation and only  grey is solved for, as described in Sec.3.4.1.
adopt the above due to its simplicity and the fact that it appears to saturate how well the temperature evolution of the S N spectra can be fit with a simplified grey-body model.
For each set of spectra, we fix θ ejecta for the SAEE model to those from the simulation, and find the value of grey opacity which minimizes equation 19.Finding the corresponding  grey for all simulated points θ ejecta we define a sparse mapping to greyopacity, i.e.,  grey (θ ejecta ).This approximates a grey opacity surface spanning the ejecta parameters of interest.See Fig. 4 for an example fit from this procedure.

Gaussian Process Emulation
To extend this mapping to arbitrary parameter combinations we interpolate this surface using Gaussian processes with the package (Ambikasaran et al. 2016).We train on the set of 230 simulations represented by (x)-marks in Fig. 2. Two simulations were removed from this set due to clear fit failures and outlier chisq.values.Given the observed variation of the inferred values on the ejecta properties, we model the covariance of each dimension independently for any chosen kernel function.After exploring a handful of standard kernels, we have chosen to use the Matern 5/2 as the final implementation for our results (Genton 2001;Rasmussen & Williams 2006).This was chosen over the other kernels as it produced a solution which largely avoids the non-physical region of negative opacity values.
This kernel as implemented in , is given by (Ambikasaran et al. 2016): where Here  2 is the squared distance, given the metric , and x represent the input data coordinates with row, column indices , .The matrix elements of , and an overall amplitude of the kernel, are the hyperparameters that are optimized to minimize the log-loss of the Gaussian process with respect to the specified mean function given the training data above, Θ train ej .In our scenario we have set the off-diagonal terms to zero and optimized the diagonal terms independently.We construct a piece-wise mean function based on the results of Tanaka et al. (2020).
We test this interpolation scheme by performing a leave-oneout cross-validation test, see Fig. 5.This assesses the predictive quality of the emulator by training the emulator on the original training data leaving out one data point at a a time, predicting the value of that held-out datum, and iterating in this manner through the entire training set.In each iteration, the residual between the Leave-one-out cross-validation analysis of our chosen kernel, Matern 5/2.We show standardized residuals from a leave-one-out posterior predictive test of the emulator in comparison to a unit Gaussian distribution with zero mean.While there is a mild skew towards under-predicting the opacity, we note the large outliers are few and none are biased across the grey opacity threshold commonly used to approximate detectability, i.e.,  grey ≈ 10 cm 2 /g.
prediction at the location of the removed data point and the ground truth value from the training data, divided by the predicted emulator uncertainty at that point, is computed.This procedure yields a distribution of residuals, as shown in Fig. 5.For a perfect emulator, this distribution of residuals should match a Gaussian distribution with unit standard deviation and zero mean (shown for comparison in the figure).We find that the emulator mildly skews towards underpredicting the opacity, i.e., more negative residuals; however, this bias occurs in the extremely lanthanide-rich region where opacities inferred from the training data can reach more than 50 cm 2 g −1 .
The bias does not cause any of the predicted opacities to cross the approximate detectability threshold of 10 cm 2 g −1 (Setzer et al. 2019).Consequently, we do not expect this to significantly impact predictions for observations, though the presence of this bias does reflect the possibility for future improvements to the Gaussian process construction for this emulation.
In order to use this in forward-modelling of BNS merger kNe we use the trained Gaussian process to predict the opacity given the arbitrary ejecta parameter combinations of each simulated source.The Gaussian process predicts a mean value and an uncertainty; thus, we sample from this distribution to obtain an opacity for each simulated kN.To avoid unphysical opacities, we reject values below  grey = 0.1 cm 2 g −1 , the minimum opacity found in the analysis by Tanaka et al. (2020).We generally find the transitional behaviour of the grey opacity at an electron fraction of  e = 0.2 − 0.25, see Fig. 6.Studies of this mapping generally agree that there is a transition in the opacity, due to the change in elements abundances of the material, around an electron fraction of  e 0.25 (Korobkin et al. 2012;Kasen et al. 2013;Lippuner & Roberts 2015), also see fig. 2 of Rosswog et al. (2018).
We find that the range of opacity values predicted for any electron fraction is mildly dependent on the ejecta mass and velocity of the material, see Fig. 7.As the grey opacity is the only free parameter in the emulator we see that some of this dimming behav-  Figure 7. Slice through the ejecta parameter space of the Gaussian process emulator for the kN grey opacity.We show a slice through the region of ejecta parameter space predicted by the mappings of our population prior.This slice illustrates the dependence of the grey opacity on the median ejecta velocity of the material.This view of the emulator is shown at a fixed electron fraction, i.e.,  e = 0.2 near the transition region, stacking the solutions for each of the training data values of ejecta mass.
ior is reflected in the emulated grey opacity, such that lower total ejecta mass and higher median ejecta velocity lead to a larger grey opacity.This arises as our SAEE kN model does not simulate the same detailed physics as S N (such as the density-dependent thermalisation prescription).Thus, as we are fixing all parameters apart from the grey opacity, the physical dependencies get pushed into the variation of the grey opacity values we infer.
This completes the set of parameters necessary to simulate kNsignals with our SAEE model.Given the parameter mappings above, we can directly generate a grey-body spectral-timeseries given the source-frame neutron star component masses, observer viewingangle, and choice of EOS. 2

RESULTS & DISCUSSION
We have simulated a population of BNS kNe from realistic priors, incorporating EOS information, viewing-angle dependence, and relations between the intrinsic binary parameters and resulting kNe ejecta.This has allowed us to construct a population model of kNe light curves consistent with a progenitor BNS population.

Characteristics of the Simulated kN Population
Given the complex distribution of the sampled kN-parameters, as illustrated in Fig. 2, we expect a rich diversity of light curves predicted by this population.This is seen in the distribution of peak magnitudes, Fig. 8, and also the distribution of the duration when the light curve is within one magnitude of the peak magnitude, see Fig. 9.As models of kNe generally exhibit their maximum peak magnitudes in the redder optical/near-infrared wavelength range, we show these distributions in the Rubin Observatory's i-band.We find the peak brightness for the population varies between −20 ≤   ≤ −11.The distribution is approximately bimodal, with central peak-magnitudes of approximately −15.5 and −12 respectively.The brighter peak is correlated with lanthanide-free, i.e., high electron-fraction, low-opacity kNe and the dimmer peak with lanthanide-rich, i.e., low electron-fraction, high-opacity kNe.We also show in Fig. 8 that the peak magnitude 2 This model will be available as a standalone pip-install-able package, with documentation, available at: https://github.com/cnsetzer/Setzer2022_BNSpopkNe, and will be updated as improvements are made.  of GW170817/AT2017gfo is compatible with our kNe population.The distribution of time spent near peak brightness in this band peaks around 1.25 days with a minimum duration of approximately 0.25 days and a long tail to higher values.We note this tail is comprised of high-opacity kNe with slow-moving ejecta.Although the long time-scales would be optimistic for detection, these kNe are inherently the dimmest population, and modelling of this part of parameter space is more uncertain.We see in Fig. 10 that the population peak brightness strongly depends on the viewing-angle of the binary.There is at least a three magnitude difference from face-on to edge-on orientations.This is expected due to the deterministic relationship of the viewing-angle to the electron fraction composition and the strong relationship between electron fraction and grey opacity, see Fig. 6.This illustrates that the peak brightness is very sensitive to the composition.We have considered the case of the mass-weighted average composition profile, with contributions only of the line-of-sight material.It will be important to investigate in future work the impact of non-line-of-sight components of varying composition, aspherical morphologies, and the effect of differences in composition between the inner and outer regions of the ejecta.Studies which have already considered some of these effects, such as Bulla (2019), predict peak brightnesses within the range predicted by our population.
In the same figure, we also see a mild trend to greater peak brightness with increasing total mass.However, once the total mass approaches the threshold for prompt black hole collapse, there is a sharp decrease of about 1 magnitude in the peak brightness.In the middle panel we see some modelling artefacts related to the Gaussian process interpolation in the range around 50 degrees.It is very difficult to model the sharp change in opacity in this region without a large increase in training data.The artefact is also related to the piece-wise prescription of the mean function specified in our grey-opacity fitting scheme, see Sec. 3.4.3.
A significant modelling uncertainty is the contribution to the total ejecta mass coming from the amount of matter unbound from the remnant accretion disc.In Fig. 11, we show the variation in peak brightness due to the modelling uncertainty of the disc unbindingfraction.Over the range of unbinding percentages considered, we   see an approximate 0.5 mag change in the peak brightness.This clearly shows that the unbinding uncertainty is subdominant to the viewing-angle contribution to peak brightness.Indeed, we find the dependence of peak brightness with respect to all other parameters of the model is subdominant to the viewing-angle/composition contribution.

Comparison to GW170817/AT2017gfo
We find that predictions from our population model are able to reproduce features observed in the light curves of AT2017gfo, despite not calibrating the model to this event.Specifically, we see that the population is able to produce kNe with the same peak brightness and peak duration as GW170817/AT2017gfo.Considering Fig. 12, we see that the smooth spectrum of our grey-body model captures the overall spectral shape of the observations-based model of AT2017gfo (Scolnic et al. 2017).To make this comparison we did not fit the SAEE model to observations of AT2017gfo, but rather sampled the component masses and viewing-angle from the binary parameter posteriors for GW170817 (Abbott et al. 2019) to construct a posterior predictive test for AT2017gfo using our model.Given the GW170817 posterior, our model predicts a range of kNe from which we construct the quantiles shown in Fig. 12.For binary parameters consistent with the gravitational wave signal, our model produces kNe that show a 2-3 orders of magnitude variation in the amplitude of the flux and luminosity, which is due to the wide range of mass-ratios and inclinations supported by the GW170817 posterior.We find that the median spectra are in good agreement with observations-calibrated modelling of the kNe emission from AT2017gfo (Scolnic et al. 2017).
The median total luminosity over the observationally relevant region is also in close agreement with predictions of the DES GW170817 spectral model (Scolnic et al. 2017).The grey-body model captures the time-scales of the smooth rise and fall of the luminosity.These results demonstrate overall self-consistency of our model with the gravitational wave and electromagnetic data on this single event; however, please note that one should not expect to find detailed agreement with second-order effects in time-evolution with the single-component grey-body modelling used in this work.
Our results indicate overall that detection prospects for kNe populations derived assuming that GW170817/AT2017gfo is typical are likely too optimistic.While our population can accommodate the features of GW170817, it is brighter than the typical kN in our population.However, with only one BNS kN confirmed with multi-messenger observations, we stress that the uncertainties in modelling make it difficult to draw definitive conclusions on the prospects of future detections.It is clear that observing the next kNe will be crucial to inform modelling the broader population of kNe.(Scolnic et al. 2017) to the range of kNe produced with the SAEE model for 5000 draws of the component masses and viewing angle from the LVC posterior for GW170817 (Abbott et al. 2019).These binary parameters are mapped to the associated kN ejecta parameters through the relations specified in Sec. 3. Though our kN model does not include spin effects, we draw samples from the LVC high-spin posterior due to its greater support for higher mass ratios compatible with the full range we consider when simulating the population of kNe produced in this work.

CONCLUSIONS
We have modelled a population of kNe by relating the parameters of neutron star binaries to the parameters of kN ejecta, enabling simulation of light curves dependent only on the properties of the binary.This utilizes a series of mappings constructed from numerical simulations of more detailed physics.Specifically, we have constructed an analytic description of the composition profile for the outermost ejecta.Using detailed radiation transport simulations we have calibrated a grey opacity model and trained a Gaussian process emulator to predict grey-opacities over a broad parameter space of kN ejecta.
We find that the resulting population of BNS merger kNe produce a diverse range of kN-signals.Considering the most observationally relevant regions of the kN evolution, we predict a range of peak brightness, measured as the peak absolute magnitude in the Rubin observatory i-band, to be in the range of −20 ≤   ≤ −11.We additionally find that kNe are indeed short-lived transients, with the duration of their light curves around peak brightness lasting ∼ 1.25 days.Our results are consistent with the cosmological kNe population model presented by Colombo et al. (2022).Both analyses predict a short-lived population of kNe, spanning at least six magnitudes in peak brightness, and peak-times lasting a few days at most, see fig. 2 and fig. 10 of Colombo et al. (2022).
Our population model is consistent with GW170817, reproducing key features of the observations.The model also predicts a bimodal light curve distribution due to the transition from high to low electron fraction ejecta, which supports the commonly used assumption of a binary opacity choice when simulating kN light curves.However, it is clear that grey-body modelling has limitations, particularly in reproducing the detailed spectral shape.Furthermore, as we have our calibrated the model to early-time, near-peak magnitudes of the light curves, there is significant uncertainty in the late-time predictions of the model.This work represents a major step toward developing a robust population model for kNe.Better linking of BNS properties to kNe properties, via empirical mappings calibrated with simulation data, will be key to building further on this work.The connection of the light curve to the intrinsic parameters of the BNS enables studies of the combined EM-GW observational prospects of such a population.Additionally, observational prospects for different binary population synthesis models or equations of state can be explored.Furthermore, light curve simulations dependent on the parameters of the binary makes possible a joint analysis of GW and kN data, where both signals are sampled from a self-consistent population prior.In future work, we will use this population model to derive the observational selection function for EM-GW events, which is a critical ingredient for cosmological inference with such populations.

Figure 1 .
Figure1.Parameter space of BNS component masses for which we simulate a population of kNe.We begin with the typical range considered for BNS mergers by the LVC searches, modifying it with additional restrictions.The feature in the lower-right corner arises from a limit on the maximum ejecta mass which correlates strongly with mass-ratio as seen in Fig.2, rather than the minimum mass-ratio limit.This results in an effective mass-ratio limit of  0.55.

Figure 2 .
Figure 2. Ejecta parameter space given by the prior on the binary masses, Fig. 1, mapping equations 4 -5, and the SFHo EOS.The points are shaded by the mass ratio of the neutron stars comprising the binary which maps to each ejecta parameter pairing.Over-plotted are the locations in the ejecta parameter space of the training data used to calibrate the grey opacity of the model to SN radiation transport simulations.While this is a complex shape given the non-linear mappings of the parameters, we can identify the swooping feature at the left of the figure to arise from merger scenarios where the remnant object undergoes prompt collapse to a black hole, thus capturing more material and preventing efficient formation of an accretion disc from which additional material can be ejected through MHD-turbulence and other secular processes.

Figure 3 .
Figure3.Fit for the mass-weighted average electron fraction,  e ( ) for the SFHo equation of state with simulation data coming from a subset ofRadice et al. (2018b) which have a more detailed prescription for the microphysics.Deviations from the fit are greatest around the merger plane, i.e., a viewing angle of /2, where the structure of the electron fraction distribution becomes more complex.The two sets of points in the figure represent simulations with different initial conditions for the constituent masses of the binary.The blue circles correspond an equal-mass neutron star binary of two 1.35  neutron stars and the orange triangles correspond to a 1.4  -1.2  binary neutron star system.

Figure 4 .
Figure 4. Snapshots of the SED timeseries and overall bolometric luminosity evolution for a single model from our S N training set data and the grey-body model fit.The ejecta parameter values for this simulation are  ej = 0.03 [ ],  ej = 0.25 [ ],  e = 0.15.This model represents a typical event drawn from the population.The upper-left panel shows the spectrum closest to peak luminosity of the grey-body model.The other epochs shown are typical of the time scales at which it is expected observational followup of kNe to be made.We see good agreement in the relative color, i.e. shape, of the spectrum, although the model tends to systematically under-predict the overall luminosity, i.e., the amplitude.The agreement of the bolometric luminosity, as seen in the lower-right panel, is good over the range of evolution most likely detectable, i.e., less than 2 days, with deviations increasing at later times.Note the grey-body fit uses the same ejecta parameters, i.e.,  ej ,  ej ,  e , as those used in the S N simulation and only  grey is solved for, as described in Sec.3.4.1.
Figure5.Leave-one-out cross-validation analysis of our chosen kernel, Matern 5/2.We show standardized residuals from a leave-one-out posterior predictive test of the emulator in comparison to a unit Gaussian distribution with zero mean.While there is a mild skew towards under-predicting the opacity, we note the large outliers are few and none are biased across the grey opacity threshold commonly used to approximate detectability, i.e.,  grey ≈ 10 cm 2 /g.

Figure 6 .
Figure 6.2D histogram of the population of BNS kNe viewed in the predicted  grey − Θ obs plane.The viewing-angle, and hence the composition, has a strong relationship with the grey opacity seen most clearly at viewingangles greater than 55 deg., i.e., when viewed closer to edge-on.This transition to higher opacities implies an upper-limit opening angle of the lanthanide-rich ejecta of this population of approximately 35 deg.about the merger plane.

Figure 8 .
Figure 8. Peak absolute magnitude in the LSST i-band for a population of BNS merger kNe.The population has two peaks in the distribution which correspond to kNe which are either lanthanide-free along the observer lineof-sight, i.e., low electron fraction and thus brighter transient, or lanthaniderich in the direction of the observer, i.e., high electron fraction leading to more opaque material and a dimmer transient.The vertical line indicates the peak i-band absolute magnitude of GW170817/AT2017gfo as inferred from observations, i.e.,   = −15.7 from Smartt et al. (2017).

Figure 9 .
Figure 9. Distribution of time spent within one magnitude of the peak brightness in the LSST i-band by our simulated population.The distribution peaks around 1.25 days, indicating that the majority of BNS spend around 0.75 − 3 days within one magnitude of the peak brightness.

Figure 10 .
Figure 10.Two-dimensional histograms showing the variation of Rubin Obs.i-band absolute magnitude against intrinsic parameters of the binaries we simulate.(Left) Variation w.r.t the total mass of the binary, the vertical line indicates the threshold mass for prompt black hole collapse.(Middle) The variation w.r.t. the observer viewing-angle.(Right) The variation w.r.t. the compactness of the neutron stars averaged over the two components of the binary.

Figure 11 .
Figure 11.Two-dimensional histogram of the variation in i-band peak brightness vs. the fraction of material unbound from the remnant disc.Considering the mean of the two populations coming from high/low opacities, we see the percentage of unbound disc matter contributes at most an approximate 0.5 magnitude change in the resulting peak brightness.

Figure 12 .
Figure 12.Similar to Fig. 4.This compares the DES spectral time series model of AT2017gfo(Scolnic et al. 2017) to the range of kNe produced with the SAEE model for 5000 draws of the component masses and viewing angle from the LVC posterior for GW170817(Abbott et al. 2019).These binary parameters are mapped to the associated kN ejecta parameters through the relations specified in Sec. 3. Though our kN model does not include spin effects, we draw samples from the LVC high-spin posterior due to its greater support for higher mass ratios compatible with the full range we consider when simulating the population of kNe produced in this work.

Table 1 .
Prior on the progenitor BNS joint component mass distribution.The starting point is the broad prior used by LVC in their searches for neutron star gravitational wave emission, approximately 1 − 3  .This prior is modified using additional criteria.

Table 2 .
Approximate regions of allowable ejecta parameter space.This is meant to parameterize, crudely, the "elephant"-shaped region of ejecta parameter combinations allowable from the prior and mappings used in this work shown in Fig.2.